src/language/stats/regression.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005, 2009 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <gsl/gsl_cdf.h>
  20 #include <gsl/gsl_matrix.h>
  21 #include <gsl/gsl_vector.h>
  22 #include <math.h>
  23 #include <stdlib.h>
  24
  25 #include <data/case.h>
  26 #include <data/casegrouper.h>
  27 #include <data/casereader.h>
  28 #include <data/category.h>
  29 #include <data/dictionary.h>
  30 #include <data/missing-values.h>
  31 #include <data/procedure.h>
  32 #include <data/transformations.h>
  33 #include <data/value-labels.h>
  34 #include <data/variable.h>
  35 #include <language/command.h>
  36 #include <language/dictionary/split-file.h>
  37 #include <language/data-io/file-handle.h>
  38 #include <language/lexer/lexer.h>
  39 #include <libpspp/compiler.h>
  40 #include <libpspp/message.h>
  41 #include <libpspp/taint.h>
  42 #include <math/design-matrix.h>
  43 #include <math/coefficient.h>
  44 #include <math/linreg.h>
  45 #include <math/moments.h>
  46 #include <output/tab.h>
  47
  48 #include "xalloc.h"
  49
  50 #include "gettext.h"
  51 #define _(msgid) gettext (msgid)
  52
  53 #define REG_LARGE_DATA 1000
  54
  55 /* (headers) */
  56
  57 /* (specification)
  58    "REGRESSION" (regression_):
  59    *variables=custom;
  60    +statistics[st_]=r,
  61                     coeff,
  62                     anova,
  63                     outs,
  64                     zpp,
  65                     label,
  66                     sha,
  67                     ci,
  68                     bcov,
  69                     ses,
  70                     xtx,
  71                     collin,
  72                     tol,
  73                     selection,
  74                     f,
  75                     defaults,
  76                     all;
  77    ^dependent=varlist;
  78    +save[sv_]=resid,pred;
  79    +method=enter.
  80 */
  81 /* (declarations) */
  82 /* (functions) */
  83 static struct cmd_regression cmd;
  84
  85 /*
  86   Moments for each of the variables used.
  87  */
  88 struct moments_var
  89 {
  90   struct moments1 *m;
  91   const struct variable *v;
  92 };
  93
  94 /*
  95   Transformations for saving predicted values
  96   and residuals, etc.
  97  */
  98 struct reg_trns
  99 {
 100   int n_trns;                   /* Number of transformations. */
 101   int trns_id;                  /* Which trns is this one? */
 102   pspp_linreg_cache *c;         /* Linear model for this trns. */
 103 };
 104 /*
 105   Variables used (both explanatory and response).
 106  */
 107 static const struct variable **v_variables;
 108
 109 /*
 110   Number of variables.
 111  */
 112 static size_t n_variables;
 113
 114 static bool run_regression (struct casereader *, struct cmd_regression *,
 115                             struct dataset *, pspp_linreg_cache **);
 116
 117 /*
 118    STATISTICS subcommand output functions.
 119  */
 120 static void reg_stats_r (pspp_linreg_cache *);
 121 static void reg_stats_coeff (pspp_linreg_cache *);
 122 static void reg_stats_anova (pspp_linreg_cache *);
 123 static void reg_stats_outs (pspp_linreg_cache *);
 124 static void reg_stats_zpp (pspp_linreg_cache *);
 125 static void reg_stats_label (pspp_linreg_cache *);
 126 static void reg_stats_sha (pspp_linreg_cache *);
 127 static void reg_stats_ci (pspp_linreg_cache *);
 128 static void reg_stats_f (pspp_linreg_cache *);
 129 static void reg_stats_bcov (pspp_linreg_cache *);
 130 static void reg_stats_ses (pspp_linreg_cache *);
 131 static void reg_stats_xtx (pspp_linreg_cache *);
 132 static void reg_stats_collin (pspp_linreg_cache *);
 133 static void reg_stats_tol (pspp_linreg_cache *);
 134 static void reg_stats_selection (pspp_linreg_cache *);
 135 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
 136                                        int, pspp_linreg_cache *);
 137
 138 static void
 139 reg_stats_r (pspp_linreg_cache * c)
 140 {
 141   struct tab_table *t;
 142   int n_rows = 2;
 143   int n_cols = 5;
 144   double rsq;
 145   double adjrsq;
 146   double std_error;
 147
 148   assert (c != NULL);
 149   rsq = c->ssm / c->sst;
 150   adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
 151   std_error = sqrt (pspp_linreg_mse (c));
 152   t = tab_create (n_cols, n_rows);
 153   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 154   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 155   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 156   tab_vline (t, TAL_0, 1, 0, 0);
 157
 158   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
 159   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
 160   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
 161   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
 162   tab_double (t, 1, 1, TAB_RIGHT, sqrt (rsq), NULL);
 163   tab_double (t, 2, 1, TAB_RIGHT, rsq, NULL);
 164   tab_double (t, 3, 1, TAB_RIGHT, adjrsq, NULL);
 165   tab_double (t, 4, 1, TAB_RIGHT, std_error, NULL);
 166   tab_title (t, _("Model Summary"));
 167   tab_submit (t);
 168 }
 169
 170 /*
 171   Table showing estimated regression coefficients.
 172  */
 173 static void
 174 reg_stats_coeff (pspp_linreg_cache * c)
 175 {
 176   size_t j;
 177   int n_cols = 7;
 178   int n_rows;
 179   int this_row;
 180   double t_stat;
 181   double pval;
 182   double std_err;
 183   double beta;
 184   const char *label;
 185
 186   const struct variable *v;
 187   const union value *val;
 188   struct tab_table *t;
 189
 190   assert (c != NULL);
 191   n_rows = c->n_coeffs + 3;
 192
 193   t = tab_create (n_cols, n_rows);
 194   tab_headers (t, 2, 0, 1, 0);
 195   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 196   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 197   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 198   tab_vline (t, TAL_0, 1, 0, 0);
 199
 200   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
 201   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 202   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
 203   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
 204   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 205   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
 206   tab_double (t, 2, 1, 0, c->intercept, NULL);
 207   std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
 208   tab_double (t, 3, 1, 0, std_err, NULL);
 209   tab_double (t, 4, 1, 0, 0.0, NULL);
 210   t_stat = c->intercept / std_err;
 211   tab_double (t, 5, 1, 0, t_stat, NULL);
 212   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 213   tab_double (t, 6, 1, 0, pval, NULL);
 214   for (j = 0; j < c->n_coeffs; j++)
 215     {
 216       struct string tstr;
 217       ds_init_empty (&tstr);
 218       this_row = j + 2;
 219
 220       v = pspp_coeff_get_var (c->coeff[j], 0);
 221       label = var_to_string (v);
 222       /* Do not overwrite the variable's name. */
 223       ds_put_cstr (&tstr, label);
 224       if (var_is_alpha (v))
 225         {
 226           /*
 227              Append the value associated with this coefficient.
 228              This makes sense only if we us the usual binary encoding
 229              for that value.
 230            */
 231
 232           val = pspp_coeff_get_value (c->coeff[j], v);
 233
 234           var_append_value_name (v, val, &tstr);
 235         }
 236
 237       tab_text (t, 1, this_row, TAB_CENTER, ds_cstr (&tstr));
 238       /*
 239          Regression coefficients.
 240        */
 241       tab_double (t, 2, this_row, 0, c->coeff[j]->estimate, NULL);
 242       /*
 243          Standard error of the coefficients.
 244        */
 245       std_err = sqrt (gsl_matrix_get (c->cov, j + 1, j + 1));
 246       tab_double (t, 3, this_row, 0, std_err, NULL);
 247       /*
 248          Standardized coefficient, i.e., regression coefficient
 249          if all variables had unit variance.
 250        */
 251       beta = pspp_coeff_get_sd (c->coeff[j]);
 252       beta *= c->coeff[j]->estimate / c->depvar_std;
 253       tab_double (t, 4, this_row, 0, beta, NULL);
 254
 255       /*
 256          Test statistic for H0: coefficient is 0.
 257        */
 258       t_stat = c->coeff[j]->estimate / std_err;
 259       tab_double (t, 5, this_row, 0, t_stat, NULL);
 260       /*
 261          P values for the test statistic above.
 262        */
 263       pval =
 264         2 * gsl_cdf_tdist_Q (fabs (t_stat),
 265                              (double) (c->n_obs - c->n_coeffs));
 266       tab_double (t, 6, this_row, 0, pval, NULL);
 267       ds_destroy (&tstr);
 268     }
 269   tab_title (t, _("Coefficients"));
 270   tab_submit (t);
 271 }
 272
 273 /*
 274   Display the ANOVA table.
 275  */
 276 static void
 277 reg_stats_anova (pspp_linreg_cache * c)
 278 {
 279   int n_cols = 7;
 280   int n_rows = 4;
 281   const double msm = c->ssm / c->dfm;
 282   const double mse = pspp_linreg_mse (c);
 283   const double F = msm / mse;
 284   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
 285
 286   struct tab_table *t;
 287
 288   assert (c != NULL);
 289   t = tab_create (n_cols, n_rows);
 290   tab_headers (t, 2, 0, 1, 0);
 291
 292   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 293
 294   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 295   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 296   tab_vline (t, TAL_0, 1, 0, 0);
 297
 298   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 299   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 300   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 301   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 302   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 303
 304   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
 305   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
 306   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
 307
 308   /* Sums of Squares */
 309   tab_double (t, 2, 1, 0, c->ssm, NULL);
 310   tab_double (t, 2, 3, 0, c->sst, NULL);
 311   tab_double (t, 2, 2, 0, c->sse, NULL);
 312
 313
 314   /* Degrees of freedom */
 315   tab_text_format (t, 3, 1, TAB_RIGHT, "%g", c->dfm);
 316   tab_text_format (t, 3, 2, TAB_RIGHT, "%g", c->dfe);
 317   tab_text_format (t, 3, 3, TAB_RIGHT, "%g", c->dft);
 318
 319   /* Mean Squares */
 320   tab_double (t, 4, 1, TAB_RIGHT, msm, NULL);
 321   tab_double (t, 4, 2, TAB_RIGHT, mse, NULL);
 322
 323   tab_double (t, 5, 1, 0, F, NULL);
 324
 325   tab_double (t, 6, 1, 0, pval, NULL);
 326
 327   tab_title (t, _("ANOVA"));
 328   tab_submit (t);
 329 }
 330
 331 static void
 332 reg_stats_outs (pspp_linreg_cache * c)
 333 {
 334   assert (c != NULL);
 335 }
 336
 337 static void
 338 reg_stats_zpp (pspp_linreg_cache * c)
 339 {
 340   assert (c != NULL);
 341 }
 342
 343 static void
 344 reg_stats_label (pspp_linreg_cache * c)
 345 {
 346   assert (c != NULL);
 347 }
 348
 349 static void
 350 reg_stats_sha (pspp_linreg_cache * c)
 351 {
 352   assert (c != NULL);
 353 }
 354 static void
 355 reg_stats_ci (pspp_linreg_cache * c)
 356 {
 357   assert (c != NULL);
 358 }
 359 static void
 360 reg_stats_f (pspp_linreg_cache * c)
 361 {
 362   assert (c != NULL);
 363 }
 364 static void
 365 reg_stats_bcov (pspp_linreg_cache * c)
 366 {
 367   int n_cols;
 368   int n_rows;
 369   int i;
 370   int k;
 371   int row;
 372   int col;
 373   const char *label;
 374   struct tab_table *t;
 375
 376   assert (c != NULL);
 377   n_cols = c->n_indeps + 1 + 2;
 378   n_rows = 2 * (c->n_indeps + 1);
 379   t = tab_create (n_cols, n_rows);
 380   tab_headers (t, 2, 0, 1, 0);
 381   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 382   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 383   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 384   tab_vline (t, TAL_0, 1, 0, 0);
 385   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 386   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 387   for (i = 0; i < c->n_coeffs; i++)
 388     {
 389       const struct variable *v = pspp_coeff_get_var (c->coeff[i], 0);
 390       label = var_to_string (v);
 391       tab_text (t, 2, i, TAB_CENTER, label);
 392       tab_text (t, i + 2, 0, TAB_CENTER, label);
 393       for (k = 1; k < c->n_coeffs; k++)
 394         {
 395           col = (i <= k) ? k : i;
 396           row = (i <= k) ? i : k;
 397           tab_double (t, k + 2, i, TAB_CENTER,
 398                      gsl_matrix_get (c->cov, row, col), NULL);
 399         }
 400     }
 401   tab_title (t, _("Coefficient Correlations"));
 402   tab_submit (t);
 403 }
 404 static void
 405 reg_stats_ses (pspp_linreg_cache * c)
 406 {
 407   assert (c != NULL);
 408 }
 409 static void
 410 reg_stats_xtx (pspp_linreg_cache * c)
 411 {
 412   assert (c != NULL);
 413 }
 414 static void
 415 reg_stats_collin (pspp_linreg_cache * c)
 416 {
 417   assert (c != NULL);
 418 }
 419 static void
 420 reg_stats_tol (pspp_linreg_cache * c)
 421 {
 422   assert (c != NULL);
 423 }
 424 static void
 425 reg_stats_selection (pspp_linreg_cache * c)
 426 {
 427   assert (c != NULL);
 428 }
 429
 430 static void
 431 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 432                            int keyword, pspp_linreg_cache * c)
 433 {
 434   if (keyword)
 435     {
 436       (*function) (c);
 437     }
 438 }
 439
 440 static void
 441 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 442 {
 443   /*
 444      The order here must match the order in which the STATISTICS
 445      keywords appear in the specification section above.
 446    */
 447   enum
 448   { r,
 449     coeff,
 450     anova,
 451     outs,
 452     zpp,
 453     label,
 454     sha,
 455     ci,
 456     bcov,
 457     ses,
 458     xtx,
 459     collin,
 460     tol,
 461     selection,
 462     f,
 463     defaults,
 464     all
 465   };
 466   int i;
 467   int d = 1;
 468
 469   if (keywords[all])
 470     {
 471       /*
 472          Set everything but F.
 473        */
 474       for (i = 0; i < f; i++)
 475         {
 476           keywords[i] = 1;
 477         }
 478     }
 479   else
 480     {
 481       for (i = 0; i < all; i++)
 482         {
 483           if (keywords[i])
 484             {
 485               d = 0;
 486             }
 487         }
 488       /*
 489          Default output: ANOVA table, parameter estimates,
 490          and statistics for variables not entered into model,
 491          if appropriate.
 492        */
 493       if (keywords[defaults] | d)
 494         {
 495           keywords[anova] = 1;
 496           keywords[outs] = 1;
 497           keywords[coeff] = 1;
 498           keywords[r] = 1;
 499         }
 500     }
 501   statistics_keyword_output (reg_stats_r, keywords[r], c);
 502   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 503   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 504   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 505   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 506   statistics_keyword_output (reg_stats_label, keywords[label], c);
 507   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 508   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 509   statistics_keyword_output (reg_stats_f, keywords[f], c);
 510   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 511   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 512   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 513   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 514   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 515   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 516 }
 517
 518 /*
 519   Free the transformation. Free its linear model if this
 520   transformation is the last one.
 521  */
 522 static bool
 523 regression_trns_free (void *t_)
 524 {
 525   bool result = true;
 526   struct reg_trns *t = t_;
 527
 528   if (t->trns_id == t->n_trns)
 529     {
 530       result = pspp_linreg_cache_free (t->c);
 531     }
 532   free (t);
 533
 534   return result;
 535 }
 536
 537 /*
 538   Gets the predicted values.
 539  */
 540 static int
 541 regression_trns_pred_proc (void *t_, struct ccase **c,
 542                            casenumber case_idx UNUSED)
 543 {
 544   size_t i;
 545   size_t n_vals;
 546   struct reg_trns *trns = t_;
 547   pspp_linreg_cache *model;
 548   union value *output = NULL;
 549   const union value **vals = NULL;
 550   const struct variable **vars = NULL;
 551
 552   assert (trns != NULL);
 553   model = trns->c;
 554   assert (model != NULL);
 555   assert (model->depvar != NULL);
 556   assert (model->pred != NULL);
 557
 558   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 559   n_vals = (*model->get_vars) (model, vars);
 560
 561   vals = xnmalloc (n_vals, sizeof (*vals));
 562   *c = case_unshare (*c);
 563   output = case_data_rw (*c, model->pred);
 564
 565   for (i = 0; i < n_vals; i++)
 566     {
 567       vals[i] = case_data (*c, vars[i]);
 568     }
 569   output->f = (*model->predict) ((const struct variable **) vars,
 570                                  vals, model, n_vals);
 571   free (vals);
 572   free (vars);
 573   return TRNS_CONTINUE;
 574 }
 575
 576 /*
 577   Gets the residuals.
 578  */
 579 static int
 580 regression_trns_resid_proc (void *t_, struct ccase **c,
 581                             casenumber case_idx UNUSED)
 582 {
 583   size_t i;
 584   size_t n_vals;
 585   struct reg_trns *trns = t_;
 586   pspp_linreg_cache *model;
 587   union value *output = NULL;
 588   const union value **vals = NULL;
 589   const union value *obs = NULL;
 590   const struct variable **vars = NULL;
 591
 592   assert (trns != NULL);
 593   model = trns->c;
 594   assert (model != NULL);
 595   assert (model->depvar != NULL);
 596   assert (model->resid != NULL);
 597
 598   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 599   n_vals = (*model->get_vars) (model, vars);
 600
 601   vals = xnmalloc (n_vals, sizeof (*vals));
 602   *c = case_unshare (*c);
 603   output = case_data_rw (*c, model->resid);
 604   assert (output != NULL);
 605
 606   for (i = 0; i < n_vals; i++)
 607     {
 608       vals[i] = case_data (*c, vars[i]);
 609     }
 610   obs = case_data (*c, model->depvar);
 611   output->f = (*model->residual) ((const struct variable **) vars,
 612                                   vals, obs, model, n_vals);
 613   free (vals);
 614   free (vars);
 615   return TRNS_CONTINUE;
 616 }
 617
 618 /*
 619    Returns false if NAME is a duplicate of any existing variable name.
 620 */
 621 static bool
 622 try_name (const struct dictionary *dict, const char *name)
 623 {
 624   if (dict_lookup_var (dict, name) != NULL)
 625     return false;
 626
 627   return true;
 628 }
 629
 630 static void
 631 reg_get_name (const struct dictionary *dict, char name[VAR_NAME_LEN],
 632               const char prefix[VAR_NAME_LEN])
 633 {
 634   int i = 1;
 635
 636   snprintf (name, VAR_NAME_LEN, "%s%d", prefix, i);
 637   while (!try_name (dict, name))
 638     {
 639       i++;
 640       snprintf (name, VAR_NAME_LEN, "%s%d", prefix, i);
 641     }
 642 }
 643
 644 static void
 645 reg_save_var (struct dataset *ds, const char *prefix, trns_proc_func * f,
 646               pspp_linreg_cache * c, struct variable **v, int n_trns)
 647 {
 648   struct dictionary *dict = dataset_dict (ds);
 649   static int trns_index = 1;
 650   char name[VAR_NAME_LEN];
 651   struct variable *new_var;
 652   struct reg_trns *t = NULL;
 653
 654   t = xmalloc (sizeof (*t));
 655   t->trns_id = trns_index;
 656   t->n_trns = n_trns;
 657   t->c = c;
 658   reg_get_name (dict, name, prefix);
 659   new_var = dict_create_var (dict, name, 0);
 660   assert (new_var != NULL);
 661   *v = new_var;
 662   add_transformation (ds, f, regression_trns_free, t);
 663   trns_index++;
 664 }
 665 static void
 666 subcommand_save (struct dataset *ds, int save, pspp_linreg_cache ** models)
 667 {
 668   pspp_linreg_cache **lc;
 669   int n_trns = 0;
 670   int i;
 671
 672   assert (models != NULL);
 673
 674   if (save)
 675     {
 676       /* Count the number of transformations we will need. */
 677       for (i = 0; i < REGRESSION_SV_count; i++)
 678         {
 679           if (cmd.a_save[i])
 680             {
 681               n_trns++;
 682             }
 683         }
 684       n_trns *= cmd.n_dependent;
 685
 686       for (lc = models; lc < models + cmd.n_dependent; lc++)
 687         {
 688           if (*lc != NULL)
 689             {
 690               if ((*lc)->depvar != NULL)
 691                 {
 692                   if (cmd.a_save[REGRESSION_SV_RESID])
 693                     {
 694                       reg_save_var (ds, "RES", regression_trns_resid_proc, *lc,
 695                                     &(*lc)->resid, n_trns);
 696                     }
 697                   if (cmd.a_save[REGRESSION_SV_PRED])
 698                     {
 699                       reg_save_var (ds, "PRED", regression_trns_pred_proc, *lc,
 700                                     &(*lc)->pred, n_trns);
 701                     }
 702                 }
 703             }
 704         }
 705     }
 706   else
 707     {
 708       for (lc = models; lc < models + cmd.n_dependent; lc++)
 709         {
 710           if (*lc != NULL)
 711             {
 712               pspp_linreg_cache_free (*lc);
 713             }
 714         }
 715     }
 716 }
 717
 718 int
 719 cmd_regression (struct lexer *lexer, struct dataset *ds)
 720 {
 721   struct casegrouper *grouper;
 722   struct casereader *group;
 723   pspp_linreg_cache **models;
 724   bool ok;
 725   size_t i;
 726
 727   if (!parse_regression (lexer, ds, &cmd, NULL))
 728     {
 729       return CMD_FAILURE;
 730     }
 731
 732   models = xnmalloc (cmd.n_dependent, sizeof *models);
 733   for (i = 0; i < cmd.n_dependent; i++)
 734     {
 735       models[i] = NULL;
 736     }
 737
 738   /* Data pass. */
 739   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 740   while (casegrouper_get_next_group (grouper, &group))
 741     run_regression (group, &cmd, ds, models);
 742   ok = casegrouper_destroy (grouper);
 743   ok = proc_commit (ds) && ok;
 744
 745   subcommand_save (ds, cmd.sbc_save, models);
 746   free (v_variables);
 747   free (models);
 748   free_regression (&cmd);
 749
 750   return ok ? CMD_SUCCESS : CMD_FAILURE;
 751 }
 752
 753 /*
 754   Is variable k the dependent variable?
 755  */
 756 static bool
 757 is_depvar (size_t k, const struct variable *v)
 758 {
 759   return v == v_variables[k];
 760 }
 761
 762 /* Parser for the variables sub command */
 763 static int
 764 regression_custom_variables (struct lexer *lexer, struct dataset *ds,
 765                              struct cmd_regression *cmd UNUSED,
 766                              void *aux UNUSED)
 767 {
 768   const struct dictionary *dict = dataset_dict (ds);
 769
 770   lex_match (lexer, '=');
 771
 772   if ((lex_token (lexer) != T_ID
 773        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 774       && lex_token (lexer) != T_ALL)
 775     return 2;
 776
 777
 778   if (!parse_variables_const
 779       (lexer, dict, &v_variables, &n_variables, PV_NONE))
 780     {
 781       free (v_variables);
 782       return 0;
 783     }
 784   assert (n_variables);
 785
 786   return 1;
 787 }
 788
 789 /* Identify the explanatory variables in v_variables.  Returns
 790    the number of independent variables. */
 791 static int
 792 identify_indep_vars (const struct variable **indep_vars,
 793                      const struct variable *depvar)
 794 {
 795   int n_indep_vars = 0;
 796   int i;
 797
 798   for (i = 0; i < n_variables; i++)
 799     if (!is_depvar (i, depvar))
 800       indep_vars[n_indep_vars++] = v_variables[i];
 801   if ((n_indep_vars < 1) && is_depvar (0, depvar))
 802     {
 803       /*
 804         There is only one independent variable, and it is the same
 805         as the dependent variable. Print a warning and continue.
 806        */
 807       msg (SE,
 808            gettext ("The dependent variable is equal to the independent variable."
 809                     "The least squares line is therefore Y=X."
 810                     "Standard errors and related statistics may be meaningless."));
 811       n_indep_vars = 1;
 812       indep_vars[0] = v_variables[0];
 813     }
 814   return n_indep_vars;
 815 }
 816
 817 /* Encode categorical variables.
 818    Returns number of valid cases. */
 819 static int
 820 prepare_categories (struct casereader *input,
 821                     const struct variable **vars, size_t n_vars,
 822                     struct moments_var *mom)
 823 {
 824   int n_data;
 825   struct ccase *c;
 826   size_t i;
 827
 828   assert (vars != NULL);
 829   assert (mom != NULL);
 830
 831   for (i = 0; i < n_vars; i++)
 832     if (var_is_alpha (vars[i]))
 833       cat_stored_values_create (vars[i]);
 834
 835   n_data = 0;
 836   for (; (c = casereader_read (input)) != NULL; case_unref (c))
 837     {
 838       /*
 839          The second condition ensures the program will run even if
 840          there is only one variable to act as both explanatory and
 841          response.
 842        */
 843       for (i = 0; i < n_vars; i++)
 844         {
 845           const union value *val = case_data (c, vars[i]);
 846           if (var_is_alpha (vars[i]))
 847             cat_value_update (vars[i], val);
 848           else
 849             moments1_add (mom[i].m, val->f, 1.0);
 850         }
 851       n_data++;
 852     }
 853   casereader_destroy (input);
 854
 855   return n_data;
 856 }
 857
 858 static void
 859 coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
 860 {
 861   c->coeff = xnmalloc (dm->m->size2, sizeof (*c->coeff));
 862   pspp_coeff_init (c->coeff, dm);
 863 }
 864
 865 static bool
 866 run_regression (struct casereader *input, struct cmd_regression *cmd,
 867                 struct dataset *ds, pspp_linreg_cache **models)
 868 {
 869   size_t i;
 870   int n_indep = 0;
 871   int k;
 872   struct ccase *c;
 873   const struct variable **indep_vars;
 874   struct design_matrix *X;
 875   struct moments_var *mom;
 876   gsl_vector *Y;
 877
 878   pspp_linreg_opts lopts;
 879
 880   assert (models != NULL);
 881
 882   c = casereader_peek (input, 0);
 883   if (c == NULL)
 884     {
 885       casereader_destroy (input);
 886       return true;
 887     }
 888   output_split_file_values (ds, c);
 889   case_unref (c);
 890
 891   if (!v_variables)
 892     {
 893       dict_get_vars (dataset_dict (ds), &v_variables, &n_variables, 0);
 894     }
 895
 896   for (i = 0; i < cmd->n_dependent; i++)
 897     {
 898       if (!var_is_numeric (cmd->v_dependent[i]))
 899         {
 900           msg (SE, _("Dependent variable must be numeric."));
 901           return false;
 902         }
 903     }
 904
 905   mom = xnmalloc (n_variables, sizeof (*mom));
 906   for (i = 0; i < n_variables; i++)
 907     {
 908       (mom + i)->m = moments1_create (MOMENT_VARIANCE);
 909       (mom + i)->v = v_variables[i];
 910     }
 911   lopts.get_depvar_mean_std = 1;
 912
 913
 914   indep_vars = xnmalloc (n_variables, sizeof *indep_vars);
 915
 916   for (k = 0; k < cmd->n_dependent; k++)
 917     {
 918       const struct variable *dep_var;
 919       struct casereader *reader;
 920       casenumber row;
 921       struct ccase *c;
 922       size_t n_data;            /* Number of valid cases. */
 923
 924       dep_var = cmd->v_dependent[k];
 925       n_indep = identify_indep_vars (indep_vars, dep_var);
 926       reader = casereader_clone (input);
 927       reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
 928                                                  MV_ANY, NULL, NULL);
 929       reader = casereader_create_filter_missing (reader, &dep_var, 1,
 930                                                  MV_ANY, NULL, NULL);
 931       n_data = prepare_categories (casereader_clone (reader),
 932                                    indep_vars, n_indep, mom);
 933
 934       if ((n_data > 0) && (n_indep > 0))
 935         {
 936           Y = gsl_vector_alloc (n_data);
 937           X =
 938             design_matrix_create (n_indep,
 939                                   (const struct variable **) indep_vars,
 940                                   n_data);
 941           lopts.get_indep_mean_std = xnmalloc (X->m->size2, sizeof (int));
 942           for (i = 0; i < X->m->size2; i++)
 943             {
 944               lopts.get_indep_mean_std[i] = 1;
 945             }
 946           models[k] = pspp_linreg_cache_alloc (dep_var, (const struct variable **) indep_vars,
 947                                                X->m->size1, n_indep);
 948           models[k]->depvar = dep_var;
 949           /*
 950              For large data sets, use QR decomposition.
 951            */
 952           if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
 953             {
 954               models[k]->method = PSPP_LINREG_QR;
 955             }
 956
 957           /*
 958              The second pass fills the design matrix.
 959            */
 960           reader = casereader_create_counter (reader, &row, -1);
 961           for (; (c = casereader_read (reader)) != NULL; case_unref (c))
 962             {
 963               for (i = 0; i < n_indep; ++i)
 964                 {
 965                   const struct variable *v = indep_vars[i];
 966                   const union value *val = case_data (c, v);
 967                   if (var_is_alpha (v))
 968                     design_matrix_set_categorical (X, row, v, val);
 969                   else
 970                     design_matrix_set_numeric (X, row, v, val);
 971                 }
 972               gsl_vector_set (Y, row, case_num (c, dep_var));
 973             }
 974           /*
 975              Now that we know the number of coefficients, allocate space
 976              and store pointers to the variables that correspond to the
 977              coefficients.
 978            */
 979           coeff_init (models[k], X);
 980
 981           /*
 982              Find the least-squares estimates and other statistics.
 983            */
 984           pspp_linreg ((const gsl_vector *) Y, X, &lopts, models[k]);
 985
 986           if (!taint_has_tainted_successor (casereader_get_taint (input)))
 987             {
 988               subcommand_statistics (cmd->a_statistics, models[k]);
 989             }
 990
 991           gsl_vector_free (Y);
 992           design_matrix_destroy (X);
 993           free (lopts.get_indep_mean_std);
 994         }
 995       else
 996         {
 997           msg (SE,
 998                gettext ("No valid data found. This command was skipped."));
 999         }
1000       casereader_destroy (reader);
1001     }
1002   for (i = 0; i < n_variables; i++)
1003     {
1004       moments1_destroy ((mom + i)->m);
1005     }
1006   free (mom);
1007   free (indep_vars);
1008   casereader_destroy (input);
1009
1010   return true;
1011 }
1012
1013 /*
1014   Local Variables:
1015   mode: c
1016   End:
1017 */