src/language/stats/regression.q

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