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/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 static 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 static struct file_handle *model_file;
 114
 115 /*
 116   Return value for the procedure.
 117  */
 118 static 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_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_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_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,
 548                            casenum_t case_idx UNUSED)
 549 {
 550   size_t i;
 551   size_t n_vals;
 552   struct reg_trns *trns = t_;
 553   pspp_linreg_cache *model;
 554   union value *output = NULL;
 555   const union value **vals = NULL;
 556   struct variable **vars = NULL;
 557
 558   assert (trns != NULL);
 559   model = trns->c;
 560   assert (model != NULL);
 561   assert (model->depvar != NULL);
 562   assert (model->pred != NULL);
 563
 564   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 565   n_vals = (*model->get_vars) (model, vars);
 566
 567   vals = xnmalloc (n_vals, sizeof (*vals));
 568   output = case_data_rw (c, model->pred->fv);
 569   assert (output != NULL);
 570
 571   for (i = 0; i < n_vals; i++)
 572     {
 573       vals[i] = case_data (c, vars[i]->fv);
 574     }
 575   output->f = (*model->predict) ((const struct variable **) vars,
 576                                  vals, model, n_vals);
 577   free (vals);
 578   free (vars);
 579   return TRNS_CONTINUE;
 580 }
 581
 582 /*
 583   Gets the residuals.
 584  */
 585 static int
 586 regression_trns_resid_proc (void *t_, struct ccase *c,
 587                             casenum_t case_idx UNUSED)
 588 {
 589   size_t i;
 590   size_t n_vals;
 591   struct reg_trns *trns = t_;
 592   pspp_linreg_cache *model;
 593   union value *output = NULL;
 594   const union value **vals = NULL;
 595   const union value *obs = NULL;
 596   struct variable **vars = NULL;
 597
 598   assert (trns != NULL);
 599   model = trns->c;
 600   assert (model != NULL);
 601   assert (model->depvar != NULL);
 602   assert (model->resid != NULL);
 603
 604   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 605   n_vals = (*model->get_vars) (model, vars);
 606
 607   vals = xnmalloc (n_vals, sizeof (*vals));
 608   output = case_data_rw (c, model->resid->fv);
 609   assert (output != NULL);
 610
 611   for (i = 0; i < n_vals; i++)
 612     {
 613       vals[i] = case_data (c, vars[i]->fv);
 614     }
 615   obs = case_data (c, model->depvar->fv);
 616   output->f = (*model->residual) ((const struct variable **) vars,
 617                                   vals, obs, model, n_vals);
 618   free (vals);
 619   free (vars);
 620   return TRNS_CONTINUE;
 621 }
 622
 623 /*
 624    Returns 0 if NAME is a duplicate of any existing variable name.
 625 */
 626 static int
 627 try_name (char *name)
 628 {
 629   if (dict_lookup_var (dataset_dict (current_dataset), name) != NULL)
 630     return 0;
 631
 632   return 1;
 633 }
 634 static void
 635 reg_get_name (char name[LONG_NAME_LEN], const char prefix[LONG_NAME_LEN])
 636 {
 637   int i = 1;
 638
 639   snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
 640   while (!try_name (name))
 641     {
 642       i++;
 643       snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
 644     }
 645 }
 646 static void
 647 reg_save_var (const char *prefix, trns_proc_func * f,
 648               pspp_linreg_cache * c, struct variable **v, int n_trns)
 649 {
 650   static int trns_index = 1;
 651   char name[LONG_NAME_LEN];
 652   struct variable *new_var;
 653   struct reg_trns *t = NULL;
 654
 655   t = xmalloc (sizeof (*t));
 656   t->trns_id = trns_index;
 657   t->n_trns = n_trns;
 658   t->c = c;
 659   reg_get_name (name, prefix);
 660   new_var = dict_create_var (dataset_dict (current_dataset), name, 0);
 661   assert (new_var != NULL);
 662   *v = new_var;
 663   add_transformation (current_dataset, f, regression_trns_free, t);
 664   trns_index++;
 665 }
 666 static void
 667 subcommand_save (int save, pspp_linreg_cache ** models)
 668 {
 669   pspp_linreg_cache **lc;
 670   int n_trns = 0;
 671   int i;
 672
 673   assert (models != NULL);
 674
 675   if (save)
 676     {
 677       /* Count the number of transformations we will need. */
 678       for (i = 0; i < REGRESSION_SV_count; i++)
 679         {
 680           if (cmd.a_save[i])
 681             {
 682               n_trns++;
 683             }
 684         }
 685       n_trns *= cmd.n_dependent;
 686
 687       for (lc = models; lc < models + cmd.n_dependent; lc++)
 688         {
 689           assert (*lc != NULL);
 690           assert ((*lc)->depvar != NULL);
 691           if (cmd.a_save[REGRESSION_SV_RESID])
 692             {
 693               reg_save_var ("RES", regression_trns_resid_proc, *lc,
 694                             &(*lc)->resid, n_trns);
 695             }
 696           if (cmd.a_save[REGRESSION_SV_PRED])
 697             {
 698               reg_save_var ("PRED", regression_trns_pred_proc, *lc,
 699                             &(*lc)->pred, n_trns);
 700             }
 701         }
 702     }
 703   else
 704     {
 705       for (lc = models; lc < models + cmd.n_dependent; lc++)
 706         {
 707           assert (*lc != NULL);
 708           pspp_linreg_cache_free (*lc);
 709         }
 710     }
 711 }
 712 static int
 713 reg_inserted (const struct variable *v, struct variable **varlist, int n_vars)
 714 {
 715   int i;
 716
 717   for (i = 0; i < n_vars; i++)
 718     {
 719       if (v->index == varlist[i]->index)
 720         {
 721           return 1;
 722         }
 723     }
 724   return 0;
 725 }
 726 static void
 727 reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
 728 {
 729   int i;
 730   size_t j;
 731   int n_vars = 0;
 732   struct variable **varlist;
 733   struct pspp_coeff *coeff;
 734   const struct variable *v;
 735   union value *val;
 736
 737   fprintf (fp, "%s", reg_export_categorical_encode_1);
 738
 739   varlist = xnmalloc (c->n_indeps, sizeof (*varlist));
 740   for (i = 1; i < c->n_indeps; i++)     /* c->coeff[0] is the intercept. */
 741     {
 742       coeff = c->coeff[i];
 743       v = pspp_coeff_get_var (coeff, 0);
 744       if (v->type == ALPHA)
 745         {
 746           if (!reg_inserted (v, varlist, n_vars))
 747             {
 748               fprintf (fp, "struct pspp_reg_categorical_variable %s;\n\t",
 749                        v->name);
 750               varlist[n_vars] = (struct variable *) v;
 751               n_vars++;
 752             }
 753         }
 754     }
 755   fprintf (fp, "int n_vars = %d;\n\t", n_vars);
 756   fprintf (fp, "struct pspp_reg_categorical_variable *varlist[%d] = {",
 757            n_vars);
 758   for (i = 0; i < n_vars - 1; i++)
 759     {
 760       fprintf (fp, "&%s,\n\t\t", varlist[i]->name);
 761     }
 762   fprintf (fp, "&%s};\n\t", varlist[i]->name);
 763
 764   for (i = 0; i < n_vars; i++)
 765     {
 766       coeff = c->coeff[i];
 767       fprintf (fp, "%s.name = \"%s\";\n\t", varlist[i]->name,
 768                varlist[i]->name);
 769       fprintf (fp, "%s.n_vals = %d;\n\t", varlist[i]->name,
 770                varlist[i]->obs_vals->n_categories);
 771
 772       for (j = 0; j < varlist[i]->obs_vals->n_categories; j++)
 773         {
 774           val = cat_subscript_to_value ((const size_t) j, varlist[i]);
 775           fprintf (fp, "%s.values[%d] = \"%s\";\n\t", varlist[i]->name, j,
 776                    value_to_string (val, varlist[i]));
 777         }
 778     }
 779   fprintf (fp, "%s", reg_export_categorical_encode_2);
 780 }
 781
 782 static void
 783 reg_print_depvars (FILE * fp, pspp_linreg_cache * c)
 784 {
 785   int i;
 786   struct pspp_coeff *coeff;
 787   const struct variable *v;
 788
 789   fprintf (fp, "char *model_depvars[%d] = {", c->n_indeps);
 790   for (i = 1; i < c->n_indeps; i++)
 791     {
 792       coeff = c->coeff[i];
 793       v = pspp_coeff_get_var (coeff, 0);
 794       fprintf (fp, "\"%s\",\n\t\t", v->name);
 795     }
 796   coeff = c->coeff[i];
 797   v = pspp_coeff_get_var (coeff, 0);
 798   fprintf (fp, "\"%s\"};\n\t", v->name);
 799 }
 800 static void
 801 reg_print_getvar (FILE * fp, pspp_linreg_cache * c)
 802 {
 803   fprintf (fp, "static int\npspp_reg_getvar (char *v_name)\n{\n\t");
 804   fprintf (fp, "int i;\n\tint n_vars = %d;\n\t", c->n_indeps);
 805   reg_print_depvars (fp, c);
 806   fprintf (fp, "for (i = 0; i < n_vars; i++)\n\t{\n\t\t");
 807   fprintf (fp,
 808            "if (strncmp (v_name, model_depvars[i], PSPP_REG_MAXLEN) == 0)\n\t\t{\n\t\t\t");
 809   fprintf (fp, "return i;\n\t\t}\n\t}\n}\n");
 810 }
 811 static int
 812 reg_has_categorical (pspp_linreg_cache * c)
 813 {
 814   int i;
 815   const struct variable *v;
 816
 817   for (i = 1; i < c->n_coeffs; i++)
 818     {
 819       v = pspp_coeff_get_var (c->coeff[i], 0);
 820       if (v->type == ALPHA)
 821         {
 822           return 1;
 823         }
 824     }
 825   return 0;
 826 }
 827
 828 static void
 829 subcommand_export (int export, pspp_linreg_cache * c)
 830 {
 831   FILE *fp;
 832   size_t i;
 833   size_t j;
 834   int n_quantiles = 100;
 835   double tmp;
 836   struct pspp_coeff *coeff;
 837
 838   if (export)
 839     {
 840       assert (c != NULL);
 841       assert (model_file != NULL);
 842       fp = fopen (fh_get_file_name (model_file), "w");
 843       assert (fp != NULL);
 844       fprintf (fp, "%s", reg_preamble);
 845       reg_print_getvar (fp, c);
 846       if (reg_has_categorical (c))
 847         {
 848           reg_print_categorical_encoding (fp, c);
 849         }
 850       fprintf (fp, "%s", reg_export_t_quantiles_1);
 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, void *aux 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, NULL))
 936     return CMD_FAILURE;
 937
 938   models = xnmalloc (cmd.n_dependent, sizeof *models);
 939   if (!multipass_procedure_with_splits (current_dataset, 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 bool
 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 true;
 959
 960   return false;
 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                              void *aux UNUSED)
1001 {
1002
1003   lex_match ('=');
1004
1005   if ((token != T_ID || dict_lookup_var (dataset_dict (current_dataset), tokid) == NULL)
1006       && token != T_ALL)
1007     return 2;
1008
1009
1010   if (!parse_variables (dataset_dict (current_dataset), &v_variables, &n_variables, PV_NONE))
1011     {
1012       free (v_variables);
1013       return 0;
1014     }
1015   assert (n_variables);
1016
1017   return 1;
1018 }
1019
1020 /*
1021   Count the explanatory variables. The user may or may
1022   not have specified a response variable in the syntax.
1023  */
1024 static int
1025 get_n_indep (const struct variable *v)
1026 {
1027   int result;
1028   int i = 0;
1029
1030   result = n_variables;
1031   while (i < n_variables)
1032     {
1033       if (is_depvar (i, v))
1034         {
1035           result--;
1036           i = n_variables;
1037         }
1038       i++;
1039     }
1040   return result;
1041 }
1042
1043 /*
1044   Read from the active file. Identify the explanatory variables in
1045   v_variables. Encode categorical variables. Drop cases with missing
1046   values.
1047 */
1048 static int
1049 prepare_data (int n_data, int is_missing_case[],
1050               struct variable **indep_vars,
1051               struct variable *depvar, const struct casefile *cf)
1052 {
1053   int i;
1054   int j;
1055
1056   assert (indep_vars != NULL);
1057   j = 0;
1058   for (i = 0; i < n_variables; i++)
1059     {
1060       if (!is_depvar (i, depvar))
1061         {
1062           indep_vars[j] = v_variables[i];
1063           j++;
1064           if (v_variables[i]->type == ALPHA)
1065             {
1066               /* Make a place to hold the binary vectors
1067                  corresponding to this variable's values. */
1068               cat_stored_values_create (v_variables[i]);
1069             }
1070           n_data =
1071             mark_missing_cases (cf, v_variables[i], is_missing_case, n_data);
1072         }
1073     }
1074   /*
1075      Mark missing cases for the dependent variable.
1076    */
1077   n_data = mark_missing_cases (cf, depvar, is_missing_case, n_data);
1078
1079   return n_data;
1080 }
1081 static void
1082 coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
1083 {
1084   c->coeff = xnmalloc (dm->m->size2 + 1, sizeof (*c->coeff));
1085   c->coeff[0] = xmalloc (sizeof (*(c->coeff[0])));      /* The first coefficient is the intercept. */
1086   c->coeff[0]->v_info = NULL;   /* Intercept has no associated variable. */
1087   pspp_coeff_init (c->coeff + 1, dm);
1088 }
1089
1090 static bool
1091 run_regression (const struct ccase *first,
1092                 const struct casefile *cf, void *cmd_ UNUSED)
1093 {
1094   size_t i;
1095   size_t n_data = 0;            /* Number of valide cases. */
1096   size_t n_cases;               /* Number of cases. */
1097   size_t row;
1098   size_t case_num;
1099   int n_indep = 0;
1100   int k;
1101   /*
1102      Keep track of the missing cases.
1103    */
1104   int *is_missing_case;
1105   const union value *val;
1106   struct casereader *r;
1107   struct ccase c;
1108   struct variable **indep_vars;
1109   struct design_matrix *X;
1110   gsl_vector *Y;
1111
1112   pspp_linreg_opts lopts;
1113
1114   assert (models != NULL);
1115
1116   output_split_file_values (first);
1117
1118   if (!v_variables)
1119     {
1120       dict_get_vars (dataset_dict (current_dataset), &v_variables, &n_variables,
1121                      1u << DC_SYSTEM);
1122     }
1123
1124   n_cases = casefile_get_case_cnt (cf);
1125
1126   for (i = 0; i < cmd.n_dependent; i++)
1127     {
1128       if (cmd.v_dependent[i]->type != NUMERIC)
1129         {
1130           msg (SE, gettext ("Dependent variable must be numeric."));
1131           pspp_reg_rc = CMD_FAILURE;
1132           return true;
1133         }
1134     }
1135
1136   is_missing_case = xnmalloc (n_cases, sizeof (*is_missing_case));
1137
1138   lopts.get_depvar_mean_std = 1;
1139
1140   for (k = 0; k < cmd.n_dependent; k++)
1141     {
1142       n_indep = get_n_indep ((const struct variable *) cmd.v_dependent[k]);
1143       lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
1144       indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
1145       assert (indep_vars != NULL);
1146
1147       for (i = 0; i < n_cases; i++)
1148         {
1149           is_missing_case[i] = 0;
1150         }
1151       n_data = prepare_data (n_cases, is_missing_case, indep_vars,
1152                              cmd.v_dependent[k],
1153                              (const struct casefile *) cf);
1154       Y = gsl_vector_alloc (n_data);
1155
1156       X =
1157         design_matrix_create (n_indep, (const struct variable **) indep_vars,
1158                               n_data);
1159       for (i = 0; i < X->m->size2; i++)
1160         {
1161           lopts.get_indep_mean_std[i] = 1;
1162         }
1163       models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
1164       models[k]->indep_means = gsl_vector_alloc (X->m->size2);
1165       models[k]->indep_std = gsl_vector_alloc (X->m->size2);
1166       models[k]->depvar = (const struct variable *) cmd.v_dependent[k];
1167       /*
1168          For large data sets, use QR decomposition.
1169        */
1170       if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
1171         {
1172           models[k]->method = PSPP_LINREG_SVD;
1173         }
1174
1175       /*
1176          The second pass fills the design matrix.
1177        */
1178       row = 0;
1179       for (r = casefile_get_reader (cf); casereader_read (r, &c);
1180            case_destroy (&c))
1181         /* Iterate over the cases. */
1182         {
1183           case_num = casereader_cnum (r) - 1;
1184           if (!is_missing_case[case_num])
1185             {
1186               for (i = 0; i < n_variables; ++i) /* Iterate over the
1187                                                    variables for the
1188                                                    current case.
1189                                                  */
1190                 {
1191                   val = case_data (&c, v_variables[i]->fv);
1192                   /*
1193                      Independent/dependent variable separation. The
1194                      'variables' subcommand specifies a varlist which contains
1195                      both dependent and independent variables. The dependent
1196                      variables are specified with the 'dependent'
1197                      subcommand, and maybe also in the 'variables' subcommand.
1198                      We need to separate the two.
1199                    */
1200                   if (!is_depvar (i, cmd.v_dependent[k]))
1201                     {
1202                       if (v_variables[i]->type == ALPHA)
1203                         {
1204                           design_matrix_set_categorical (X, row,
1205                                                          v_variables[i], val);
1206                         }
1207                       else if (v_variables[i]->type == NUMERIC)
1208                         {
1209                           design_matrix_set_numeric (X, row, v_variables[i],
1210                                                      val);
1211                         }
1212                     }
1213                 }
1214               val = case_data (&c, cmd.v_dependent[k]->fv);
1215               gsl_vector_set (Y, row, val->f);
1216               row++;
1217             }
1218         }
1219       /*
1220          Now that we know the number of coefficients, allocate space
1221          and store pointers to the variables that correspond to the
1222          coefficients.
1223        */
1224       coeff_init (models[k], X);
1225
1226       /*
1227          Find the least-squares estimates and other statistics.
1228        */
1229       pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
1230       subcommand_statistics (cmd.a_statistics, models[k]);
1231       subcommand_export (cmd.sbc_export, models[k]);
1232
1233       gsl_vector_free (Y);
1234       design_matrix_destroy (X);
1235       free (indep_vars);
1236       free (lopts.get_indep_mean_std);
1237       casereader_destroy (r);
1238     }
1239
1240   free (is_missing_case);
1241
1242   return true;
1243 }
1244
1245 /*
1246   Local Variables:
1247   mode: c
1248   End:
1249 */