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