src/language/stats/regression.q

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