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 = 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 = gsl_vector_get (c->param_estimates, 0);
 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 = 0; j < c->n_indeps; j++)
 181     {
 182       i = indep_vars[j];
 183       struct variable *v = cmd.v_variables[i];
 184       label = var_to_string (v);
 185       tab_text (t, 1, j + 2, TAB_CENTER, label);
 186       /*
 187          Regression coefficients.
 188        */
 189       coeff = gsl_vector_get (c->param_estimates, j + 1);
 190       tab_float (t, 2, j + 2, 0, coeff, 10, 2);
 191       /*
 192          Standard error of the coefficients.
 193        */
 194       std_err = sqrt (gsl_matrix_get (c->cov, j + 1, j + 1));
 195       tab_float (t, 3, j + 2, 0, std_err, 10, 2);
 196       /*
 197          'Standardized' coefficient, i.e., regression coefficient
 198          if all variables had unit variance.
 199        */
 200       beta = gsl_vector_get (c->indep_std, j + 1);
 201       beta *= coeff / c->depvar_std;
 202       tab_float (t, 4, j + 2, 0, beta, 10, 2);
 203
 204       /*
 205          Test statistic for H0: coefficient is 0.
 206        */
 207       t_stat = coeff / std_err;
 208       tab_float (t, 5, j + 2, 0, t_stat, 10, 2);
 209       /*
 210          P values for the test statistic above.
 211        */
 212       pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 213       tab_float (t, 6, j + 2, 0, pval, 10, 2);
 214     }
 215   tab_title (t, 0, _("Coefficients"));
 216   tab_submit (t);
 217 }
 218
 219 /*
 220   Display the ANOVA table.
 221  */
 222 static void
 223 reg_stats_anova (pspp_linreg_cache * c)
 224 {
 225   int n_cols = 7;
 226   int n_rows = 4;
 227   const double msm = c->ssm / c->dfm;
 228   const double mse = c->sse / c->dfe;
 229   const double F = msm / mse;
 230   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
 231
 232   struct tab_table *t;
 233
 234   assert (c != NULL);
 235   t = tab_create (n_cols, n_rows, 0);
 236   tab_headers (t, 2, 0, 1, 0);
 237   tab_dim (t, tab_natural_dimensions);
 238
 239   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 240
 241   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 242   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 243   tab_vline (t, TAL_0, 1, 0, 0);
 244
 245   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 246   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 247   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 248   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 249   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 250
 251   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
 252   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
 253   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
 254
 255   /* Sums of Squares */
 256   tab_float (t, 2, 1, 0, c->ssm, 10, 2);
 257   tab_float (t, 2, 3, 0, c->sst, 10, 2);
 258   tab_float (t, 2, 2, 0, c->sse, 10, 2);
 259
 260
 261   /* Degrees of freedom */
 262   tab_float (t, 3, 1, 0, c->dfm, 4, 0);
 263   tab_float (t, 3, 2, 0, c->dfe, 4, 0);
 264   tab_float (t, 3, 3, 0, c->dft, 4, 0);
 265
 266   /* Mean Squares */
 267
 268   tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
 269   tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
 270
 271   tab_float (t, 5, 1, 0, F, 8, 3);
 272
 273   tab_float (t, 6, 1, 0, pval, 8, 3);
 274
 275   tab_title (t, 0, _("ANOVA"));
 276   tab_submit (t);
 277 }
 278 static void
 279 reg_stats_outs (pspp_linreg_cache * c)
 280 {
 281   assert (c != NULL);
 282 }
 283 static void
 284 reg_stats_zpp (pspp_linreg_cache * c)
 285 {
 286   assert (c != NULL);
 287 }
 288 static void
 289 reg_stats_label (pspp_linreg_cache * c)
 290 {
 291   assert (c != NULL);
 292 }
 293 static void
 294 reg_stats_sha (pspp_linreg_cache * c)
 295 {
 296   assert (c != NULL);
 297 }
 298 static void
 299 reg_stats_ci (pspp_linreg_cache * c)
 300 {
 301   assert (c != NULL);
 302 }
 303 static void
 304 reg_stats_f (pspp_linreg_cache * c)
 305 {
 306   assert (c != NULL);
 307 }
 308 static void
 309 reg_stats_bcov (pspp_linreg_cache * c)
 310 {
 311   int n_cols;
 312   int n_rows;
 313   int i;
 314   int j;
 315   int k;
 316   int row;
 317   int col;
 318   const char *label;
 319   struct tab_table *t;
 320
 321   assert (c != NULL);
 322   n_cols = c->n_indeps + 1 + 2;
 323   n_rows = 2 * (c->n_indeps + 1);
 324   t = tab_create (n_cols, n_rows, 0);
 325   tab_headers (t, 2, 0, 1, 0);
 326   tab_dim (t, tab_natural_dimensions);
 327   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 328   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 329   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 330   tab_vline (t, TAL_0, 1, 0, 0);
 331   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 332   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 333   for (i = 1; i < c->n_indeps + 1; i++)
 334     {
 335       j = indep_vars[(i - 1)];
 336       struct variable *v = cmd.v_variables[j];
 337       label = var_to_string (v);
 338       tab_text (t, 2, i, TAB_CENTER, label);
 339       tab_text (t, i + 2, 0, TAB_CENTER, label);
 340       for (k = 1; k < c->n_indeps + 1; k++)
 341         {
 342           col = (i <= k) ? k : i;
 343           row = (i <= k) ? i : k;
 344           tab_float (t, k + 2, i, TAB_CENTER,
 345                      gsl_matrix_get (c->cov, row, col), 8, 3);
 346         }
 347     }
 348   tab_title (t, 0, _("Coefficient Correlations"));
 349   tab_submit (t);
 350 }
 351 static void
 352 reg_stats_ses (pspp_linreg_cache * c)
 353 {
 354   assert (c != NULL);
 355 }
 356 static void
 357 reg_stats_xtx (pspp_linreg_cache * c)
 358 {
 359   assert (c != NULL);
 360 }
 361 static void
 362 reg_stats_collin (pspp_linreg_cache * c)
 363 {
 364   assert (c != NULL);
 365 }
 366 static void
 367 reg_stats_tol (pspp_linreg_cache * c)
 368 {
 369   assert (c != NULL);
 370 }
 371 static void
 372 reg_stats_selection (pspp_linreg_cache * c)
 373 {
 374   assert (c != NULL);
 375 }
 376
 377 static void
 378 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 379                            int keyword, pspp_linreg_cache * c)
 380 {
 381   if (keyword)
 382     {
 383       (*function) (c);
 384     }
 385 }
 386
 387 static void
 388 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 389 {
 390   /*
 391      The order here must match the order in which the STATISTICS
 392      keywords appear in the specification section above.
 393    */
 394   enum
 395   { r,
 396     coeff,
 397     anova,
 398     outs,
 399     zpp,
 400     label,
 401     sha,
 402     ci,
 403     bcov,
 404     ses,
 405     xtx,
 406     collin,
 407     tol,
 408     selection,
 409     f,
 410     defaults,
 411     all
 412   };
 413   int i;
 414   int d = 1;
 415
 416   if (keywords[all])
 417     {
 418       /*
 419          Set everything but F.
 420        */
 421       for (i = 0; i < f; i++)
 422         {
 423           keywords[i] = 1;
 424         }
 425     }
 426   else
 427     {
 428       for (i = 0; i < all; i++)
 429         {
 430           if (keywords[i])
 431             {
 432               d = 0;
 433             }
 434         }
 435       /*
 436          Default output: ANOVA table, parameter estimates,
 437          and statistics for variables not entered into model,
 438          if appropriate.
 439        */
 440       if (keywords[defaults] | d)
 441         {
 442           keywords[anova] = 1;
 443           keywords[outs] = 1;
 444           keywords[coeff] = 1;
 445           keywords[r] = 1;
 446         }
 447     }
 448   statistics_keyword_output (reg_stats_r, keywords[r], c);
 449   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 450   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 451   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 452   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 453   statistics_keyword_output (reg_stats_label, keywords[label], c);
 454   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 455   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 456   statistics_keyword_output (reg_stats_f, keywords[f], c);
 457   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 458   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 459   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 460   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 461   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 462   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 463 }
 464
 465 int
 466 cmd_regression (void)
 467 {
 468   if (!parse_regression (&cmd))
 469     {
 470       return CMD_FAILURE;
 471     }
 472   multipass_procedure_with_splits (run_regression, &cmd);
 473
 474   return pspp_reg_rc;
 475 }
 476
 477 /*
 478   Is variable k one of the dependent variables?
 479  */
 480 static int
 481 is_depvar (size_t k)
 482 {
 483   size_t j = 0;
 484   for (j = 0; j < cmd.n_dependent; j++)
 485     {
 486       /*
 487          compare_var_names returns 0 if the variable
 488          names match.
 489        */
 490       if (!compare_var_names (cmd.v_dependent[j], cmd.v_variables[k], NULL))
 491         return 1;
 492     }
 493   return 0;
 494 }
 495
 496 static void
 497 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
 498 {
 499   size_t i;
 500   size_t k;
 501   size_t n_data = 0;
 502   size_t row;
 503   int n_indep;
 504   int j = 0;
 505   const union value *val;
 506   struct casereader *r;
 507   struct casereader *r2;
 508   struct ccase c;
 509   const struct variable *v;
 510   struct recoded_categorical_array *ca;
 511   struct recoded_categorical *rc;
 512   struct design_matrix *X;
 513   gsl_vector *Y;
 514   pspp_linreg_cache *lcache;
 515   pspp_linreg_opts lopts;
 516
 517   n_data = casefile_get_case_cnt (cf);
 518   n_indep = cmd.n_variables - cmd.n_dependent;
 519   indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
 520
 521   Y = gsl_vector_alloc (n_data);
 522   lopts.get_depvar_mean_std = 1;
 523   lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
 524
 525   lcache = pspp_linreg_cache_alloc (n_data, n_indep);
 526   lcache->indep_means = gsl_vector_alloc (n_indep);
 527   lcache->indep_std = gsl_vector_alloc (n_indep);
 528
 529   /*
 530      Read from the active file. The first pass encodes categorical
 531      variables.
 532    */
 533   ca = cr_recoded_cat_ar_create (cmd.n_variables, cmd.v_variables);
 534   for (r = casefile_get_reader (cf);
 535        casereader_read (r, &c); case_destroy (&c))
 536     {
 537       for (i = 0; i < ca->n_vars; i++)
 538         {
 539           v = (*(ca->a + i))->v;
 540           val = case_data (&c, v->fv);
 541           cr_value_update (*(ca->a + i), val);
 542         }
 543     }
 544   cr_create_value_matrices (ca);
 545   X =
 546     design_matrix_create (n_indep, (const struct variable **) cmd.v_variables,
 547                           ca, n_data);
 548
 549   /*
 550      The second pass creates the design matrix.
 551    */
 552   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
 553        case_destroy (&c))
 554     /* Iterate over the cases. */
 555     {
 556       k = 0;
 557       row = casereader_cnum (r2) - 1;
 558       for (i = 0; i < cmd.n_variables; ++i)     /* Iterate over the variables
 559                                                    for the current case.
 560                                                  */
 561         {
 562           v = cmd.v_variables[i];
 563           val = case_data (&c, v->fv);
 564           /*
 565              Independent/dependent variable separation. The
 566              'variables' subcommand specifies a varlist which contains
 567              both dependent and independent variables. The dependent
 568              variables are specified with the 'dependent'
 569              subcommand. We need to separate the two.
 570            */
 571           if (is_depvar (i))
 572             {
 573               if (v->type != NUMERIC)
 574                 {
 575                   msg (SE, gettext ("Dependent variable must be numeric."));
 576                   pspp_reg_rc = CMD_FAILURE;
 577                   return;
 578                 }
 579               lcache->depvar = (const struct var *) v;
 580               gsl_vector_set (Y, row, val->f);
 581             }
 582           else
 583             {
 584               if (v->type == ALPHA)
 585                 {
 586                   rc = cr_var_to_recoded_categorical (v, ca);
 587                   design_matrix_set_categorical (X, row, v, val, rc);
 588                 }
 589               else if (v->type == NUMERIC)
 590                 {
 591                   design_matrix_set_numeric (X, row, v, val);
 592                 }
 593
 594               indep_vars[k] = i;
 595               k++;
 596               lopts.get_indep_mean_std[i] = 1;
 597             }
 598         }
 599     }
 600   /*
 601      Now that we know the number of coefficients, allocate space
 602      and store pointers to the variables that correspond to the
 603      coefficients.
 604    */
 605   lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
 606   for (i = 0; i < X->m->size2; i++)
 607     {
 608       j = i + 1;                /* The first coeff is the intercept. */
 609       lcache->coeff[j].v =
 610         (const struct variable *) design_matrix_col_to_var (X, i);
 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 */