pintos-os.org Git - pspp/blob - 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
 796 static double
 797 fill_covariance (gsl_matrix *cov, struct covariance *all_cov,
 798                  const struct variable **vars,
 799                  size_t n_vars, const struct variable *dep_var,
 800                  const struct variable **all_vars, size_t n_all_vars,
 801                  double *means)
 802 {
 803   size_t i;
 804   size_t j;
 805   size_t dep_subscript;
 806   size_t *rows;
 807   const gsl_matrix *ssizes;
 808   const gsl_matrix *mean_matrix;
 809   const gsl_matrix *ssize_matrix;
 810   double result = 0.0;
 811
 812   gsl_matrix *cm = covariance_calculate_unnormalized (all_cov);
 813
 814   if ( cm == NULL)
 815     return 0;
 816
 817   rows = xnmalloc (cov->size1 - 1, sizeof (*rows));
 818
 819   for (i = 0; i < n_all_vars; i++)
 820     {
 821       for (j = 0; j < n_vars; j++)
 822         {
 823           if (vars[j] == all_vars[i])
 824             {
 825               rows[j] = i;
 826             }
 827         }
 828       if (all_vars[i] == dep_var)
 829         {
 830           dep_subscript = i;
 831         }
 832     }
 833   mean_matrix = covariance_moments (all_cov, MOMENT_MEAN);
 834   ssize_matrix = covariance_moments (all_cov, MOMENT_NONE);
 835   for (i = 0; i < cov->size1 - 1; i++)
 836     {
 837       means[i] = gsl_matrix_get (mean_matrix, rows[i], 0)
 838         / gsl_matrix_get (ssize_matrix, rows[i], 0);
 839       for (j = 0; j < cov->size2 - 1; j++)
 840         {
 841           gsl_matrix_set (cov, i, j, gsl_matrix_get (cm, rows[i], rows[j]));
 842           gsl_matrix_set (cov, j, i, gsl_matrix_get (cm, rows[j], rows[i]));
 843         }
 844     }
 845   means[cov->size1 - 1] = gsl_matrix_get (mean_matrix, dep_subscript, 0)
 846     / gsl_matrix_get (ssize_matrix, dep_subscript, 0);
 847   ssizes = covariance_moments (all_cov, MOMENT_NONE);
 848   result = gsl_matrix_get (ssizes, dep_subscript, rows[0]);
 849   for (i = 0; i < cov->size1 - 1; i++)
 850     {
 851       gsl_matrix_set (cov, i, cov->size1 - 1,
 852                       gsl_matrix_get (cm, rows[i], dep_subscript));
 853       gsl_matrix_set (cov, cov->size1 - 1, i,
 854                       gsl_matrix_get (cm, rows[i], dep_subscript));
 855       if (result > gsl_matrix_get (ssizes, rows[i], dep_subscript))
 856         {
 857           result = gsl_matrix_get (ssizes, rows[i], dep_subscript);
 858         }
 859     }
 860   gsl_matrix_set (cov, cov->size1 - 1, cov->size1 - 1,
 861                   gsl_matrix_get (cm, dep_subscript, dep_subscript));
 862   free (rows);
 863   gsl_matrix_free (cm);
 864   return result;
 865 }
 866 static size_t
 867 get_n_all_vars (struct cmd_regression *cmd)
 868 {
 869   size_t result = n_variables;
 870   size_t i;
 871   size_t j;
 872
 873   result += cmd->n_dependent;
 874   for (i = 0; i < cmd->n_dependent; i++)
 875     {
 876       for (j = 0; j < n_variables; j++)
 877         {
 878           if (v_variables[j] == cmd->v_dependent[i])
 879             {
 880               result--;
 881             }
 882         }
 883     }
 884   return result;
 885 }
 886 static void
 887 fill_all_vars (const struct variable **vars, struct cmd_regression *cmd)
 888 {
 889   size_t i;
 890   size_t j;
 891   bool absent;
 892
 893   for (i = 0; i < n_variables; i++)
 894     {
 895       vars[i] = v_variables[i];
 896     }
 897   for (i = 0; i < cmd->n_dependent; i++)
 898     {
 899       absent = true;
 900       for (j = 0; j < n_variables; j++)
 901         {
 902           if (cmd->v_dependent[i] == v_variables[j])
 903             {
 904               absent = false;
 905               break;
 906             }
 907         }
 908       if (absent)
 909         {
 910           vars[i + n_variables] = cmd->v_dependent[i];
 911         }
 912     }
 913 }
 914 static bool
 915 run_regression (struct casereader *input, struct cmd_regression *cmd,
 916                 struct dataset *ds, linreg **models)
 917 {
 918   size_t i;
 919   int n_indep = 0;
 920   int k;
 921   double n_data;
 922   double *means;
 923   struct ccase *c;
 924   struct covariance *cov;
 925   const struct variable **vars;
 926   const struct variable **all_vars;
 927   const struct variable *dep_var;
 928   struct casereader *reader;
 929   const struct dictionary *dict;
 930   size_t n_all_vars;
 931
 932   assert (models != NULL);
 933
 934   for (i = 0; i < n_variables; i++)
 935     {
 936       if (!var_is_numeric (v_variables[i]))
 937         {
 938           msg (SE, _("REGRESSION requires numeric variables."));
 939           return false;
 940         }
 941     }
 942
 943   c = casereader_peek (input, 0);
 944   if (c == NULL)
 945     {
 946       casereader_destroy (input);
 947       return true;
 948     }
 949   output_split_file_values (ds, c);
 950   case_unref (c);
 951
 952   dict = dataset_dict (ds);
 953   if (!v_variables)
 954     {
 955       dict_get_vars (dict, &v_variables, &n_variables, 0);
 956     }
 957   n_all_vars = get_n_all_vars (cmd);
 958   all_vars = xnmalloc (n_all_vars, sizeof (*all_vars));
 959   fill_all_vars (all_vars, cmd);
 960   vars = xnmalloc (n_variables, sizeof (*vars));
 961   means  = xnmalloc (n_all_vars, sizeof (*means));
 962   cov = covariance_1pass_create (n_all_vars, all_vars,
 963                                  dict_get_weight (dict), MV_ANY);
 964
 965   reader = casereader_clone (input);
 966   reader = casereader_create_filter_missing (reader, v_variables, n_variables,
 967                                              MV_ANY, NULL, NULL);
 968   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
 969     {
 970       covariance_accumulate (cov, c);
 971     }
 972
 973   for (k = 0; k < cmd->n_dependent; k++)
 974     {
 975       gsl_matrix *this_cm;
 976       dep_var = cmd->v_dependent[k];
 977       n_indep = identify_indep_vars (vars, dep_var);
 978
 979       this_cm = gsl_matrix_alloc (n_indep + 1, n_indep + 1);
 980       n_data = fill_covariance (this_cm, cov, vars, n_indep,
 981                                 dep_var, all_vars, n_all_vars, means);
 982       models[k] = linreg_alloc (dep_var, (const struct variable **) vars,
 983                                 n_data, n_indep);
 984       models[k]->depvar = dep_var;
 985       for (i = 0; i < n_indep; i++)
 986         {
 987           linreg_set_indep_variable_mean (models[k], i, means[i]);
 988         }
 989       linreg_set_depvar_mean (models[k], means[i]);
 990       /*
 991         For large data sets, use QR decomposition.
 992       */
 993       if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
 994         {
 995           models[k]->method = LINREG_QR;
 996         }
 997
 998       if (n_data > 0)
 999         {
1000           /*
1001             Find the least-squares estimates and other statistics.
1002           */
1003           linreg_fit (this_cm, models[k]);
1004
1005           if (!taint_has_tainted_successor (casereader_get_taint (input)))
1006             {
1007               subcommand_statistics (cmd->a_statistics, models[k], this_cm);
1008             }
1009         }
1010       else
1011         {
1012           msg (SE,
1013                gettext ("No valid data found. This command was skipped."));
1014           linreg_free (models[k]);
1015           models[k] = NULL;
1016         }
1017       gsl_matrix_free (this_cm);
1018     }
1019
1020   casereader_destroy (reader);
1021   free (vars);
1022   free (all_vars);
1023   free (means);
1024   casereader_destroy (input);
1025   covariance_destroy (cov);
1026
1027   return true;
1028 }
1029
1030 /*
1031   Local Variables:
1032   mode: c
1033   End:
1034 */