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 table *t;
  98   int n_rows = 2;
  99   int n_cols = 5;
 100   double rsq;
 101   double adjrsq;
 102   double stderr;
 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   stderr = sqrt (2 * (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, stderr, 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   assert (c != NULL);
 303 }
 304 static void
 305 reg_stats_ses (pspp_linreg_cache * c)
 306 {
 307   assert (c != NULL);
 308 }
 309 static void
 310 reg_stats_xtx (pspp_linreg_cache * c)
 311 {
 312   assert (c != NULL);
 313 }
 314 static void
 315 reg_stats_collin (pspp_linreg_cache * c)
 316 {
 317   assert (c != NULL);
 318 }
 319 static void
 320 reg_stats_tol (pspp_linreg_cache * c)
 321 {
 322   assert (c != NULL);
 323 }
 324 static void
 325 reg_stats_selection (pspp_linreg_cache * c)
 326 {
 327   assert (c != NULL);
 328 }
 329
 330 static void
 331 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 332                            int keyword, pspp_linreg_cache * c)
 333 {
 334   if (keyword)
 335     {
 336       (*function) (c);
 337     }
 338 }
 339
 340 static void
 341 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 342 {
 343   /*
 344      The order here must match the order in which the STATISTICS
 345      keywords appear in the specification section above.
 346    */
 347   enum
 348   { r,
 349     coeff,
 350     anova,
 351     outs,
 352     zpp,
 353     label,
 354     sha,
 355     ci,
 356     bcov,
 357     ses,
 358     xtx,
 359     collin,
 360     tol,
 361     selection,
 362     f,
 363     defaults,
 364     all
 365   };
 366   int i;
 367   int d = 1;
 368
 369   if (keywords[all])
 370     {
 371       /*
 372          Set everything but F.
 373        */
 374       for (i = 0; i < f; i++)
 375         {
 376           *(keywords + i) = 1;
 377         }
 378     }
 379   else
 380     {
 381       for (i = 0; i < all; i++)
 382         {
 383           if (keywords[i])
 384             {
 385               d = 0;
 386             }
 387         }
 388       /*
 389          Default output: ANOVA table, parameter estimates,
 390          and statistics for variables not entered into model,
 391          if appropriate.
 392        */
 393       if (keywords[defaults] | d)
 394         {
 395           *(keywords + anova) = 1;
 396           *(keywords + outs) = 1;
 397           *(keywords + coeff) = 1;
 398           *(keywords + r) = 1;
 399         }
 400     }
 401   statistics_keyword_output (reg_stats_r, keywords[r], c);
 402   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 403   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 404   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 405   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 406   statistics_keyword_output (reg_stats_label, keywords[label], c);
 407   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 408   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 409   statistics_keyword_output (reg_stats_f, keywords[f], c);
 410   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 411   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 412   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 413   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 414   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 415   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 416 }
 417
 418 int
 419 cmd_regression (void)
 420 {
 421   if (!parse_regression (&cmd))
 422     {
 423       return CMD_FAILURE;
 424     }
 425   multipass_procedure_with_splits (run_regression, &cmd);
 426
 427   return CMD_SUCCESS;
 428 }
 429
 430 /*
 431   Is variable k one of the dependent variables?
 432  */
 433 static int
 434 is_depvar (size_t k)
 435 {
 436   size_t j = 0;
 437   for (j = 0; j < cmd.n_dependent; j++)
 438     {
 439       /*
 440          compare_var_names returns 0 if the variable
 441          names match.
 442        */
 443       if (!compare_var_names (cmd.v_dependent[j], cmd.v_variables[k], NULL))
 444         return 1;
 445     }
 446   return 0;
 447 }
 448
 449 static void
 450 run_regression (const struct casefile *cf, void *cmd_)
 451 {
 452   size_t i;
 453   size_t k;
 454   size_t n_data = 0;
 455   size_t row;
 456   int n_indep;
 457   const union value *val;
 458   struct casereader *r;
 459   struct casereader *r2;
 460   struct ccase c;
 461   const struct variable *v;
 462   struct recoded_categorical_array *ca;
 463   struct recoded_categorical *rc;
 464   struct design_matrix *X;
 465   gsl_vector *Y;
 466   pspp_linreg_cache *lcache;
 467   pspp_linreg_opts lopts;
 468
 469   n_data = casefile_get_case_cnt (cf);
 470   n_indep = cmd.n_variables - cmd.n_dependent;
 471   indep_vars = (size_t *) malloc (n_indep * sizeof (*indep_vars));
 472
 473   Y = gsl_vector_alloc (n_data);
 474   lopts.get_depvar_mean_std = 1;
 475   lopts.get_indep_mean_std = (int *) malloc (n_indep * sizeof (int));
 476
 477   lcache = pspp_linreg_cache_alloc (n_data, n_indep);
 478   lcache->indep_means = gsl_vector_alloc (n_indep);
 479   lcache->indep_std = gsl_vector_alloc (n_indep);
 480
 481   /*
 482      Read from the active file. The first pass encodes categorical
 483      variables.
 484    */
 485   ca = cr_recoded_cat_ar_create (cmd.n_variables, cmd.v_variables);
 486   for (r = casefile_get_reader (cf);
 487        casereader_read (r, &c); case_destroy (&c))
 488     {
 489       for (i = 0; i < ca->n_vars; i++)
 490         {
 491           v = (*(ca->a + i))->v;
 492           val = case_data (&c, v->fv);
 493           cr_value_update (*(ca->a + i), val);
 494         }
 495     }
 496   cr_create_value_matrices (ca);
 497   X =
 498     design_matrix_create (n_indep, (const struct variable **) cmd.v_variables,
 499                           ca, n_data);
 500
 501   /*
 502      The second pass creates the design matrix.
 503    */
 504   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
 505        case_destroy (&c))
 506     /* Iterate over the cases. */
 507     {
 508       k = 0;
 509       row = casereader_cnum (r2) - 1;
 510       for (i = 0; i < cmd.n_variables; ++i)     /* Iterate over the variables
 511                                                    for the current case.
 512                                                  */
 513         {
 514           v = cmd.v_variables[i];
 515           val = case_data (&c, v->fv);
 516           /*
 517              Independent/dependent variable separation. The
 518              'variables' subcommand specifies a varlist which contains
 519              both dependent and independent variables. The dependent
 520              variables are specified with the 'dependent'
 521              subcommand. We need to separate the two.
 522            */
 523           if (is_depvar (i))
 524             {
 525               assert(v->type == NUMERIC);
 526               gsl_vector_set (Y, row, val->f);
 527             }
 528           else
 529             {
 530               if (v->type == ALPHA)
 531                 {
 532                   rc = cr_var_to_recoded_categorical (v, ca);
 533                   design_matrix_set_categorical (X, row, v, val, rc);
 534                 }
 535               else if (v->type == NUMERIC)
 536                 {
 537                   design_matrix_set_numeric (X, row, v, val);
 538                 }
 539
 540               indep_vars[k] = i;
 541               k++;
 542               lopts.get_indep_mean_std[i] = 1;
 543             }
 544         }
 545     }
 546   /*
 547      Find the least-squares estimates and other statistics.
 548    */
 549   pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
 550   subcommand_statistics (cmd.a_statistics, lcache);
 551   gsl_vector_free (Y);
 552   design_matrix_destroy (X);
 553   pspp_linreg_cache_free (lcache);
 554   free (lopts.get_indep_mean_std);
 555   free (indep_vars);
 556   casereader_destroy (r);
 557 }
 558
 559 /*
 560   Local Variables:
 561   mode: c
 562   End:
 563 */