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