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