src/language/stats/regression.q

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