src/language/stats/regression.q

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