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