src/regression.q

   1 /* PSPP - linear regression.
   2    Copyright (C) 2005 Free Software Foundation, Inc.
   3    Written by Jason H Stover <jason@sakla.net>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18    02110-1301, USA. */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_cdf.h>
  23 #include <gsl/gsl_vector.h>
  24 #include <gsl/gsl_matrix.h>
  25 #include "alloc.h"
  26 #include "case.h"
  27 #include "casefile.h"
  28 #include "cat.h"
  29 #include "command.h"
  30 #include "dictionary.h"
  31 #include "error.h"
  32 #include "file-handle.h"
  33 #include "gettext.h"
  34 #include "lexer.h"
  35 #include <linreg/pspp_linreg.h>
  36 #include "tab.h"
  37 #include "var.h"
  38 #include "vfm.h"
  39
  40 /* (headers) */
  41
  42
  43 /* (specification)
  44    "REGRESSION" (regression_):
  45    *variables=varlist;
  46    statistics[st_]=r,
  47    coeff,
  48    anova,
  49    outs,
  50    zpp,
  51    label,
  52    sha,
  53    ci,
  54    bcov,
  55    ses,
  56    xtx,
  57    collin,
  58    tol,
  59    selection,
  60    f,
  61    defaults,
  62    all;
  63    ^dependent=varlist;
  64    ^method=enter.
  65 */
  66 /* (declarations) */
  67 /* (functions) */
  68 static struct cmd_regression cmd;
  69
  70 /*
  71   Array holding the subscripts of the independent variables.
  72  */
  73 size_t *indep_vars;
  74
  75 /*
  76   Return value for the procedure.
  77  */
  78 int pspp_reg_rc = CMD_SUCCESS;
  79
  80 static void run_regression (const struct casefile *, void *);
  81 /*
  82    STATISTICS subcommand output functions.
  83  */
  84 static void reg_stats_r (pspp_linreg_cache *);
  85 static void reg_stats_coeff (pspp_linreg_cache *);
  86 static void reg_stats_anova (pspp_linreg_cache *);
  87 static void reg_stats_outs (pspp_linreg_cache *);
  88 static void reg_stats_zpp (pspp_linreg_cache *);
  89 static void reg_stats_label (pspp_linreg_cache *);
  90 static void reg_stats_sha (pspp_linreg_cache *);
  91 static void reg_stats_ci (pspp_linreg_cache *);
  92 static void reg_stats_f (pspp_linreg_cache *);
  93 static void reg_stats_bcov (pspp_linreg_cache *);
  94 static void reg_stats_ses (pspp_linreg_cache *);
  95 static void reg_stats_xtx (pspp_linreg_cache *);
  96 static void reg_stats_collin (pspp_linreg_cache *);
  97 static void reg_stats_tol (pspp_linreg_cache *);
  98 static void reg_stats_selection (pspp_linreg_cache *);
  99 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
 100                                        int, pspp_linreg_cache *);
 101
 102 static void
 103 reg_stats_r (pspp_linreg_cache * c)
 104 {
 105   struct tab_table *t;
 106   int n_rows = 2;
 107   int n_cols = 5;
 108   double rsq;
 109   double adjrsq;
 110   double std_error;
 111
 112   assert (c != NULL);
 113   rsq = c->ssm / c->sst;
 114   adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
 115   std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
 116   t = tab_create (n_cols, n_rows, 0);
 117   tab_dim (t, tab_natural_dimensions);
 118   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 119   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 120   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 121   tab_vline (t, TAL_0, 1, 0, 0);
 122
 123   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
 124   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
 125   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
 126   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
 127   tab_float (t, 1, 1, TAB_RIGHT, sqrt (rsq), 10, 2);
 128   tab_float (t, 2, 1, TAB_RIGHT, rsq, 10, 2);
 129   tab_float (t, 3, 1, TAB_RIGHT, adjrsq, 10, 2);
 130   tab_float (t, 4, 1, TAB_RIGHT, std_error, 10, 2);
 131   tab_title (t, 0, _("Model Summary"));
 132   tab_submit (t);
 133 }
 134
 135 /*
 136   Table showing estimated regression coefficients.
 137  */
 138 static void
 139 reg_stats_coeff (pspp_linreg_cache * c)
 140 {
 141   size_t i;
 142   size_t j;
 143   int n_cols = 7;
 144   int n_rows;
 145   double t_stat;
 146   double pval;
 147   double coeff;
 148   double std_err;
 149   double beta;
 150   const char *label;
 151   struct tab_table *t;
 152
 153   assert (c != NULL);
 154   n_rows = c->n_coeffs + 2;
 155
 156   t = tab_create (n_cols, n_rows, 0);
 157   tab_headers (t, 2, 0, 1, 0);
 158   tab_dim (t, tab_natural_dimensions);
 159   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 160   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 161   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 162   tab_vline (t, TAL_0, 1, 0, 0);
 163
 164   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
 165   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 166   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
 167   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
 168   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 169   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
 170   coeff = c->coeff[0].estimate;
 171   tab_float (t, 2, 1, 0, coeff, 10, 2);
 172   std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
 173   tab_float (t, 3, 1, 0, std_err, 10, 2);
 174   beta = coeff / c->depvar_std;
 175   tab_float (t, 4, 1, 0, beta, 10, 2);
 176   t_stat = coeff / std_err;
 177   tab_float (t, 5, 1, 0, t_stat, 10, 2);
 178   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 179   tab_float (t, 6, 1, 0, pval, 10, 2);
 180   for (j = 1; j <= c->n_indeps; j++)
 181     {
 182       i = indep_vars[j];
 183       label = var_to_string (c->coeff[j].v);
 184       tab_text (t, 1, j + 1, TAB_CENTER, label);
 185       /*
 186          Regression coefficients.
 187        */
 188       coeff = c->coeff[j].estimate;
 189       tab_float (t, 2, j + 1, 0, coeff, 10, 2);
 190       /*
 191          Standard error of the coefficients.
 192        */
 193       std_err = sqrt (gsl_matrix_get (c->cov, j, j));
 194       tab_float (t, 3, j + 1, 0, std_err, 10, 2);
 195       /*
 196          'Standardized' coefficient, i.e., regression coefficient
 197          if all variables had unit variance.
 198        */
 199       beta = gsl_vector_get (c->indep_std, j);
 200       beta *= coeff / c->depvar_std;
 201       tab_float (t, 4, j + 1, 0, beta, 10, 2);
 202
 203       /*
 204          Test statistic for H0: coefficient is 0.
 205        */
 206       t_stat = coeff / std_err;
 207       tab_float (t, 5, j + 1, 0, t_stat, 10, 2);
 208       /*
 209          P values for the test statistic above.
 210        */
 211       pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 212       tab_float (t, 6, j + 1, 0, pval, 10, 2);
 213     }
 214   tab_title (t, 0, _("Coefficients"));
 215   tab_submit (t);
 216 }
 217
 218 /*
 219   Display the ANOVA table.
 220  */
 221 static void
 222 reg_stats_anova (pspp_linreg_cache * c)
 223 {
 224   int n_cols = 7;
 225   int n_rows = 4;
 226   const double msm = c->ssm / c->dfm;
 227   const double mse = c->sse / c->dfe;
 228   const double F = msm / mse;
 229   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
 230
 231   struct tab_table *t;
 232
 233   assert (c != NULL);
 234   t = tab_create (n_cols, n_rows, 0);
 235   tab_headers (t, 2, 0, 1, 0);
 236   tab_dim (t, tab_natural_dimensions);
 237
 238   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 239
 240   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 241   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 242   tab_vline (t, TAL_0, 1, 0, 0);
 243
 244   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 245   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 246   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 247   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 248   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 249
 250   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
 251   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
 252   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
 253
 254   /* Sums of Squares */
 255   tab_float (t, 2, 1, 0, c->ssm, 10, 2);
 256   tab_float (t, 2, 3, 0, c->sst, 10, 2);
 257   tab_float (t, 2, 2, 0, c->sse, 10, 2);
 258
 259
 260   /* Degrees of freedom */
 261   tab_float (t, 3, 1, 0, c->dfm, 4, 0);
 262   tab_float (t, 3, 2, 0, c->dfe, 4, 0);
 263   tab_float (t, 3, 3, 0, c->dft, 4, 0);
 264
 265   /* Mean Squares */
 266
 267   tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
 268   tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
 269
 270   tab_float (t, 5, 1, 0, F, 8, 3);
 271
 272   tab_float (t, 6, 1, 0, pval, 8, 3);
 273
 274   tab_title (t, 0, _("ANOVA"));
 275   tab_submit (t);
 276 }
 277 static void
 278 reg_stats_outs (pspp_linreg_cache * c)
 279 {
 280   assert (c != NULL);
 281 }
 282 static void
 283 reg_stats_zpp (pspp_linreg_cache * c)
 284 {
 285   assert (c != NULL);
 286 }
 287 static void
 288 reg_stats_label (pspp_linreg_cache * c)
 289 {
 290   assert (c != NULL);
 291 }
 292 static void
 293 reg_stats_sha (pspp_linreg_cache * c)
 294 {
 295   assert (c != NULL);
 296 }
 297 static void
 298 reg_stats_ci (pspp_linreg_cache * c)
 299 {
 300   assert (c != NULL);
 301 }
 302 static void
 303 reg_stats_f (pspp_linreg_cache * c)
 304 {
 305   assert (c != NULL);
 306 }
 307 static void
 308 reg_stats_bcov (pspp_linreg_cache * c)
 309 {
 310   int n_cols;
 311   int n_rows;
 312   int i;
 313   int j;
 314   int k;
 315   int row;
 316   int col;
 317   const char *label;
 318   struct tab_table *t;
 319
 320   assert (c != NULL);
 321   n_cols = c->n_indeps + 1 + 2;
 322   n_rows = 2 * (c->n_indeps + 1);
 323   t = tab_create (n_cols, n_rows, 0);
 324   tab_headers (t, 2, 0, 1, 0);
 325   tab_dim (t, tab_natural_dimensions);
 326   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 327   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 328   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 329   tab_vline (t, TAL_0, 1, 0, 0);
 330   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 331   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 332   for (i = 1; i < c->n_indeps + 1; i++)
 333     {
 334       j = indep_vars[(i - 1)];
 335       struct variable *v = cmd.v_variables[j];
 336       label = var_to_string (v);
 337       tab_text (t, 2, i, TAB_CENTER, label);
 338       tab_text (t, i + 2, 0, TAB_CENTER, label);
 339       for (k = 1; k < c->n_indeps + 1; k++)
 340         {
 341           col = (i <= k) ? k : i;
 342           row = (i <= k) ? i : k;
 343           tab_float (t, k + 2, i, TAB_CENTER,
 344                      gsl_matrix_get (c->cov, row, col), 8, 3);
 345         }
 346     }
 347   tab_title (t, 0, _("Coefficient Correlations"));
 348   tab_submit (t);
 349 }
 350 static void
 351 reg_stats_ses (pspp_linreg_cache * c)
 352 {
 353   assert (c != NULL);
 354 }
 355 static void
 356 reg_stats_xtx (pspp_linreg_cache * c)
 357 {
 358   assert (c != NULL);
 359 }
 360 static void
 361 reg_stats_collin (pspp_linreg_cache * c)
 362 {
 363   assert (c != NULL);
 364 }
 365 static void
 366 reg_stats_tol (pspp_linreg_cache * c)
 367 {
 368   assert (c != NULL);
 369 }
 370 static void
 371 reg_stats_selection (pspp_linreg_cache * c)
 372 {
 373   assert (c != NULL);
 374 }
 375
 376 static void
 377 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 378                            int keyword, pspp_linreg_cache * c)
 379 {
 380   if (keyword)
 381     {
 382       (*function) (c);
 383     }
 384 }
 385
 386 static void
 387 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 388 {
 389   /*
 390      The order here must match the order in which the STATISTICS
 391      keywords appear in the specification section above.
 392    */
 393   enum
 394   { r,
 395     coeff,
 396     anova,
 397     outs,
 398     zpp,
 399     label,
 400     sha,
 401     ci,
 402     bcov,
 403     ses,
 404     xtx,
 405     collin,
 406     tol,
 407     selection,
 408     f,
 409     defaults,
 410     all
 411   };
 412   int i;
 413   int d = 1;
 414
 415   if (keywords[all])
 416     {
 417       /*
 418          Set everything but F.
 419        */
 420       for (i = 0; i < f; i++)
 421         {
 422           keywords[i] = 1;
 423         }
 424     }
 425   else
 426     {
 427       for (i = 0; i < all; i++)
 428         {
 429           if (keywords[i])
 430             {
 431               d = 0;
 432             }
 433         }
 434       /*
 435          Default output: ANOVA table, parameter estimates,
 436          and statistics for variables not entered into model,
 437          if appropriate.
 438        */
 439       if (keywords[defaults] | d)
 440         {
 441           keywords[anova] = 1;
 442           keywords[outs] = 1;
 443           keywords[coeff] = 1;
 444           keywords[r] = 1;
 445         }
 446     }
 447   statistics_keyword_output (reg_stats_r, keywords[r], c);
 448   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 449   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 450   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 451   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 452   statistics_keyword_output (reg_stats_label, keywords[label], c);
 453   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 454   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 455   statistics_keyword_output (reg_stats_f, keywords[f], c);
 456   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 457   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 458   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 459   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 460   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 461   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 462 }
 463
 464 int
 465 cmd_regression (void)
 466 {
 467   if (!parse_regression (&cmd))
 468     {
 469       return CMD_FAILURE;
 470     }
 471   multipass_procedure_with_splits (run_regression, &cmd);
 472
 473   return pspp_reg_rc;
 474 }
 475
 476 /*
 477   Is variable k one of the dependent variables?
 478  */
 479 static int
 480 is_depvar (size_t k)
 481 {
 482   size_t j = 0;
 483   for (j = 0; j < cmd.n_dependent; j++)
 484     {
 485       /*
 486          compare_var_names returns 0 if the variable
 487          names match.
 488        */
 489       if (!compare_var_names (cmd.v_dependent[j], cmd.v_variables[k], NULL))
 490         return 1;
 491     }
 492   return 0;
 493 }
 494
 495 static void
 496 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
 497 {
 498   size_t i;
 499   size_t k;
 500   size_t n_data = 0;
 501   size_t row;
 502   int n_indep;
 503   int j = 0;
 504   const union value *val;
 505   struct casereader *r;
 506   struct casereader *r2;
 507   struct ccase c;
 508   const struct variable *v;
 509   struct recoded_categorical_array *ca;
 510   struct recoded_categorical *rc;
 511   struct design_matrix *X;
 512   gsl_vector *Y;
 513   pspp_linreg_cache *lcache;
 514   pspp_linreg_opts lopts;
 515
 516   n_data = casefile_get_case_cnt (cf);
 517   n_indep = cmd.n_variables - cmd.n_dependent;
 518   indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
 519
 520   Y = gsl_vector_alloc (n_data);
 521   lopts.get_depvar_mean_std = 1;
 522   lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
 523
 524   lcache = pspp_linreg_cache_alloc (n_data, n_indep);
 525   lcache->indep_means = gsl_vector_alloc (n_indep);
 526   lcache->indep_std = gsl_vector_alloc (n_indep);
 527
 528   /*
 529      Read from the active file. The first pass encodes categorical
 530      variables.
 531    */
 532   ca = cr_recoded_cat_ar_create (cmd.n_variables, cmd.v_variables);
 533   for (r = casefile_get_reader (cf);
 534        casereader_read (r, &c); case_destroy (&c))
 535     {
 536       for (i = 0; i < ca->n_vars; i++)
 537         {
 538           v = (*(ca->a + i))->v;
 539           val = case_data (&c, v->fv);
 540           cr_value_update (*(ca->a + i), val);
 541         }
 542     }
 543   cr_create_value_matrices (ca);
 544   X =
 545     design_matrix_create (n_indep, (const struct variable **) cmd.v_variables,
 546                           ca, n_data);
 547
 548   /*
 549      The second pass creates the design matrix.
 550    */
 551   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
 552        case_destroy (&c))
 553     /* Iterate over the cases. */
 554     {
 555       k = 0;
 556       row = casereader_cnum (r2) - 1;
 557       for (i = 0; i < cmd.n_variables; ++i)     /* Iterate over the variables
 558                                                    for the current case.
 559                                                  */
 560         {
 561           v = cmd.v_variables[i];
 562           val = case_data (&c, v->fv);
 563           /*
 564              Independent/dependent variable separation. The
 565              'variables' subcommand specifies a varlist which contains
 566              both dependent and independent variables. The dependent
 567              variables are specified with the 'dependent'
 568              subcommand. We need to separate the two.
 569            */
 570           if (is_depvar (i))
 571             {
 572               if (v->type != NUMERIC)
 573                 {
 574                   msg (SE, gettext ("Dependent variable must be numeric."));
 575                   pspp_reg_rc = CMD_FAILURE;
 576                   return;
 577                 }
 578               lcache->depvar = (const struct var *) v;
 579               gsl_vector_set (Y, row, val->f);
 580             }
 581           else
 582             {
 583               if (v->type == ALPHA)
 584                 {
 585                   rc = cr_var_to_recoded_categorical (v, ca);
 586                   design_matrix_set_categorical (X, row, v, val, rc);
 587                 }
 588               else if (v->type == NUMERIC)
 589                 {
 590                   design_matrix_set_numeric (X, row, v, val);
 591                 }
 592
 593               indep_vars[k] = i;
 594               k++;
 595               lopts.get_indep_mean_std[i] = 1;
 596             }
 597         }
 598     }
 599   /*
 600      Now that we know the number of coefficients, allocate space
 601      and store pointers to the variables that correspond to the
 602      coefficients.
 603    */
 604   lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
 605   for (i = 0; i < X->m->size2; i++)
 606     {
 607       j = i + 1;                /* The first coeff is the intercept. */
 608       lcache->coeff[j].v =
 609         (const struct variable *) design_matrix_col_to_var (X, i);
 610       assert (lcache->coeff[j].v != NULL);
 611     }
 612   /*
 613      Find the least-squares estimates and other statistics.
 614    */
 615   pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
 616   subcommand_statistics (cmd.a_statistics, lcache);
 617   gsl_vector_free (Y);
 618   design_matrix_destroy (X);
 619   pspp_linreg_cache_free (lcache);
 620   free (lopts.get_indep_mean_std);
 621   free (indep_vars);
 622   casereader_destroy (r);
 623   return;
 624 }
 625
 626 /*
 627   Local Variables:
 628   mode: c
 629   End:
 630 */