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 "regression-export.h"
  26 #include <data/case.h>
  27 #include <data/casegrouper.h>
  28 #include <data/casereader.h>
  29 #include <data/category.h>
  30 #include <data/dictionary.h>
  31 #include <data/missing-values.h>
  32 #include <data/procedure.h>
  33 #include <data/transformations.h>
  34 #include <data/value-labels.h>
  35 #include <data/variable.h>
  36 #include <language/command.h>
  37 #include <language/dictionary/split-file.h>
  38 #include <language/data-io/file-handle.h>
  39 #include <language/lexer/lexer.h>
  40 #include <libpspp/compiler.h>
  41 #include <libpspp/message.h>
  42 #include <libpspp/taint.h>
  43 #include <math/design-matrix.h>
  44 #include <math/coefficient.h>
  45 #include <math/linreg/linreg.h>
  46 #include <math/moments.h>
  47 #include <output/table.h>
  48
  49 #include "xalloc.h"
  50
  51 #include "gettext.h"
  52 #define _(msgid) gettext (msgid)
  53
  54 #define REG_LARGE_DATA 1000
  55
  56 /* (headers) */
  57
  58 /* (specification)
  59    "REGRESSION" (regression_):
  60    *variables=custom;
  61    +statistics[st_]=r,
  62                     coeff,
  63                     anova,
  64                     outs,
  65                     zpp,
  66                     label,
  67                     sha,
  68                     ci,
  69                     bcov,
  70                     ses,
  71                     xtx,
  72                     collin,
  73                     tol,
  74                     selection,
  75                     f,
  76                     defaults,
  77                     all;
  78    export=custom;
  79    ^dependent=varlist;
  80    +save[sv_]=resid,pred;
  81    +method=enter.
  82 */
  83 /* (declarations) */
  84 /* (functions) */
  85 static struct cmd_regression cmd;
  86
  87 /*
  88   Moments for each of the variables used.
  89  */
  90 struct moments_var
  91 {
  92   struct moments1 *m;
  93   const struct variable *v;
  94 };
  95
  96 /*
  97   Transformations for saving predicted values
  98   and residuals, etc.
  99  */
 100 struct reg_trns
 101 {
 102   int n_trns;                   /* Number of transformations. */
 103   int trns_id;                  /* Which trns is this one? */
 104   pspp_linreg_cache *c;         /* Linear model for this trns. */
 105 };
 106 /*
 107   Variables used (both explanatory and response).
 108  */
 109 static const struct variable **v_variables;
 110
 111 /*
 112   Number of variables.
 113  */
 114 static size_t n_variables;
 115
 116 /*
 117   File where the model will be saved if the EXPORT subcommand
 118   is given.
 119  */
 120 static struct file_handle *model_file;
 121
 122 static bool run_regression (struct casereader *, struct cmd_regression *,
 123                             struct dataset *, pspp_linreg_cache **);
 124
 125 /*
 126    STATISTICS subcommand output functions.
 127  */
 128 static void reg_stats_r (pspp_linreg_cache *);
 129 static void reg_stats_coeff (pspp_linreg_cache *);
 130 static void reg_stats_anova (pspp_linreg_cache *);
 131 static void reg_stats_outs (pspp_linreg_cache *);
 132 static void reg_stats_zpp (pspp_linreg_cache *);
 133 static void reg_stats_label (pspp_linreg_cache *);
 134 static void reg_stats_sha (pspp_linreg_cache *);
 135 static void reg_stats_ci (pspp_linreg_cache *);
 136 static void reg_stats_f (pspp_linreg_cache *);
 137 static void reg_stats_bcov (pspp_linreg_cache *);
 138 static void reg_stats_ses (pspp_linreg_cache *);
 139 static void reg_stats_xtx (pspp_linreg_cache *);
 140 static void reg_stats_collin (pspp_linreg_cache *);
 141 static void reg_stats_tol (pspp_linreg_cache *);
 142 static void reg_stats_selection (pspp_linreg_cache *);
 143 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
 144                                        int, pspp_linreg_cache *);
 145
 146 static void
 147 reg_stats_r (pspp_linreg_cache * c)
 148 {
 149   struct tab_table *t;
 150   int n_rows = 2;
 151   int n_cols = 5;
 152   double rsq;
 153   double adjrsq;
 154   double std_error;
 155
 156   assert (c != NULL);
 157   rsq = c->ssm / c->sst;
 158   adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
 159   std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
 160   t = tab_create (n_cols, n_rows, 0);
 161   tab_dim (t, tab_natural_dimensions);
 162   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 163   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 164   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 165   tab_vline (t, TAL_0, 1, 0, 0);
 166
 167   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
 168   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
 169   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
 170   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
 171   tab_float (t, 1, 1, TAB_RIGHT, sqrt (rsq), 10, 2);
 172   tab_float (t, 2, 1, TAB_RIGHT, rsq, 10, 2);
 173   tab_float (t, 3, 1, TAB_RIGHT, adjrsq, 10, 2);
 174   tab_float (t, 4, 1, TAB_RIGHT, std_error, 10, 2);
 175   tab_title (t, _("Model Summary"));
 176   tab_submit (t);
 177 }
 178
 179 /*
 180   Table showing estimated regression coefficients.
 181  */
 182 static void
 183 reg_stats_coeff (pspp_linreg_cache * c)
 184 {
 185   size_t j;
 186   int n_cols = 7;
 187   int n_rows;
 188   double t_stat;
 189   double pval;
 190   double coeff;
 191   double std_err;
 192   double beta;
 193   const char *label;
 194   char *tmp;
 195   const struct variable *v;
 196   const union value *val;
 197   const char *val_s;
 198   struct tab_table *t;
 199
 200   assert (c != NULL);
 201   tmp = xnmalloc (MAX_STRING, sizeof (*tmp));
 202   n_rows = c->n_coeffs + 2;
 203
 204   t = tab_create (n_cols, n_rows, 0);
 205   tab_headers (t, 2, 0, 1, 0);
 206   tab_dim (t, tab_natural_dimensions);
 207   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 208   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 209   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 210   tab_vline (t, TAL_0, 1, 0, 0);
 211
 212   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
 213   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 214   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
 215   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
 216   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 217   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
 218   coeff = c->coeff[0]->estimate;
 219   tab_float (t, 2, 1, 0, coeff, 10, 2);
 220   std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
 221   tab_float (t, 3, 1, 0, std_err, 10, 2);
 222   beta = coeff / c->depvar_std;
 223   tab_float (t, 4, 1, 0, beta, 10, 2);
 224   t_stat = coeff / std_err;
 225   tab_float (t, 5, 1, 0, t_stat, 10, 2);
 226   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 227   tab_float (t, 6, 1, 0, pval, 10, 2);
 228   for (j = 1; j <= c->n_indeps; j++)
 229     {
 230       v = pspp_coeff_get_var (c->coeff[j], 0);
 231       label = var_to_string (v);
 232       /* Do not overwrite the variable's name. */
 233       strncpy (tmp, label, MAX_STRING);
 234       if (var_is_alpha (v))
 235         {
 236           /*
 237              Append the value associated with this coefficient.
 238              This makes sense only if we us the usual binary encoding
 239              for that value.
 240            */
 241
 242           val = pspp_coeff_get_value (c->coeff[j], v);
 243           val_s = var_get_value_name (v, val);
 244           strncat (tmp, val_s, MAX_STRING);
 245         }
 246
 247       tab_text (t, 1, j + 1, TAB_CENTER, tmp);
 248       /*
 249          Regression coefficients.
 250        */
 251       coeff = c->coeff[j]->estimate;
 252       tab_float (t, 2, j + 1, 0, coeff, 10, 2);
 253       /*
 254          Standard error of the coefficients.
 255        */
 256       std_err = sqrt (gsl_matrix_get (c->cov, j, j));
 257       tab_float (t, 3, j + 1, 0, std_err, 10, 2);
 258       /*
 259          'Standardized' coefficient, i.e., regression coefficient
 260          if all variables had unit variance.
 261        */
 262       beta = gsl_vector_get (c->indep_std, j);
 263       beta *= coeff / c->depvar_std;
 264       tab_float (t, 4, j + 1, 0, beta, 10, 2);
 265
 266       /*
 267          Test statistic for H0: coefficient is 0.
 268        */
 269       t_stat = coeff / std_err;
 270       tab_float (t, 5, j + 1, 0, t_stat, 10, 2);
 271       /*
 272          P values for the test statistic above.
 273        */
 274       pval =
 275         2 * gsl_cdf_tdist_Q (fabs (t_stat),
 276                              (double) (c->n_obs - c->n_coeffs));
 277       tab_float (t, 6, j + 1, 0, pval, 10, 2);
 278     }
 279   tab_title (t, _("Coefficients"));
 280   tab_submit (t);
 281   free (tmp);
 282 }
 283
 284 /*
 285   Display the ANOVA table.
 286  */
 287 static void
 288 reg_stats_anova (pspp_linreg_cache * c)
 289 {
 290   int n_cols = 7;
 291   int n_rows = 4;
 292   const double msm = c->ssm / c->dfm;
 293   const double mse = c->sse / c->dfe;
 294   const double F = msm / mse;
 295   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
 296
 297   struct tab_table *t;
 298
 299   assert (c != NULL);
 300   t = tab_create (n_cols, n_rows, 0);
 301   tab_headers (t, 2, 0, 1, 0);
 302   tab_dim (t, tab_natural_dimensions);
 303
 304   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 305
 306   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 307   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 308   tab_vline (t, TAL_0, 1, 0, 0);
 309
 310   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 311   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 312   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 313   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 314   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 315
 316   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
 317   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
 318   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
 319
 320   /* Sums of Squares */
 321   tab_float (t, 2, 1, 0, c->ssm, 10, 2);
 322   tab_float (t, 2, 3, 0, c->sst, 10, 2);
 323   tab_float (t, 2, 2, 0, c->sse, 10, 2);
 324
 325
 326   /* Degrees of freedom */
 327   tab_float (t, 3, 1, 0, c->dfm, 4, 0);
 328   tab_float (t, 3, 2, 0, c->dfe, 4, 0);
 329   tab_float (t, 3, 3, 0, c->dft, 4, 0);
 330
 331   /* Mean Squares */
 332
 333   tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
 334   tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
 335
 336   tab_float (t, 5, 1, 0, F, 8, 3);
 337
 338   tab_float (t, 6, 1, 0, pval, 8, 3);
 339
 340   tab_title (t, _("ANOVA"));
 341   tab_submit (t);
 342 }
 343 static void
 344 reg_stats_outs (pspp_linreg_cache * c)
 345 {
 346   assert (c != NULL);
 347 }
 348 static void
 349 reg_stats_zpp (pspp_linreg_cache * c)
 350 {
 351   assert (c != NULL);
 352 }
 353 static void
 354 reg_stats_label (pspp_linreg_cache * c)
 355 {
 356   assert (c != NULL);
 357 }
 358 static void
 359 reg_stats_sha (pspp_linreg_cache * c)
 360 {
 361   assert (c != NULL);
 362 }
 363 static void
 364 reg_stats_ci (pspp_linreg_cache * c)
 365 {
 366   assert (c != NULL);
 367 }
 368 static void
 369 reg_stats_f (pspp_linreg_cache * c)
 370 {
 371   assert (c != NULL);
 372 }
 373 static void
 374 reg_stats_bcov (pspp_linreg_cache * c)
 375 {
 376   int n_cols;
 377   int n_rows;
 378   int i;
 379   int k;
 380   int row;
 381   int col;
 382   const char *label;
 383   struct tab_table *t;
 384
 385   assert (c != NULL);
 386   n_cols = c->n_indeps + 1 + 2;
 387   n_rows = 2 * (c->n_indeps + 1);
 388   t = tab_create (n_cols, n_rows, 0);
 389   tab_headers (t, 2, 0, 1, 0);
 390   tab_dim (t, tab_natural_dimensions);
 391   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 392   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 393   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 394   tab_vline (t, TAL_0, 1, 0, 0);
 395   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 396   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 397   for (i = 1; i < c->n_coeffs; i++)
 398     {
 399       const struct variable *v = pspp_coeff_get_var (c->coeff[i], 0);
 400       label = var_to_string (v);
 401       tab_text (t, 2, i, TAB_CENTER, label);
 402       tab_text (t, i + 2, 0, TAB_CENTER, label);
 403       for (k = 1; k < c->n_coeffs; k++)
 404         {
 405           col = (i <= k) ? k : i;
 406           row = (i <= k) ? i : k;
 407           tab_float (t, k + 2, i, TAB_CENTER,
 408                      gsl_matrix_get (c->cov, row, col), 8, 3);
 409         }
 410     }
 411   tab_title (t, _("Coefficient Correlations"));
 412   tab_submit (t);
 413 }
 414 static void
 415 reg_stats_ses (pspp_linreg_cache * c)
 416 {
 417   assert (c != NULL);
 418 }
 419 static void
 420 reg_stats_xtx (pspp_linreg_cache * c)
 421 {
 422   assert (c != NULL);
 423 }
 424 static void
 425 reg_stats_collin (pspp_linreg_cache * c)
 426 {
 427   assert (c != NULL);
 428 }
 429 static void
 430 reg_stats_tol (pspp_linreg_cache * c)
 431 {
 432   assert (c != NULL);
 433 }
 434 static void
 435 reg_stats_selection (pspp_linreg_cache * c)
 436 {
 437   assert (c != NULL);
 438 }
 439
 440 static void
 441 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 442                            int keyword, pspp_linreg_cache * c)
 443 {
 444   if (keyword)
 445     {
 446       (*function) (c);
 447     }
 448 }
 449
 450 static void
 451 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 452 {
 453   /*
 454      The order here must match the order in which the STATISTICS
 455      keywords appear in the specification section above.
 456    */
 457   enum
 458   { r,
 459     coeff,
 460     anova,
 461     outs,
 462     zpp,
 463     label,
 464     sha,
 465     ci,
 466     bcov,
 467     ses,
 468     xtx,
 469     collin,
 470     tol,
 471     selection,
 472     f,
 473     defaults,
 474     all
 475   };
 476   int i;
 477   int d = 1;
 478
 479   if (keywords[all])
 480     {
 481       /*
 482          Set everything but F.
 483        */
 484       for (i = 0; i < f; i++)
 485         {
 486           keywords[i] = 1;
 487         }
 488     }
 489   else
 490     {
 491       for (i = 0; i < all; i++)
 492         {
 493           if (keywords[i])
 494             {
 495               d = 0;
 496             }
 497         }
 498       /*
 499          Default output: ANOVA table, parameter estimates,
 500          and statistics for variables not entered into model,
 501          if appropriate.
 502        */
 503       if (keywords[defaults] | d)
 504         {
 505           keywords[anova] = 1;
 506           keywords[outs] = 1;
 507           keywords[coeff] = 1;
 508           keywords[r] = 1;
 509         }
 510     }
 511   statistics_keyword_output (reg_stats_r, keywords[r], c);
 512   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 513   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 514   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 515   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 516   statistics_keyword_output (reg_stats_label, keywords[label], c);
 517   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 518   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 519   statistics_keyword_output (reg_stats_f, keywords[f], c);
 520   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 521   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 522   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 523   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 524   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 525   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 526 }
 527
 528 /*
 529   Free the transformation. Free its linear model if this
 530   transformation is the last one.
 531  */
 532 static bool
 533 regression_trns_free (void *t_)
 534 {
 535   bool result = true;
 536   struct reg_trns *t = t_;
 537
 538   if (t->trns_id == t->n_trns)
 539     {
 540       result = pspp_linreg_cache_free (t->c);
 541     }
 542   free (t);
 543
 544   return result;
 545 }
 546
 547 /*
 548   Gets the predicted values.
 549  */
 550 static int
 551 regression_trns_pred_proc (void *t_, struct ccase *c,
 552                            casenumber case_idx UNUSED)
 553 {
 554   size_t i;
 555   size_t n_vals;
 556   struct reg_trns *trns = t_;
 557   pspp_linreg_cache *model;
 558   union value *output = NULL;
 559   const union value **vals = NULL;
 560   const struct variable **vars = NULL;
 561
 562   assert (trns != NULL);
 563   model = trns->c;
 564   assert (model != NULL);
 565   assert (model->depvar != NULL);
 566   assert (model->pred != NULL);
 567
 568   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 569   n_vals = (*model->get_vars) (model, vars);
 570
 571   vals = xnmalloc (n_vals, sizeof (*vals));
 572   output = case_data_rw (c, model->pred);
 573   assert (output != NULL);
 574
 575   for (i = 0; i < n_vals; i++)
 576     {
 577       vals[i] = case_data (c, vars[i]);
 578     }
 579   output->f = (*model->predict) ((const struct variable **) vars,
 580                                  vals, model, n_vals);
 581   free (vals);
 582   free (vars);
 583   return TRNS_CONTINUE;
 584 }
 585
 586 /*
 587   Gets the residuals.
 588  */
 589 static int
 590 regression_trns_resid_proc (void *t_, struct ccase *c,
 591                             casenumber case_idx UNUSED)
 592 {
 593   size_t i;
 594   size_t n_vals;
 595   struct reg_trns *trns = t_;
 596   pspp_linreg_cache *model;
 597   union value *output = NULL;
 598   const union value **vals = NULL;
 599   const union value *obs = NULL;
 600   const struct variable **vars = NULL;
 601
 602   assert (trns != NULL);
 603   model = trns->c;
 604   assert (model != NULL);
 605   assert (model->depvar != NULL);
 606   assert (model->resid != NULL);
 607
 608   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 609   n_vals = (*model->get_vars) (model, vars);
 610
 611   vals = xnmalloc (n_vals, sizeof (*vals));
 612   output = case_data_rw (c, model->resid);
 613   assert (output != NULL);
 614
 615   for (i = 0; i < n_vals; i++)
 616     {
 617       vals[i] = case_data (c, vars[i]);
 618     }
 619   obs = case_data (c, model->depvar);
 620   output->f = (*model->residual) ((const struct variable **) vars,
 621                                   vals, obs, model, n_vals);
 622   free (vals);
 623   free (vars);
 624   return TRNS_CONTINUE;
 625 }
 626
 627 /*
 628    Returns false if NAME is a duplicate of any existing variable name.
 629 */
 630 static bool
 631 try_name (const struct dictionary *dict, const char *name)
 632 {
 633   if (dict_lookup_var (dict, name) != NULL)
 634     return false;
 635
 636   return true;
 637 }
 638
 639 static void
 640 reg_get_name (const struct dictionary *dict, char name[LONG_NAME_LEN],
 641               const char prefix[LONG_NAME_LEN])
 642 {
 643   int i = 1;
 644
 645   snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
 646   while (!try_name (dict, name))
 647     {
 648       i++;
 649       snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
 650     }
 651 }
 652
 653 static void
 654 reg_save_var (struct dataset *ds, const char *prefix, trns_proc_func * f,
 655               pspp_linreg_cache * c, struct variable **v, int n_trns)
 656 {
 657   struct dictionary *dict = dataset_dict (ds);
 658   static int trns_index = 1;
 659   char name[LONG_NAME_LEN];
 660   struct variable *new_var;
 661   struct reg_trns *t = NULL;
 662
 663   t = xmalloc (sizeof (*t));
 664   t->trns_id = trns_index;
 665   t->n_trns = n_trns;
 666   t->c = c;
 667   reg_get_name (dict, name, prefix);
 668   new_var = dict_create_var (dict, name, 0);
 669   assert (new_var != NULL);
 670   *v = new_var;
 671   add_transformation (ds, f, regression_trns_free, t);
 672   trns_index++;
 673 }
 674
 675 static void
 676 subcommand_save (struct dataset *ds, int save, pspp_linreg_cache ** models)
 677 {
 678   pspp_linreg_cache **lc;
 679   int n_trns = 0;
 680   int i;
 681
 682   assert (models != NULL);
 683
 684   if (save)
 685     {
 686       /* Count the number of transformations we will need. */
 687       for (i = 0; i < REGRESSION_SV_count; i++)
 688         {
 689           if (cmd.a_save[i])
 690             {
 691               n_trns++;
 692             }
 693         }
 694       n_trns *= cmd.n_dependent;
 695
 696       for (lc = models; lc < models + cmd.n_dependent; lc++)
 697         {
 698           assert (*lc != NULL);
 699           assert ((*lc)->depvar != NULL);
 700           if (cmd.a_save[REGRESSION_SV_RESID])
 701             {
 702               reg_save_var (ds, "RES", regression_trns_resid_proc, *lc,
 703                             &(*lc)->resid, n_trns);
 704             }
 705           if (cmd.a_save[REGRESSION_SV_PRED])
 706             {
 707               reg_save_var (ds, "PRED", regression_trns_pred_proc, *lc,
 708                             &(*lc)->pred, n_trns);
 709             }
 710         }
 711     }
 712   else
 713     {
 714       for (lc = models; lc < models + cmd.n_dependent; lc++)
 715         {
 716           if (*lc != NULL)
 717             {
 718               pspp_linreg_cache_free (*lc);
 719             }
 720         }
 721     }
 722 }
 723
 724 static int
 725 reg_inserted (const struct variable *v, struct variable **varlist, int n_vars)
 726 {
 727   int i;
 728
 729   for (i = 0; i < n_vars; i++)
 730     {
 731       if (v == varlist[i])
 732         {
 733           return 1;
 734         }
 735     }
 736   return 0;
 737 }
 738
 739 static void
 740 reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
 741 {
 742   int i;
 743   int n_vars = 0;
 744   struct variable **varlist;
 745
 746   fprintf (fp, "%s", reg_export_categorical_encode_1);
 747
 748   varlist = xnmalloc (c->n_indeps, sizeof (*varlist));
 749   for (i = 1; i < c->n_indeps; i++)     /* c->coeff[0] is the intercept. */
 750     {
 751       struct pspp_coeff *coeff = c->coeff[i];
 752       const struct variable *v = pspp_coeff_get_var (coeff, 0);
 753       if (var_is_alpha (v))
 754         {
 755           if (!reg_inserted (v, varlist, n_vars))
 756             {
 757               fprintf (fp, "struct pspp_reg_categorical_variable %s;\n\t",
 758                        var_get_name (v));
 759               varlist[n_vars] = (struct variable *) v;
 760               n_vars++;
 761             }
 762         }
 763     }
 764   fprintf (fp, "int n_vars = %d;\n\t", n_vars);
 765   fprintf (fp, "struct pspp_reg_categorical_variable *varlist[%d] = {",
 766            n_vars);
 767   for (i = 0; i < n_vars - 1; i++)
 768     {
 769       fprintf (fp, "&%s,\n\t\t", var_get_name (varlist[i]));
 770     }
 771   fprintf (fp, "&%s};\n\t", var_get_name (varlist[i]));
 772
 773   for (i = 0; i < n_vars; i++)
 774     {
 775       int n_categories = cat_get_n_categories (varlist[i]);
 776       int j;
 777
 778       fprintf (fp, "%s.name = \"%s\";\n\t",
 779                var_get_name (varlist[i]), var_get_name (varlist[i]));
 780       fprintf (fp, "%s.n_vals = %d;\n\t",
 781                var_get_name (varlist[i]), n_categories);
 782
 783       for (j = 0; j < n_categories; j++)
 784         {
 785           const union value *val = cat_subscript_to_value (j, varlist[i]);
 786           fprintf (fp, "%s.values[%d] = \"%s\";\n\t",
 787                    var_get_name (varlist[i]), j,
 788                    var_get_value_name (varlist[i], val));
 789         }
 790     }
 791   fprintf (fp, "%s", reg_export_categorical_encode_2);
 792 }
 793
 794 static void
 795 reg_print_depvars (FILE * fp, pspp_linreg_cache * c)
 796 {
 797   int i;
 798   struct pspp_coeff *coeff;
 799   const struct variable *v;
 800
 801   fprintf (fp, "char *model_depvars[%d] = {", c->n_indeps);
 802   for (i = 1; i < c->n_indeps; i++)
 803     {
 804       coeff = c->coeff[i];
 805       v = pspp_coeff_get_var (coeff, 0);
 806       fprintf (fp, "\"%s\",\n\t\t", var_get_name (v));
 807     }
 808   coeff = c->coeff[i];
 809   v = pspp_coeff_get_var (coeff, 0);
 810   fprintf (fp, "\"%s\"};\n\t", var_get_name (v));
 811 }
 812 static void
 813 reg_print_getvar (FILE * fp, pspp_linreg_cache * c)
 814 {
 815   fprintf (fp, "static int\npspp_reg_getvar (char *v_name)\n{\n\t");
 816   fprintf (fp, "int i;\n\tint n_vars = %d;\n\t", c->n_indeps);
 817   reg_print_depvars (fp, c);
 818   fprintf (fp, "for (i = 0; i < n_vars; i++)\n\t{\n\t\t");
 819   fprintf (fp,
 820            "if (strncmp (v_name, model_depvars[i], PSPP_REG_MAXLEN) == 0)\n\t\t{\n\t\t\t");
 821   fprintf (fp, "return i;\n\t\t}\n\t}\n}\n");
 822 }
 823 static int
 824 reg_has_categorical (pspp_linreg_cache * c)
 825 {
 826   int i;
 827   const struct variable *v;
 828
 829   for (i = 1; i < c->n_coeffs; i++)
 830     {
 831       v = pspp_coeff_get_var (c->coeff[i], 0);
 832       if (var_is_alpha (v))
 833         return 1;
 834     }
 835   return 0;
 836 }
 837
 838 static void
 839 subcommand_export (int export, pspp_linreg_cache * c)
 840 {
 841   FILE *fp;
 842   size_t i;
 843   size_t j;
 844   int n_quantiles = 100;
 845   double tmp;
 846   struct pspp_coeff *coeff;
 847
 848   if (export)
 849     {
 850       assert (c != NULL);
 851       assert (model_file != NULL);
 852       fp = fopen (fh_get_file_name (model_file), "w");
 853       assert (fp != NULL);
 854       fprintf (fp, "%s", reg_preamble);
 855       reg_print_getvar (fp, c);
 856       if (reg_has_categorical (c))
 857         {
 858           reg_print_categorical_encoding (fp, c);
 859         }
 860       fprintf (fp, "%s", reg_export_t_quantiles_1);
 861       for (i = 0; i < n_quantiles - 1; i++)
 862         {
 863           tmp = 0.5 + 0.005 * (double) i;
 864           fprintf (fp, "%.15e,\n\t\t",
 865                    gsl_cdf_tdist_Pinv (tmp, c->n_obs - c->n_indeps));
 866         }
 867       fprintf (fp, "%.15e};\n\t",
 868                gsl_cdf_tdist_Pinv (.9995, c->n_obs - c->n_indeps));
 869       fprintf (fp, "%s", reg_export_t_quantiles_2);
 870       fprintf (fp, "%s", reg_mean_cmt);
 871       fprintf (fp, "double\npspp_reg_estimate (const double *var_vals,");
 872       fprintf (fp, "const char *var_names[])\n{\n\t");
 873       fprintf (fp, "double model_coeffs[%d] = {", c->n_indeps);
 874       for (i = 1; i < c->n_indeps; i++)
 875         {
 876           coeff = c->coeff[i];
 877           fprintf (fp, "%.15e,\n\t\t", coeff->estimate);
 878         }
 879       coeff = c->coeff[i];
 880       fprintf (fp, "%.15e};\n\t", coeff->estimate);
 881       coeff = c->coeff[0];
 882       fprintf (fp, "double estimate = %.15e;\n\t", coeff->estimate);
 883       fprintf (fp, "int i;\n\tint j;\n\n\t");
 884       fprintf (fp, "for (i = 0; i < %d; i++)\n\t", c->n_indeps);
 885       fprintf (fp, "%s", reg_getvar);
 886       fprintf (fp, "const double cov[%d][%d] = {\n\t", c->n_coeffs,
 887                c->n_coeffs);
 888       for (i = 0; i < c->cov->size1 - 1; i++)
 889         {
 890           fprintf (fp, "{");
 891           for (j = 0; j < c->cov->size2 - 1; j++)
 892             {
 893               fprintf (fp, "%.15e, ", gsl_matrix_get (c->cov, i, j));
 894             }
 895           fprintf (fp, "%.15e},\n\t", gsl_matrix_get (c->cov, i, j));
 896         }
 897       fprintf (fp, "{");
 898       for (j = 0; j < c->cov->size2 - 1; j++)
 899         {
 900           fprintf (fp, "%.15e, ",
 901                    gsl_matrix_get (c->cov, c->cov->size1 - 1, j));
 902         }
 903       fprintf (fp, "%.15e}\n\t",
 904                gsl_matrix_get (c->cov, c->cov->size1 - 1, c->cov->size2 - 1));
 905       fprintf (fp, "};\n\tint n_vars = %d;\n\tint i;\n\tint j;\n\t",
 906                c->n_indeps);
 907       fprintf (fp, "double unshuffled_vals[%d];\n\t", c->n_indeps);
 908       fprintf (fp, "%s", reg_variance);
 909       fprintf (fp, "%s", reg_export_confidence_interval);
 910       tmp = c->mse * c->mse;
 911       fprintf (fp, "%s %.15e", reg_export_prediction_interval_1, tmp);
 912       fprintf (fp, "%s %.15e", reg_export_prediction_interval_2, tmp);
 913       fprintf (fp, "%s", reg_export_prediction_interval_3);
 914       fclose (fp);
 915       fp = fopen ("pspp_model_reg.h", "w");
 916       fprintf (fp, "%s", reg_header);
 917       fclose (fp);
 918     }
 919 }
 920
 921 static int
 922 regression_custom_export (struct lexer *lexer, struct dataset *ds UNUSED,
 923                           struct cmd_regression *cmd UNUSED, void *aux UNUSED)
 924 {
 925   /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
 926   if (!lex_force_match (lexer, '('))
 927     return 0;
 928
 929   if (lex_match (lexer, '*'))
 930     model_file = NULL;
 931   else
 932     {
 933       fh_unref (model_file);
 934       model_file = fh_parse (lexer, FH_REF_FILE);
 935       if (model_file == NULL)
 936         return 0;
 937     }
 938
 939   if (!lex_force_match (lexer, ')'))
 940     return 0;
 941
 942   return 1;
 943 }
 944
 945 int
 946 cmd_regression (struct lexer *lexer, struct dataset *ds)
 947 {
 948   struct casegrouper *grouper;
 949   struct casereader *group;
 950   pspp_linreg_cache **models;
 951   bool ok;
 952   size_t i;
 953
 954   model_file = NULL;
 955   if (!parse_regression (lexer, ds, &cmd, NULL))
 956     {
 957       fh_unref (model_file);
 958       return CMD_FAILURE;
 959     }
 960
 961   models = xnmalloc (cmd.n_dependent, sizeof *models);
 962   for (i = 0; i < cmd.n_dependent; i++)
 963     {
 964       models[i] = NULL;
 965     }
 966
 967   /* Data pass. */
 968   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 969   while (casegrouper_get_next_group (grouper, &group))
 970     run_regression (group, &cmd, ds, models);
 971   ok = casegrouper_destroy (grouper);
 972   ok = proc_commit (ds) && ok;
 973
 974   subcommand_save (ds, cmd.sbc_save, models);
 975   free (v_variables);
 976   free (models);
 977   free_regression (&cmd);
 978   fh_unref (model_file);
 979
 980   return ok ? CMD_SUCCESS : CMD_FAILURE;
 981 }
 982
 983 /*
 984   Is variable k the dependent variable?
 985  */
 986 static bool
 987 is_depvar (size_t k, const struct variable *v)
 988 {
 989   return v == v_variables[k];
 990 }
 991
 992 /* Parser for the variables sub command */
 993 static int
 994 regression_custom_variables (struct lexer *lexer, struct dataset *ds,
 995                              struct cmd_regression *cmd UNUSED,
 996                              void *aux UNUSED)
 997 {
 998   const struct dictionary *dict = dataset_dict (ds);
 999
1000   lex_match (lexer, '=');
1001
1002   if ((lex_token (lexer) != T_ID
1003        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
1004       && lex_token (lexer) != T_ALL)
1005     return 2;
1006
1007
1008   if (!parse_variables_const
1009       (lexer, dict, &v_variables, &n_variables, PV_NONE))
1010     {
1011       free (v_variables);
1012       return 0;
1013     }
1014   assert (n_variables);
1015
1016   return 1;
1017 }
1018
1019 /* Identify the explanatory variables in v_variables.  Returns
1020    the number of independent variables. */
1021 static int
1022 identify_indep_vars (const struct variable **indep_vars,
1023                      const struct variable *depvar)
1024 {
1025   int n_indep_vars = 0;
1026   int i;
1027
1028   for (i = 0; i < n_variables; i++)
1029     if (!is_depvar (i, depvar))
1030       indep_vars[n_indep_vars++] = v_variables[i];
1031   if ((n_indep_vars < 2) && is_depvar (0, depvar))
1032     {
1033       /*
1034         There is only one independent variable, and it is the same
1035         as the dependent variable. Print a warning and continue.
1036        */
1037       msg (SE,
1038            gettext ("The dependent variable is equal to the independent variable."
1039                     "The least squares line is therefore Y=X."
1040                     "Standard errors and related statistics may be meaningless."));
1041       n_indep_vars = 1;
1042       indep_vars[0] = v_variables[0];
1043     }
1044   return n_indep_vars;
1045 }
1046
1047 /* Encode categorical variables.
1048    Returns number of valid cases. */
1049 static int
1050 prepare_categories (struct casereader *input,
1051                     const struct variable **vars, size_t n_vars,
1052                     struct moments_var *mom)
1053 {
1054   int n_data;
1055   struct ccase c;
1056   size_t i;
1057
1058   assert (vars != NULL);
1059   assert (mom != NULL);
1060
1061   for (i = 0; i < n_vars; i++)
1062     if (var_is_alpha (vars[i]))
1063       cat_stored_values_create (vars[i]);
1064
1065   n_data = 0;
1066   for (; casereader_read (input, &c); case_destroy (&c))
1067     {
1068       /*
1069          The second condition ensures the program will run even if
1070          there is only one variable to act as both explanatory and
1071          response.
1072        */
1073       for (i = 0; i < n_vars; i++)
1074         {
1075           const union value *val = case_data (&c, vars[i]);
1076           if (var_is_alpha (vars[i]))
1077             cat_value_update (vars[i], val);
1078           else
1079             moments1_add (mom[i].m, val->f, 1.0);
1080         }
1081       n_data++;
1082     }
1083   casereader_destroy (input);
1084
1085   return n_data;
1086 }
1087
1088 static void
1089 coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
1090 {
1091   c->coeff = xnmalloc (dm->m->size2 + 1, sizeof (*c->coeff));
1092   c->coeff[0] = xmalloc (sizeof (*(c->coeff[0])));      /* The first coefficient is the intercept. */
1093   c->coeff[0]->v_info = NULL;   /* Intercept has no associated variable. */
1094   pspp_coeff_init (c->coeff + 1, dm);
1095 }
1096
1097 /*
1098   Put the moments in the linreg cache.
1099  */
1100 static void
1101 compute_moments (pspp_linreg_cache * c, struct moments_var *mom,
1102                  struct design_matrix *dm, size_t n)
1103 {
1104   size_t i;
1105   size_t j;
1106   double weight;
1107   double mean;
1108   double variance;
1109   double skewness;
1110   double kurtosis;
1111   /*
1112      Scan the variable names in the columns of the design matrix.
1113      When we find the variable we need, insert its mean in the cache.
1114    */
1115   for (i = 0; i < dm->m->size2; i++)
1116     {
1117       for (j = 0; j < n; j++)
1118         {
1119           if (design_matrix_col_to_var (dm, i) == (mom + j)->v)
1120             {
1121               moments1_calculate ((mom + j)->m, &weight, &mean, &variance,
1122                                   &skewness, &kurtosis);
1123               gsl_vector_set (c->indep_means, i, mean);
1124               gsl_vector_set (c->indep_std, i, sqrt (variance));
1125             }
1126         }
1127     }
1128 }
1129
1130 static bool
1131 run_regression (struct casereader *input, struct cmd_regression *cmd,
1132                 struct dataset *ds, pspp_linreg_cache **models)
1133 {
1134   size_t i;
1135   int n_indep = 0;
1136   int k;
1137   struct ccase c;
1138   const struct variable **indep_vars;
1139   struct design_matrix *X;
1140   struct moments_var *mom;
1141   gsl_vector *Y;
1142
1143   pspp_linreg_opts lopts;
1144
1145   assert (models != NULL);
1146
1147   if (!casereader_peek (input, 0, &c))
1148     {
1149       casereader_destroy (input);
1150       return true;
1151     }
1152   output_split_file_values (ds, &c);
1153   case_destroy (&c);
1154
1155   if (!v_variables)
1156     {
1157       dict_get_vars (dataset_dict (ds), &v_variables, &n_variables,
1158                      1u << DC_SYSTEM);
1159     }
1160
1161   for (i = 0; i < cmd->n_dependent; i++)
1162     {
1163       if (!var_is_numeric (cmd->v_dependent[i]))
1164         {
1165           msg (SE, _("Dependent variable must be numeric."));
1166           return false;
1167         }
1168     }
1169
1170   mom = xnmalloc (n_variables, sizeof (*mom));
1171   for (i = 0; i < n_variables; i++)
1172     {
1173       (mom + i)->m = moments1_create (MOMENT_VARIANCE);
1174       (mom + i)->v = v_variables[i];
1175     }
1176   lopts.get_depvar_mean_std = 1;
1177
1178   lopts.get_indep_mean_std = xnmalloc (n_variables, sizeof (int));
1179   indep_vars = xnmalloc (n_variables, sizeof *indep_vars);
1180
1181   for (k = 0; k < cmd->n_dependent; k++)
1182     {
1183       const struct variable *dep_var;
1184       struct casereader *reader;
1185       casenumber row;
1186       struct ccase c;
1187       size_t n_data;            /* Number of valid cases. */
1188
1189       dep_var = cmd->v_dependent[k];
1190       n_indep = identify_indep_vars (indep_vars, dep_var);
1191       reader = casereader_clone (input);
1192       reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
1193                                                  MV_ANY, NULL);
1194       reader = casereader_create_filter_missing (reader, &dep_var, 1,
1195                                                  MV_ANY, NULL);
1196       n_data = prepare_categories (casereader_clone (reader),
1197                                    indep_vars, n_indep, mom);
1198
1199       if ((n_data > 0) && (n_indep > 0))
1200         {
1201           Y = gsl_vector_alloc (n_data);
1202           X =
1203             design_matrix_create (n_indep,
1204                                   (const struct variable **) indep_vars,
1205                                   n_data);
1206           for (i = 0; i < X->m->size2; i++)
1207             {
1208               lopts.get_indep_mean_std[i] = 1;
1209             }
1210           models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
1211           models[k]->depvar = dep_var;
1212           /*
1213              For large data sets, use QR decomposition.
1214            */
1215           if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
1216             {
1217               models[k]->method = PSPP_LINREG_QR;
1218             }
1219
1220           /*
1221              The second pass fills the design matrix.
1222            */
1223           reader = casereader_create_counter (reader, &row, -1);
1224           for (; casereader_read (reader, &c); case_destroy (&c))
1225             {
1226               for (i = 0; i < n_indep; ++i)
1227                 {
1228                   const struct variable *v = indep_vars[i];
1229                   const union value *val = case_data (&c, v);
1230                   if (var_is_alpha (v))
1231                     design_matrix_set_categorical (X, row, v, val);
1232                   else
1233                     design_matrix_set_numeric (X, row, v, val);
1234                 }
1235               gsl_vector_set (Y, row, case_num (&c, dep_var));
1236             }
1237           /*
1238              Now that we know the number of coefficients, allocate space
1239              and store pointers to the variables that correspond to the
1240              coefficients.
1241            */
1242           coeff_init (models[k], X);
1243
1244           /*
1245              Find the least-squares estimates and other statistics.
1246            */
1247           pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
1248           compute_moments (models[k], mom, X, n_variables);
1249
1250           if (!taint_has_tainted_successor (casereader_get_taint (input)))
1251             {
1252               subcommand_statistics (cmd->a_statistics, models[k]);
1253               subcommand_export (cmd->sbc_export, models[k]);
1254             }
1255
1256           gsl_vector_free (Y);
1257           design_matrix_destroy (X);
1258         }
1259       else
1260         {
1261           msg (SE,
1262                gettext ("No valid data found. This command was skipped."));
1263         }
1264       casereader_destroy (reader);
1265     }
1266   for (i = 0; i < n_variables; i++)
1267     {
1268       moments1_destroy ((mom + i)->m);
1269     }
1270   free (mom);
1271   free (indep_vars);
1272   free (lopts.get_indep_mean_std);
1273   casereader_destroy (input);
1274
1275   return true;
1276 }
1277
1278 /*
1279   Local Variables:
1280   mode: c
1281   End:
1282 */