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