src/language/stats/regression.q

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