src/language/stats/regression.q

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