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       struct string tstr;
 221       ds_init_empty (&tstr);
 222       this_row = j + 2;
 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_text (t, 3, 1, TAB_RIGHT | TAT_PRINTF, "%g", c->dfm);
 322   tab_text (t, 3, 2, TAB_RIGHT | TAT_PRINTF, "%g", c->dfe);
 323   tab_text (t, 3, 3, TAB_RIGHT | TAT_PRINTF, "%g", c->dft);
 324
 325   /* Mean Squares */
 326   tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
 327   tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
 328
 329   tab_float (t, 5, 1, 0, F, 8, 3);
 330
 331   tab_float (t, 6, 1, 0, pval, 8, 3);
 332
 333   tab_title (t, _("ANOVA"));
 334   tab_submit (t);
 335 }
 336
 337 static void
 338 reg_stats_outs (pspp_linreg_cache * c)
 339 {
 340   assert (c != NULL);
 341 }
 342
 343 static void
 344 reg_stats_zpp (pspp_linreg_cache * c)
 345 {
 346   assert (c != NULL);
 347 }
 348
 349 static void
 350 reg_stats_label (pspp_linreg_cache * c)
 351 {
 352   assert (c != NULL);
 353 }
 354
 355 static void
 356 reg_stats_sha (pspp_linreg_cache * c)
 357 {
 358   assert (c != NULL);
 359 }
 360 static void
 361 reg_stats_ci (pspp_linreg_cache * c)
 362 {
 363   assert (c != NULL);
 364 }
 365 static void
 366 reg_stats_f (pspp_linreg_cache * c)
 367 {
 368   assert (c != NULL);
 369 }
 370 static void
 371 reg_stats_bcov (pspp_linreg_cache * c)
 372 {
 373   int n_cols;
 374   int n_rows;
 375   int i;
 376   int k;
 377   int row;
 378   int col;
 379   const char *label;
 380   struct tab_table *t;
 381
 382   assert (c != NULL);
 383   n_cols = c->n_indeps + 1 + 2;
 384   n_rows = 2 * (c->n_indeps + 1);
 385   t = tab_create (n_cols, n_rows, 0);
 386   tab_headers (t, 2, 0, 1, 0);
 387   tab_dim (t, tab_natural_dimensions);
 388   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 389   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 390   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 391   tab_vline (t, TAL_0, 1, 0, 0);
 392   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 393   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 394   for (i = 0; i < c->n_coeffs; i++)
 395     {
 396       const struct variable *v = pspp_coeff_get_var (c->coeff[i], 0);
 397       label = var_to_string (v);
 398       tab_text (t, 2, i, TAB_CENTER, label);
 399       tab_text (t, i + 2, 0, TAB_CENTER, label);
 400       for (k = 1; k < c->n_coeffs; k++)
 401         {
 402           col = (i <= k) ? k : i;
 403           row = (i <= k) ? i : k;
 404           tab_float (t, k + 2, i, TAB_CENTER,
 405                      gsl_matrix_get (c->cov, row, col), 8, 3);
 406         }
 407     }
 408   tab_title (t, _("Coefficient Correlations"));
 409   tab_submit (t);
 410 }
 411 static void
 412 reg_stats_ses (pspp_linreg_cache * c)
 413 {
 414   assert (c != NULL);
 415 }
 416 static void
 417 reg_stats_xtx (pspp_linreg_cache * c)
 418 {
 419   assert (c != NULL);
 420 }
 421 static void
 422 reg_stats_collin (pspp_linreg_cache * c)
 423 {
 424   assert (c != NULL);
 425 }
 426 static void
 427 reg_stats_tol (pspp_linreg_cache * c)
 428 {
 429   assert (c != NULL);
 430 }
 431 static void
 432 reg_stats_selection (pspp_linreg_cache * c)
 433 {
 434   assert (c != NULL);
 435 }
 436
 437 static void
 438 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 439                            int keyword, pspp_linreg_cache * c)
 440 {
 441   if (keyword)
 442     {
 443       (*function) (c);
 444     }
 445 }
 446
 447 static void
 448 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 449 {
 450   /*
 451      The order here must match the order in which the STATISTICS
 452      keywords appear in the specification section above.
 453    */
 454   enum
 455   { r,
 456     coeff,
 457     anova,
 458     outs,
 459     zpp,
 460     label,
 461     sha,
 462     ci,
 463     bcov,
 464     ses,
 465     xtx,
 466     collin,
 467     tol,
 468     selection,
 469     f,
 470     defaults,
 471     all
 472   };
 473   int i;
 474   int d = 1;
 475
 476   if (keywords[all])
 477     {
 478       /*
 479          Set everything but F.
 480        */
 481       for (i = 0; i < f; i++)
 482         {
 483           keywords[i] = 1;
 484         }
 485     }
 486   else
 487     {
 488       for (i = 0; i < all; i++)
 489         {
 490           if (keywords[i])
 491             {
 492               d = 0;
 493             }
 494         }
 495       /*
 496          Default output: ANOVA table, parameter estimates,
 497          and statistics for variables not entered into model,
 498          if appropriate.
 499        */
 500       if (keywords[defaults] | d)
 501         {
 502           keywords[anova] = 1;
 503           keywords[outs] = 1;
 504           keywords[coeff] = 1;
 505           keywords[r] = 1;
 506         }
 507     }
 508   statistics_keyword_output (reg_stats_r, keywords[r], c);
 509   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 510   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 511   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 512   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 513   statistics_keyword_output (reg_stats_label, keywords[label], c);
 514   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 515   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 516   statistics_keyword_output (reg_stats_f, keywords[f], c);
 517   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 518   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 519   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 520   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 521   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 522   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 523 }
 524
 525 /*
 526   Free the transformation. Free its linear model if this
 527   transformation is the last one.
 528  */
 529 static bool
 530 regression_trns_free (void *t_)
 531 {
 532   bool result = true;
 533   struct reg_trns *t = t_;
 534
 535   if (t->trns_id == t->n_trns)
 536     {
 537       result = pspp_linreg_cache_free (t->c);
 538     }
 539   free (t);
 540
 541   return result;
 542 }
 543
 544 /*
 545   Gets the predicted values.
 546  */
 547 static int
 548 regression_trns_pred_proc (void *t_, struct ccase *c,
 549                            casenumber case_idx UNUSED)
 550 {
 551   size_t i;
 552   size_t n_vals;
 553   struct reg_trns *trns = t_;
 554   pspp_linreg_cache *model;
 555   union value *output = NULL;
 556   const union value **vals = NULL;
 557   const struct variable **vars = NULL;
 558
 559   assert (trns != NULL);
 560   model = trns->c;
 561   assert (model != NULL);
 562   assert (model->depvar != NULL);
 563   assert (model->pred != NULL);
 564
 565   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 566   n_vals = (*model->get_vars) (model, vars);
 567
 568   vals = xnmalloc (n_vals, sizeof (*vals));
 569   output = case_data_rw (c, model->pred);
 570   assert (output != NULL);
 571
 572   for (i = 0; i < n_vals; i++)
 573     {
 574       vals[i] = case_data (c, vars[i]);
 575     }
 576   output->f = (*model->predict) ((const struct variable **) vars,
 577                                  vals, model, n_vals);
 578   free (vals);
 579   free (vars);
 580   return TRNS_CONTINUE;
 581 }
 582
 583 /*
 584   Gets the residuals.
 585  */
 586 static int
 587 regression_trns_resid_proc (void *t_, struct ccase *c,
 588                             casenumber case_idx UNUSED)
 589 {
 590   size_t i;
 591   size_t n_vals;
 592   struct reg_trns *trns = t_;
 593   pspp_linreg_cache *model;
 594   union value *output = NULL;
 595   const union value **vals = NULL;
 596   const union value *obs = NULL;
 597   const struct variable **vars = NULL;
 598
 599   assert (trns != NULL);
 600   model = trns->c;
 601   assert (model != NULL);
 602   assert (model->depvar != NULL);
 603   assert (model->resid != NULL);
 604
 605   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 606   n_vals = (*model->get_vars) (model, vars);
 607
 608   vals = xnmalloc (n_vals, sizeof (*vals));
 609   output = case_data_rw (c, model->resid);
 610   assert (output != NULL);
 611
 612   for (i = 0; i < n_vals; i++)
 613     {
 614       vals[i] = case_data (c, vars[i]);
 615     }
 616   obs = case_data (c, model->depvar);
 617   output->f = (*model->residual) ((const struct variable **) vars,
 618                                   vals, obs, model, n_vals);
 619   free (vals);
 620   free (vars);
 621   return TRNS_CONTINUE;
 622 }
 623
 624 /*
 625    Returns false if NAME is a duplicate of any existing variable name.
 626 */
 627 static bool
 628 try_name (const struct dictionary *dict, const char *name)
 629 {
 630   if (dict_lookup_var (dict, name) != NULL)
 631     return false;
 632
 633   return true;
 634 }
 635
 636 static void
 637 reg_get_name (const struct dictionary *dict, char name[VAR_NAME_LEN],
 638               const char prefix[VAR_NAME_LEN])
 639 {
 640   int i = 1;
 641
 642   snprintf (name, VAR_NAME_LEN, "%s%d", prefix, i);
 643   while (!try_name (dict, name))
 644     {
 645       i++;
 646       snprintf (name, VAR_NAME_LEN, "%s%d", prefix, i);
 647     }
 648 }
 649
 650 static void
 651 reg_save_var (struct dataset *ds, const char *prefix, trns_proc_func * f,
 652               pspp_linreg_cache * c, struct variable **v, int n_trns)
 653 {
 654   struct dictionary *dict = dataset_dict (ds);
 655   static int trns_index = 1;
 656   char name[VAR_NAME_LEN];
 657   struct variable *new_var;
 658   struct reg_trns *t = NULL;
 659
 660   t = xmalloc (sizeof (*t));
 661   t->trns_id = trns_index;
 662   t->n_trns = n_trns;
 663   t->c = c;
 664   reg_get_name (dict, name, prefix);
 665   new_var = dict_create_var (dict, name, 0);
 666   assert (new_var != NULL);
 667   *v = new_var;
 668   add_transformation (ds, f, regression_trns_free, t);
 669   trns_index++;
 670 }
 671 static void
 672 subcommand_save (struct dataset *ds, int save, pspp_linreg_cache ** models)
 673 {
 674   pspp_linreg_cache **lc;
 675   int n_trns = 0;
 676   int i;
 677
 678   assert (models != NULL);
 679
 680   if (save)
 681     {
 682       /* Count the number of transformations we will need. */
 683       for (i = 0; i < REGRESSION_SV_count; i++)
 684         {
 685           if (cmd.a_save[i])
 686             {
 687               n_trns++;
 688             }
 689         }
 690       n_trns *= cmd.n_dependent;
 691
 692       for (lc = models; lc < models + cmd.n_dependent; lc++)
 693         {
 694           assert (*lc != NULL);
 695           assert ((*lc)->depvar != NULL);
 696           if (cmd.a_save[REGRESSION_SV_RESID])
 697             {
 698               reg_save_var (ds, "RES", regression_trns_resid_proc, *lc,
 699                             &(*lc)->resid, n_trns);
 700             }
 701           if (cmd.a_save[REGRESSION_SV_PRED])
 702             {
 703               reg_save_var (ds, "PRED", regression_trns_pred_proc, *lc,
 704                             &(*lc)->pred, n_trns);
 705             }
 706         }
 707     }
 708   else
 709     {
 710       for (lc = models; lc < models + cmd.n_dependent; lc++)
 711         {
 712           if (*lc != NULL)
 713             {
 714               pspp_linreg_cache_free (*lc);
 715             }
 716         }
 717     }
 718 }
 719
 720 int
 721 cmd_regression (struct lexer *lexer, struct dataset *ds)
 722 {
 723   struct casegrouper *grouper;
 724   struct casereader *group;
 725   pspp_linreg_cache **models;
 726   bool ok;
 727   size_t i;
 728
 729   if (!parse_regression (lexer, ds, &cmd, NULL))
 730     {
 731       return CMD_FAILURE;
 732     }
 733
 734   models = xnmalloc (cmd.n_dependent, sizeof *models);
 735   for (i = 0; i < cmd.n_dependent; i++)
 736     {
 737       models[i] = NULL;
 738     }
 739
 740   /* Data pass. */
 741   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 742   while (casegrouper_get_next_group (grouper, &group))
 743     run_regression (group, &cmd, ds, models);
 744   ok = casegrouper_destroy (grouper);
 745   ok = proc_commit (ds) && ok;
 746
 747   subcommand_save (ds, cmd.sbc_save, models);
 748   free (v_variables);
 749   free (models);
 750   free_regression (&cmd);
 751
 752   return ok ? CMD_SUCCESS : CMD_FAILURE;
 753 }
 754
 755 /*
 756   Is variable k the dependent variable?
 757  */
 758 static bool
 759 is_depvar (size_t k, const struct variable *v)
 760 {
 761   return v == v_variables[k];
 762 }
 763
 764 /* Parser for the variables sub command */
 765 static int
 766 regression_custom_variables (struct lexer *lexer, struct dataset *ds,
 767                              struct cmd_regression *cmd UNUSED,
 768                              void *aux UNUSED)
 769 {
 770   const struct dictionary *dict = dataset_dict (ds);
 771
 772   lex_match (lexer, '=');
 773
 774   if ((lex_token (lexer) != T_ID
 775        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 776       && lex_token (lexer) != T_ALL)
 777     return 2;
 778
 779
 780   if (!parse_variables_const
 781       (lexer, dict, &v_variables, &n_variables, PV_NONE))
 782     {
 783       free (v_variables);
 784       return 0;
 785     }
 786   assert (n_variables);
 787
 788   return 1;
 789 }
 790
 791 /* Identify the explanatory variables in v_variables.  Returns
 792    the number of independent variables. */
 793 static int
 794 identify_indep_vars (const struct variable **indep_vars,
 795                      const struct variable *depvar)
 796 {
 797   int n_indep_vars = 0;
 798   int i;
 799
 800   for (i = 0; i < n_variables; i++)
 801     if (!is_depvar (i, depvar))
 802       indep_vars[n_indep_vars++] = v_variables[i];
 803   if ((n_indep_vars < 1) && is_depvar (0, depvar))
 804     {
 805       /*
 806         There is only one independent variable, and it is the same
 807         as the dependent variable. Print a warning and continue.
 808        */
 809       msg (SE,
 810            gettext ("The dependent variable is equal to the independent variable."
 811                     "The least squares line is therefore Y=X."
 812                     "Standard errors and related statistics may be meaningless."));
 813       n_indep_vars = 1;
 814       indep_vars[0] = v_variables[0];
 815     }
 816   return n_indep_vars;
 817 }
 818
 819 /* Encode categorical variables.
 820    Returns number of valid cases. */
 821 static int
 822 prepare_categories (struct casereader *input,
 823                     const struct variable **vars, size_t n_vars,
 824                     struct moments_var *mom)
 825 {
 826   int n_data;
 827   struct ccase c;
 828   size_t i;
 829
 830   assert (vars != NULL);
 831   assert (mom != NULL);
 832
 833   for (i = 0; i < n_vars; i++)
 834     if (var_is_alpha (vars[i]))
 835       cat_stored_values_create (vars[i]);
 836
 837   n_data = 0;
 838   for (; casereader_read (input, &c); case_destroy (&c))
 839     {
 840       /*
 841          The second condition ensures the program will run even if
 842          there is only one variable to act as both explanatory and
 843          response.
 844        */
 845       for (i = 0; i < n_vars; i++)
 846         {
 847           const union value *val = case_data (&c, vars[i]);
 848           if (var_is_alpha (vars[i]))
 849             cat_value_update (vars[i], val);
 850           else
 851             moments1_add (mom[i].m, val->f, 1.0);
 852         }
 853       n_data++;
 854     }
 855   casereader_destroy (input);
 856
 857   return n_data;
 858 }
 859
 860 static void
 861 coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
 862 {
 863   c->coeff = xnmalloc (dm->m->size2, sizeof (*c->coeff));
 864   pspp_coeff_init (c->coeff, dm);
 865 }
 866
 867 /*
 868   Put the moments in the linreg cache.
 869  */
 870 static void
 871 compute_moments (pspp_linreg_cache * c, struct moments_var *mom,
 872                  struct design_matrix *dm, size_t n)
 873 {
 874   size_t i;
 875   size_t j;
 876   double weight;
 877   double mean;
 878   double variance;
 879   double skewness;
 880   double kurtosis;
 881   /*
 882      Scan the variable names in the columns of the design matrix.
 883      When we find the variable we need, insert its mean in the cache.
 884    */
 885   for (i = 0; i < dm->m->size2; i++)
 886     {
 887       for (j = 0; j < n; j++)
 888         {
 889           if (design_matrix_col_to_var (dm, i) == (mom + j)->v)
 890             {
 891               moments1_calculate ((mom + j)->m, &weight, &mean, &variance,
 892                                   &skewness, &kurtosis);
 893               gsl_vector_set (c->indep_means, i, mean);
 894               gsl_vector_set (c->indep_std, i, sqrt (variance));
 895             }
 896         }
 897     }
 898 }
 899
 900 static bool
 901 run_regression (struct casereader *input, struct cmd_regression *cmd,
 902                 struct dataset *ds, pspp_linreg_cache **models)
 903 {
 904   size_t i;
 905   int n_indep = 0;
 906   int k;
 907   struct ccase c;
 908   const struct variable **indep_vars;
 909   struct design_matrix *X;
 910   struct moments_var *mom;
 911   gsl_vector *Y;
 912
 913   pspp_linreg_opts lopts;
 914
 915   assert (models != NULL);
 916
 917   if (!casereader_peek (input, 0, &c))
 918     {
 919       casereader_destroy (input);
 920       return true;
 921     }
 922   output_split_file_values (ds, &c);
 923   case_destroy (&c);
 924
 925   if (!v_variables)
 926     {
 927       dict_get_vars (dataset_dict (ds), &v_variables, &n_variables, 0);
 928     }
 929
 930   for (i = 0; i < cmd->n_dependent; i++)
 931     {
 932       if (!var_is_numeric (cmd->v_dependent[i]))
 933         {
 934           msg (SE, _("Dependent variable must be numeric."));
 935           return false;
 936         }
 937     }
 938
 939   mom = xnmalloc (n_variables, sizeof (*mom));
 940   for (i = 0; i < n_variables; i++)
 941     {
 942       (mom + i)->m = moments1_create (MOMENT_VARIANCE);
 943       (mom + i)->v = v_variables[i];
 944     }
 945   lopts.get_depvar_mean_std = 1;
 946
 947   lopts.get_indep_mean_std = xnmalloc (n_variables, sizeof (int));
 948   indep_vars = xnmalloc (n_variables, sizeof *indep_vars);
 949
 950   for (k = 0; k < cmd->n_dependent; k++)
 951     {
 952       const struct variable *dep_var;
 953       struct casereader *reader;
 954       casenumber row;
 955       struct ccase c;
 956       size_t n_data;            /* Number of valid cases. */
 957
 958       dep_var = cmd->v_dependent[k];
 959       n_indep = identify_indep_vars (indep_vars, dep_var);
 960       reader = casereader_clone (input);
 961       reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
 962                                                  MV_ANY, NULL);
 963       reader = casereader_create_filter_missing (reader, &dep_var, 1,
 964                                                  MV_ANY, NULL);
 965       n_data = prepare_categories (casereader_clone (reader),
 966                                    indep_vars, n_indep, mom);
 967
 968       if ((n_data > 0) && (n_indep > 0))
 969         {
 970           Y = gsl_vector_alloc (n_data);
 971           X =
 972             design_matrix_create (n_indep,
 973                                   (const struct variable **) indep_vars,
 974                                   n_data);
 975           for (i = 0; i < X->m->size2; i++)
 976             {
 977               lopts.get_indep_mean_std[i] = 1;
 978             }
 979           models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
 980           models[k]->depvar = dep_var;
 981           /*
 982              For large data sets, use QR decomposition.
 983            */
 984           if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
 985             {
 986               models[k]->method = PSPP_LINREG_QR;
 987             }
 988
 989           /*
 990              The second pass fills the design matrix.
 991            */
 992           reader = casereader_create_counter (reader, &row, -1);
 993           for (; casereader_read (reader, &c); case_destroy (&c))
 994             {
 995               for (i = 0; i < n_indep; ++i)
 996                 {
 997                   const struct variable *v = indep_vars[i];
 998                   const union value *val = case_data (&c, v);
 999                   if (var_is_alpha (v))
1000                     design_matrix_set_categorical (X, row, v, val);
1001                   else
1002                     design_matrix_set_numeric (X, row, v, val);
1003                 }
1004               gsl_vector_set (Y, row, case_num (&c, dep_var));
1005             }
1006           /*
1007              Now that we know the number of coefficients, allocate space
1008              and store pointers to the variables that correspond to the
1009              coefficients.
1010            */
1011           coeff_init (models[k], X);
1012
1013           /*
1014              Find the least-squares estimates and other statistics.
1015            */
1016           pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
1017           compute_moments (models[k], mom, X, n_variables);
1018
1019           if (!taint_has_tainted_successor (casereader_get_taint (input)))
1020             {
1021               subcommand_statistics (cmd->a_statistics, models[k]);
1022             }
1023
1024           gsl_vector_free (Y);
1025           design_matrix_destroy (X);
1026         }
1027       else
1028         {
1029           msg (SE,
1030                gettext ("No valid data found. This command was skipped."));
1031         }
1032       casereader_destroy (reader);
1033     }
1034   for (i = 0; i < n_variables; i++)
1035     {
1036       moments1_destroy ((mom + i)->m);
1037     }
1038   free (mom);
1039   free (indep_vars);
1040   free (lopts.get_indep_mean_std);
1041   casereader_destroy (input);
1042
1043   return true;
1044 }
1045
1046 /*
1047   Local Variables:
1048   mode: c
1049   End:
1050 */