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