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       model_file = fh_parse (lexer, FH_REF_FILE);
 934       if (model_file == NULL)
 935         return 0;
 936     }
 937
 938   if (!lex_force_match (lexer, ')'))
 939     return 0;
 940
 941   return 1;
 942 }
 943
 944 int
 945 cmd_regression (struct lexer *lexer, struct dataset *ds)
 946 {
 947   struct casegrouper *grouper;
 948   struct casereader *group;
 949   pspp_linreg_cache **models;
 950   bool ok;
 951   size_t i;
 952
 953   if (!parse_regression (lexer, ds, &cmd, NULL))
 954     return CMD_FAILURE;
 955
 956   models = xnmalloc (cmd.n_dependent, sizeof *models);
 957   for (i = 0; i < cmd.n_dependent; i++)
 958     {
 959       models[i] = NULL;
 960     }
 961
 962   /* Data pass. */
 963   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 964   while (casegrouper_get_next_group (grouper, &group))
 965     run_regression (group, &cmd, ds, models);
 966   ok = casegrouper_destroy (grouper);
 967   ok = proc_commit (ds) && ok;
 968
 969   subcommand_save (ds, cmd.sbc_save, models);
 970   free (v_variables);
 971   free (models);
 972   free_regression (&cmd);
 973
 974   return ok ? CMD_SUCCESS : CMD_FAILURE;
 975 }
 976
 977 /*
 978   Is variable k the dependent variable?
 979  */
 980 static bool
 981 is_depvar (size_t k, const struct variable *v)
 982 {
 983   return v == v_variables[k];
 984 }
 985
 986 /* Parser for the variables sub command */
 987 static int
 988 regression_custom_variables (struct lexer *lexer, struct dataset *ds,
 989                              struct cmd_regression *cmd UNUSED,
 990                              void *aux UNUSED)
 991 {
 992   const struct dictionary *dict = dataset_dict (ds);
 993
 994   lex_match (lexer, '=');
 995
 996   if ((lex_token (lexer) != T_ID
 997        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 998       && lex_token (lexer) != T_ALL)
 999     return 2;
1000
1001
1002   if (!parse_variables_const
1003       (lexer, dict, &v_variables, &n_variables, PV_NONE))
1004     {
1005       free (v_variables);
1006       return 0;
1007     }
1008   assert (n_variables);
1009
1010   return 1;
1011 }
1012
1013 /* Identify the explanatory variables in v_variables.  Returns
1014    the number of independent variables. */
1015 static int
1016 identify_indep_vars (const struct variable **indep_vars,
1017                      const struct variable *depvar)
1018 {
1019   int n_indep_vars = 0;
1020   int i;
1021
1022   for (i = 0; i < n_variables; i++)
1023     if (!is_depvar (i, depvar))
1024       indep_vars[n_indep_vars++] = v_variables[i];
1025   if ((n_indep_vars < 2) && is_depvar (0, depvar))
1026     {
1027       /*
1028         There is only one independent variable, and it is the same
1029         as the dependent variable. Print a warning and continue.
1030        */
1031       msg (SE,
1032            gettext ("The dependent variable is equal to the independent variable."
1033                     "The least squares line is therefore Y=X."
1034                     "Standard errors and related statistics may be meaningless."));
1035       n_indep_vars = 1;
1036       indep_vars[0] = v_variables[0];
1037     }
1038   return n_indep_vars;
1039 }
1040
1041 /* Encode categorical variables.
1042    Returns number of valid cases. */
1043 static int
1044 prepare_categories (struct casereader *input,
1045                     const struct variable **vars, size_t n_vars,
1046                     struct moments_var *mom)
1047 {
1048   int n_data;
1049   struct ccase c;
1050   size_t i;
1051
1052   assert (vars != NULL);
1053   assert (mom != NULL);
1054
1055   for (i = 0; i < n_vars; i++)
1056     if (var_is_alpha (vars[i]))
1057       cat_stored_values_create (vars[i]);
1058
1059   n_data = 0;
1060   for (; casereader_read (input, &c); case_destroy (&c))
1061     {
1062       /*
1063          The second condition ensures the program will run even if
1064          there is only one variable to act as both explanatory and
1065          response.
1066        */
1067       for (i = 0; i < n_vars; i++)
1068         {
1069           const union value *val = case_data (&c, vars[i]);
1070           if (var_is_alpha (vars[i]))
1071             cat_value_update (vars[i], val);
1072           else
1073             moments1_add (mom[i].m, val->f, 1.0);
1074         }
1075       n_data++;
1076     }
1077   casereader_destroy (input);
1078
1079   return n_data;
1080 }
1081
1082 static void
1083 coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
1084 {
1085   c->coeff = xnmalloc (dm->m->size2 + 1, sizeof (*c->coeff));
1086   c->coeff[0] = xmalloc (sizeof (*(c->coeff[0])));      /* The first coefficient is the intercept. */
1087   c->coeff[0]->v_info = NULL;   /* Intercept has no associated variable. */
1088   pspp_coeff_init (c->coeff + 1, dm);
1089 }
1090
1091 /*
1092   Put the moments in the linreg cache.
1093  */
1094 static void
1095 compute_moments (pspp_linreg_cache * c, struct moments_var *mom,
1096                  struct design_matrix *dm, size_t n)
1097 {
1098   size_t i;
1099   size_t j;
1100   double weight;
1101   double mean;
1102   double variance;
1103   double skewness;
1104   double kurtosis;
1105   /*
1106      Scan the variable names in the columns of the design matrix.
1107      When we find the variable we need, insert its mean in the cache.
1108    */
1109   for (i = 0; i < dm->m->size2; i++)
1110     {
1111       for (j = 0; j < n; j++)
1112         {
1113           if (design_matrix_col_to_var (dm, i) == (mom + j)->v)
1114             {
1115               moments1_calculate ((mom + j)->m, &weight, &mean, &variance,
1116                                   &skewness, &kurtosis);
1117               gsl_vector_set (c->indep_means, i, mean);
1118               gsl_vector_set (c->indep_std, i, sqrt (variance));
1119             }
1120         }
1121     }
1122 }
1123
1124 static bool
1125 run_regression (struct casereader *input, struct cmd_regression *cmd,
1126                 struct dataset *ds, pspp_linreg_cache **models)
1127 {
1128   size_t i;
1129   int n_indep = 0;
1130   int k;
1131   struct ccase c;
1132   const struct variable **indep_vars;
1133   struct design_matrix *X;
1134   struct moments_var *mom;
1135   gsl_vector *Y;
1136
1137   pspp_linreg_opts lopts;
1138
1139   assert (models != NULL);
1140
1141   if (!casereader_peek (input, 0, &c))
1142     {
1143       casereader_destroy (input);
1144       return true;
1145     }
1146   output_split_file_values (ds, &c);
1147   case_destroy (&c);
1148
1149   if (!v_variables)
1150     {
1151       dict_get_vars (dataset_dict (ds), &v_variables, &n_variables,
1152                      1u << DC_SYSTEM);
1153     }
1154
1155   for (i = 0; i < cmd->n_dependent; i++)
1156     {
1157       if (!var_is_numeric (cmd->v_dependent[i]))
1158         {
1159           msg (SE, _("Dependent variable must be numeric."));
1160           return false;
1161         }
1162     }
1163
1164   mom = xnmalloc (n_variables, sizeof (*mom));
1165   for (i = 0; i < n_variables; i++)
1166     {
1167       (mom + i)->m = moments1_create (MOMENT_VARIANCE);
1168       (mom + i)->v = v_variables[i];
1169     }
1170   lopts.get_depvar_mean_std = 1;
1171
1172   lopts.get_indep_mean_std = xnmalloc (n_variables, sizeof (int));
1173   indep_vars = xnmalloc (n_variables, sizeof *indep_vars);
1174
1175   for (k = 0; k < cmd->n_dependent; k++)
1176     {
1177       const struct variable *dep_var;
1178       struct casereader *reader;
1179       casenumber row;
1180       struct ccase c;
1181       size_t n_data;            /* Number of valid cases. */
1182
1183       dep_var = cmd->v_dependent[k];
1184       n_indep = identify_indep_vars (indep_vars, dep_var);
1185       reader = casereader_clone (input);
1186       reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
1187                                                  MV_ANY, NULL);
1188       reader = casereader_create_filter_missing (reader, &dep_var, 1,
1189                                                  MV_ANY, NULL);
1190       n_data = prepare_categories (casereader_clone (reader),
1191                                    indep_vars, n_indep, mom);
1192
1193       if ((n_data > 0) && (n_indep > 0))
1194         {
1195           Y = gsl_vector_alloc (n_data);
1196           X =
1197             design_matrix_create (n_indep,
1198                                   (const struct variable **) indep_vars,
1199                                   n_data);
1200           for (i = 0; i < X->m->size2; i++)
1201             {
1202               lopts.get_indep_mean_std[i] = 1;
1203             }
1204           models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
1205           models[k]->depvar = dep_var;
1206           /*
1207              For large data sets, use QR decomposition.
1208            */
1209           if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
1210             {
1211               models[k]->method = PSPP_LINREG_QR;
1212             }
1213
1214           /*
1215              The second pass fills the design matrix.
1216            */
1217           reader = casereader_create_counter (reader, &row, -1);
1218           for (; casereader_read (reader, &c); case_destroy (&c))
1219             {
1220               for (i = 0; i < n_indep; ++i)
1221                 {
1222                   const struct variable *v = indep_vars[i];
1223                   const union value *val = case_data (&c, v);
1224                   if (var_is_alpha (v))
1225                     design_matrix_set_categorical (X, row, v, val);
1226                   else
1227                     design_matrix_set_numeric (X, row, v, val);
1228                 }
1229               gsl_vector_set (Y, row, case_num (&c, dep_var));
1230             }
1231           /*
1232              Now that we know the number of coefficients, allocate space
1233              and store pointers to the variables that correspond to the
1234              coefficients.
1235            */
1236           coeff_init (models[k], X);
1237
1238           /*
1239              Find the least-squares estimates and other statistics.
1240            */
1241           pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
1242           compute_moments (models[k], mom, X, n_variables);
1243
1244           if (!taint_has_tainted_successor (casereader_get_taint (input)))
1245             {
1246               subcommand_statistics (cmd->a_statistics, models[k]);
1247               subcommand_export (cmd->sbc_export, models[k]);
1248             }
1249
1250           gsl_vector_free (Y);
1251           design_matrix_destroy (X);
1252         }
1253       else
1254         {
1255           msg (SE,
1256                gettext ("No valid data found. This command was skipped."));
1257         }
1258       casereader_destroy (reader);
1259     }
1260   for (i = 0; i < n_variables; i++)
1261     {
1262       moments1_destroy ((mom + i)->m);
1263     }
1264   free (mom);
1265   free (indep_vars);
1266   free (lopts.get_indep_mean_std);
1267   casereader_destroy (input);
1268
1269   return true;
1270 }
1271
1272 /*
1273   Local Variables:
1274   mode: c
1275   End:
1276 */