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_vals;
 604   struct reg_trns *trns = t_;
 605   pspp_linreg_cache *model;
 606   union value *output = NULL;
 607   const union value **vals = NULL;
 608   const union value *obs = NULL;
 609   struct variable **vars = NULL;
 610
 611   assert (trns!= NULL);
 612   model = trns->c;
 613   assert (model != NULL);
 614   assert (model->depvar != NULL);
 615   assert (model->resid != NULL);
 616
 617   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
 618   n_vals = (*model->get_vars) (model, vars);
 619
 620   vals = xnmalloc (n_vals, sizeof (*vals));
 621   output = case_data_rw (c, model->resid->fv);
 622   assert (output != NULL);
 623
 624   for (i = 0; i < n_vals; i++)
 625     {
 626       vals[i] = case_data (c, vars[i]->fv);
 627     }
 628   obs = case_data (c, model->depvar->fv);
 629   output->f = (*model->residual) ((const struct variable **) vars,
 630                                   vals, obs, model, n_vals);
 631   free (vals);
 632   free (vars);
 633   return TRNS_CONTINUE;
 634 }
 635 /*
 636    Returns 0 if NAME is a duplicate of any existing variable name.
 637 */
 638 static int
 639 try_name (char *name)
 640 {
 641   if (dict_lookup_var (default_dict, name) != NULL)
 642     return 0;
 643
 644   return 1;
 645 }
 646 static
 647 void reg_get_name (char name[LONG_NAME_LEN], const char prefix[LONG_NAME_LEN])
 648 {
 649   int i = 1;
 650
 651   snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
 652   while (!try_name(name))
 653     {
 654       i++;
 655       snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
 656     }
 657 }
 658 static void
 659 reg_save_var (const char *prefix, trns_proc_func *f,
 660               pspp_linreg_cache *c, struct variable **v,
 661               int n_trns)
 662 {
 663   static int trns_index = 1;
 664   char name[LONG_NAME_LEN];
 665   struct variable *new_var;
 666   struct reg_trns *t = NULL;
 667
 668   t = xmalloc (sizeof (*t));
 669   t->trns_id = trns_index;
 670   t->n_trns = n_trns;
 671   t->c = c;
 672   reg_get_name (name, prefix);
 673   new_var = dict_create_var (default_dict, name, 0);
 674   assert (new_var != NULL);
 675   *v = new_var;
 676   add_transformation (f, regression_trns_free, t);
 677   trns_index++;
 678 }
 679 static void
 680 subcommand_save (int save, pspp_linreg_cache **models)
 681 {
 682   pspp_linreg_cache **lc;
 683   int n_trns = 0;
 684   int i;
 685
 686   assert (models != NULL);
 687
 688   if (save)
 689     {
 690       /* Count the number of transformations we will need. */
 691       for (i = 0; i < REGRESSION_SV_count; i++)
 692         {
 693           if (cmd.a_save[i])
 694             {
 695               n_trns++;
 696             }
 697         }
 698       n_trns *= cmd.n_dependent;
 699
 700       for (lc = models; lc < models + cmd.n_dependent; lc++)
 701         {
 702           assert (*lc != NULL);
 703           assert ((*lc)->depvar != NULL);
 704           if (cmd.a_save[REGRESSION_SV_RESID])
 705             {
 706               reg_save_var ("RES", regression_trns_resid_proc, *lc, &(*lc)->resid, n_trns);
 707             }
 708           if (cmd.a_save[REGRESSION_SV_PRED])
 709             {
 710               reg_save_var ("PRED", regression_trns_pred_proc, *lc, &(*lc)->pred, n_trns);
 711             }
 712         }
 713     }
 714   else
 715     {
 716       for (lc = models; lc < models + cmd.n_dependent; lc++)
 717         {
 718           assert (*lc != NULL);
 719           pspp_linreg_cache_free (*lc);
 720         }
 721     }
 722 }
 723 static int
 724 reg_inserted (const struct variable *v, struct variable **varlist, int n_vars)
 725 {
 726   int i;
 727
 728   for (i = 0; i < n_vars; i++)
 729     {
 730       if (v->index == varlist[i]->index)
 731         {
 732           return 1;
 733         }
 734     }
 735   return 0;
 736 }
 737 static void
 738 reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
 739 {
 740   int i;
 741   size_t j;
 742   int n_vars = 0;
 743   struct variable **varlist;
 744   struct pspp_linreg_coeff *coeff;
 745   const struct variable *v;
 746   union value *val;
 747
 748   fprintf (fp, "%s", reg_export_categorical_encode_1);
 749
 750   varlist = xnmalloc (c->n_indeps, sizeof (*varlist));
 751   for (i = 1; i < c->n_indeps; i++)     /* c->coeff[0] is the intercept. */
 752     {
 753       coeff = c->coeff + i;
 754       v = pspp_linreg_coeff_get_var (coeff, 0);
 755       if (v->type == ALPHA)
 756         {
 757           if (!reg_inserted (v, varlist, n_vars))
 758             {
 759               fprintf (fp, "struct pspp_reg_categorical_variable %s;\n\t",
 760                        v->name);
 761               varlist[n_vars] = (struct variable *) v;
 762               n_vars++;
 763             }
 764         }
 765     }
 766   fprintf (fp, "int n_vars = %d;\n\t", n_vars);
 767   fprintf (fp, "struct pspp_reg_categorical_variable *varlist[%d] = {",
 768            n_vars);
 769   for (i = 0; i < n_vars - 1; i++)
 770     {
 771       fprintf (fp, "&%s,\n\t\t", varlist[i]->name);
 772     }
 773   fprintf (fp, "&%s};\n\t", varlist[i]->name);
 774
 775   for (i = 0; i < n_vars; i++)
 776     {
 777       coeff = c->coeff + i;
 778       fprintf (fp, "%s.name = \"%s\";\n\t", varlist[i]->name,
 779                varlist[i]->name);
 780       fprintf (fp, "%s.n_vals = %d;\n\t", varlist[i]->name,
 781                varlist[i]->obs_vals->n_categories);
 782
 783       for (j = 0; j < varlist[i]->obs_vals->n_categories; j++)
 784         {
 785           val = cat_subscript_to_value ((const size_t) j, varlist[i]);
 786           fprintf (fp, "%s.values[%d] = \"%s\";\n\t", varlist[i]->name, j,
 787                    value_to_string (val, varlist[i]));
 788         }
 789     }
 790   fprintf (fp, "%s", reg_export_categorical_encode_2);
 791 }
 792
 793 static void
 794 reg_print_depvars (FILE * fp, pspp_linreg_cache * c)
 795 {
 796   int i;
 797   struct pspp_linreg_coeff *coeff;
 798   const struct variable *v;
 799
 800   fprintf (fp, "char *model_depvars[%d] = {", c->n_indeps);
 801   for (i = 1; i < c->n_indeps; i++)
 802     {
 803       coeff = c->coeff + i;
 804       v = pspp_linreg_coeff_get_var (coeff, 0);
 805       fprintf (fp, "\"%s\",\n\t\t", v->name);
 806     }
 807   coeff = c->coeff + i;
 808   v = pspp_linreg_coeff_get_var (coeff, 0);
 809   fprintf (fp, "\"%s\"};\n\t", v->name);
 810 }
 811 static void
 812 reg_print_getvar (FILE * fp, pspp_linreg_cache * c)
 813 {
 814   fprintf (fp, "static int\npspp_reg_getvar (char *v_name)\n{\n\t");
 815   fprintf (fp, "int i;\n\tint n_vars = %d;\n\t", c->n_indeps);
 816   reg_print_depvars (fp, c);
 817   fprintf (fp, "for (i = 0; i < n_vars; i++)\n\t{\n\t\t");
 818   fprintf (fp,
 819            "if (strncmp (v_name, model_depvars[i], PSPP_REG_MAXLEN) == 0)\n\t\t{\n\t\t\t");
 820   fprintf (fp, "return i;\n\t\t}\n\t}\n}\n");
 821 }
 822 static int
 823 reg_has_categorical (pspp_linreg_cache * c)
 824 {
 825   int i;
 826   const struct variable *v;
 827
 828   for (i = 1; i < c->n_coeffs; i++)
 829     {
 830       v = pspp_linreg_coeff_get_var (c->coeff + i, 0);
 831       if (v->type == ALPHA)
 832         {
 833           return 1;
 834         }
 835     }
 836   return 0;
 837 }
 838
 839 static void
 840 subcommand_export (int export, pspp_linreg_cache * c)
 841 {
 842   FILE *fp;
 843   size_t i;
 844   size_t j;
 845   int n_quantiles = 100;
 846   double increment;
 847   double tmp;
 848   struct pspp_linreg_coeff coeff;
 849
 850   if (export)
 851     {
 852       assert (c != NULL);
 853       assert (model_file != NULL);
 854       fp = fopen (fh_get_file_name (model_file), "w");
 855       assert (fp != NULL);
 856       fprintf (fp, "%s", reg_preamble);
 857       reg_print_getvar (fp, c);
 858       if (reg_has_categorical (c))
 859         {
 860           reg_print_categorical_encoding (fp, c);
 861         }
 862       fprintf (fp, "%s", reg_export_t_quantiles_1);
 863       increment = 0.5 / (double) increment;
 864       for (i = 0; i < n_quantiles - 1; i++)
 865         {
 866           tmp = 0.5 + 0.005 * (double) i;
 867           fprintf (fp, "%.15e,\n\t\t",
 868                    gsl_cdf_tdist_Pinv (tmp, c->n_obs - c->n_indeps));
 869         }
 870       fprintf (fp, "%.15e};\n\t",
 871                gsl_cdf_tdist_Pinv (.9995, c->n_obs - c->n_indeps));
 872       fprintf (fp, "%s", reg_export_t_quantiles_2);
 873       fprintf (fp, "%s", reg_mean_cmt);
 874       fprintf (fp, "double\npspp_reg_estimate (const double *var_vals,");
 875       fprintf (fp, "const char *var_names[])\n{\n\t");
 876       fprintf (fp, "double model_coeffs[%d] = {", c->n_indeps);
 877       for (i = 1; i < c->n_indeps; i++)
 878         {
 879           coeff = c->coeff[i];
 880           fprintf (fp, "%.15e,\n\t\t", coeff.estimate);
 881         }
 882       coeff = c->coeff[i];
 883       fprintf (fp, "%.15e};\n\t", coeff.estimate);
 884       coeff = c->coeff[0];
 885       fprintf (fp, "double estimate = %.15e;\n\t", coeff.estimate);
 886       fprintf (fp, "int i;\n\tint j;\n\n\t");
 887       fprintf (fp, "for (i = 0; i < %d; i++)\n\t", c->n_indeps);
 888       fprintf (fp, "%s", reg_getvar);
 889       fprintf (fp, "const double cov[%d][%d] = {\n\t", c->n_coeffs,
 890                c->n_coeffs);
 891       for (i = 0; i < c->cov->size1 - 1; i++)
 892         {
 893           fprintf (fp, "{");
 894           for (j = 0; j < c->cov->size2 - 1; j++)
 895             {
 896               fprintf (fp, "%.15e, ", gsl_matrix_get (c->cov, i, j));
 897             }
 898           fprintf (fp, "%.15e},\n\t", gsl_matrix_get (c->cov, i, j));
 899         }
 900       fprintf (fp, "{");
 901       for (j = 0; j < c->cov->size2 - 1; j++)
 902         {
 903           fprintf (fp, "%.15e, ",
 904                    gsl_matrix_get (c->cov, c->cov->size1 - 1, j));
 905         }
 906       fprintf (fp, "%.15e}\n\t",
 907                gsl_matrix_get (c->cov, c->cov->size1 - 1, c->cov->size2 - 1));
 908       fprintf (fp, "};\n\tint n_vars = %d;\n\tint i;\n\tint j;\n\t",
 909                c->n_indeps);
 910       fprintf (fp, "double unshuffled_vals[%d];\n\t", c->n_indeps);
 911       fprintf (fp, "%s", reg_variance);
 912       fprintf (fp, "%s", reg_export_confidence_interval);
 913       tmp = c->mse * c->mse;
 914       fprintf (fp, "%s %.15e", reg_export_prediction_interval_1, tmp);
 915       fprintf (fp, "%s %.15e", reg_export_prediction_interval_2, tmp);
 916       fprintf (fp, "%s", reg_export_prediction_interval_3);
 917       fclose (fp);
 918       fp = fopen ("pspp_model_reg.h", "w");
 919       fprintf (fp, "%s", reg_header);
 920       fclose (fp);
 921     }
 922 }
 923 static int
 924 regression_custom_export (struct cmd_regression *cmd UNUSED)
 925 {
 926   /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
 927   if (!lex_force_match ('('))
 928     return 0;
 929
 930   if (lex_match ('*'))
 931     model_file = NULL;
 932   else
 933     {
 934       model_file = fh_parse (FH_REF_FILE);
 935       if (model_file == NULL)
 936         return 0;
 937     }
 938
 939   if (!lex_force_match (')'))
 940     return 0;
 941
 942   return 1;
 943 }
 944
 945 int
 946 cmd_regression (void)
 947 {
 948   if (!parse_regression (&cmd))
 949     return CMD_FAILURE;
 950
 951   models = xnmalloc (cmd.n_dependent, sizeof *models);
 952   if (!multipass_procedure_with_splits (run_regression, &cmd))
 953     return CMD_CASCADING_FAILURE;
 954   subcommand_save (cmd.sbc_save, models);
 955   free (v_variables);
 956   free (models);
 957   return pspp_reg_rc;
 958 }
 959
 960 /*
 961   Is variable k the dependent variable?
 962  */
 963 static int
 964 is_depvar (size_t k, const struct variable *v)
 965 {
 966   /*
 967     compare_var_names returns 0 if the variable
 968     names match.
 969   */
 970   if (!compare_var_names (v, v_variables[k], NULL))
 971     return 1;
 972
 973   return 0;
 974 }
 975
 976 /*
 977   Mark missing cases. Return the number of non-missing cases.
 978  */
 979 static size_t
 980 mark_missing_cases (const struct casefile *cf, struct variable *v,
 981                     int *is_missing_case, double n_data)
 982 {
 983   struct casereader *r;
 984   struct ccase c;
 985   size_t row;
 986   const union value *val;
 987
 988   for (r = casefile_get_reader (cf);
 989        casereader_read (r, &c); case_destroy (&c))
 990     {
 991       row = casereader_cnum (r) - 1;
 992
 993       val = case_data (&c, v->fv);
 994       cat_value_update (v, val);
 995       if (mv_is_value_missing (&v->miss, val))
 996         {
 997           if (!is_missing_case[row])
 998             {
 999               /* Now it is missing. */
1000               n_data--;
1001               is_missing_case[row] = 1;
1002             }
1003         }
1004     }
1005   casereader_destroy (r);
1006
1007   return n_data;
1008 }
1009
1010 /* Parser for the variables sub command */
1011 static int
1012 regression_custom_variables(struct cmd_regression *cmd UNUSED)
1013 {
1014
1015   lex_match('=');
1016
1017   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
1018       && token != T_ALL)
1019     return 2;
1020
1021
1022   if (!parse_variables (default_dict, &v_variables, &n_variables,
1023                         PV_NONE ))
1024     {
1025       free (v_variables);
1026       return 0;
1027     }
1028   assert(n_variables);
1029
1030   return 1;
1031 }
1032 /*
1033   Count the explanatory variables. The user may or may
1034   not have specified a response variable in the syntax.
1035  */
1036 static
1037 int get_n_indep (const struct variable *v)
1038 {
1039   int result;
1040   int i = 0;
1041
1042   result = n_variables;
1043   while (i < n_variables)
1044     {
1045       if (is_depvar (i, v))
1046         {
1047           result--;
1048           i = n_variables;
1049         }
1050       i++;
1051     }
1052   return result;
1053 }
1054 /*
1055   Read from the active file. Identify the explanatory variables in
1056   v_variables. Encode categorical variables. Drop cases with missing
1057   values.
1058 */
1059 static
1060 int prepare_data (int n_data, int is_missing_case[],
1061                   struct variable **indep_vars,
1062                   struct variable *depvar,
1063                   const struct casefile *cf)
1064 {
1065   int i;
1066   int j;
1067
1068   assert (indep_vars != NULL);
1069   j = 0;
1070   for (i = 0; i < n_variables; i++)
1071     {
1072       if (!is_depvar (i, depvar))
1073         {
1074           indep_vars[j] = v_variables[i];
1075           j++;
1076           if (v_variables[i]->type == ALPHA)
1077             {
1078               /* Make a place to hold the binary vectors
1079                  corresponding to this variable's values. */
1080               cat_stored_values_create (v_variables[i]);
1081             }
1082           n_data = mark_missing_cases (cf, v_variables[i], is_missing_case, n_data);
1083         }
1084     }
1085   /*
1086     Mark missing cases for the dependent variable.
1087    */
1088   n_data = mark_missing_cases (cf, depvar, is_missing_case, n_data);
1089
1090   return n_data;
1091 }
1092 static bool
1093 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
1094 {
1095   size_t i;
1096   size_t n_data = 0; /* Number of valide cases. */
1097   size_t n_cases; /* Number of cases. */
1098   size_t row;
1099   size_t case_num;
1100   int n_indep = 0;
1101   int k;
1102   /*
1103      Keep track of the missing cases.
1104    */
1105   int *is_missing_case;
1106   const union value *val;
1107   struct casereader *r;
1108   struct ccase c;
1109   struct variable **indep_vars;
1110   struct design_matrix *X;
1111   gsl_vector *Y;
1112
1113   pspp_linreg_opts lopts;
1114
1115   assert (models != NULL);
1116   if (!v_variables)
1117     {
1118       dict_get_vars (default_dict, &v_variables, &n_variables,
1119                      1u << DC_SYSTEM);
1120     }
1121
1122   n_cases = casefile_get_case_cnt (cf);
1123
1124   for (i = 0; i < cmd.n_dependent; i++)
1125     {
1126       if (cmd.v_dependent[i]->type != NUMERIC)
1127         {
1128           msg (SE, gettext ("Dependent variable must be numeric."));
1129           pspp_reg_rc = CMD_FAILURE;
1130           return true;
1131         }
1132     }
1133
1134   is_missing_case = xnmalloc (n_cases, sizeof (*is_missing_case));
1135
1136   lopts.get_depvar_mean_std = 1;
1137
1138   for (k = 0; k < cmd.n_dependent; k++)
1139     {
1140       n_indep = get_n_indep ((const struct variable *) cmd.v_dependent[k]);
1141       lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
1142       indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
1143       assert (indep_vars != NULL);
1144
1145       for (i = 0; i < n_cases; i++)
1146         {
1147           is_missing_case[i] = 0;
1148         }
1149       n_data = prepare_data (n_cases, is_missing_case, indep_vars,
1150                              cmd.v_dependent[k],
1151                              (const struct casefile *) cf);
1152       Y = gsl_vector_alloc (n_data);
1153
1154       X =
1155         design_matrix_create (n_indep, (const struct variable **) indep_vars,
1156                               n_data);
1157       for (i = 0; i < X->m->size2; i++)
1158         {
1159           lopts.get_indep_mean_std[i] = 1;
1160         }
1161       models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
1162       models[k]->indep_means = gsl_vector_alloc (X->m->size2);
1163       models[k]->indep_std = gsl_vector_alloc (X->m->size2);
1164       models[k]->depvar = (const struct variable *) cmd.v_dependent[k];
1165       /*
1166          For large data sets, use QR decomposition.
1167        */
1168       if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
1169         {
1170           models[k]->method = PSPP_LINREG_SVD;
1171         }
1172
1173       /*
1174          The second pass fills the design matrix.
1175        */
1176       row = 0;
1177       for (r = casefile_get_reader (cf); casereader_read (r, &c);
1178            case_destroy (&c))
1179         /* Iterate over the cases. */
1180         {
1181           case_num = casereader_cnum (r) - 1;
1182           if (!is_missing_case[case_num])
1183             {
1184               for (i = 0; i < n_variables; ++i) /* Iterate over the
1185                                                    variables for the
1186                                                    current case.
1187                                                 */
1188                 {
1189                   val = case_data (&c, v_variables[i]->fv);
1190                   /*
1191                      Independent/dependent variable separation. The
1192                      'variables' subcommand specifies a varlist which contains
1193                      both dependent and independent variables. The dependent
1194                      variables are specified with the 'dependent'
1195                      subcommand, and maybe also in the 'variables' subcommand.
1196                      We need to separate the two.
1197                    */
1198                   if (!is_depvar (i, cmd.v_dependent[k]))
1199                     {
1200                       if (v_variables[i]->type == ALPHA)
1201                         {
1202                           design_matrix_set_categorical (X, row, v_variables[i], val);
1203                         }
1204                       else if (v_variables[i]->type == NUMERIC)
1205                         {
1206                           design_matrix_set_numeric (X, row, v_variables[i], val);
1207                         }
1208                     }
1209                 }
1210               val = case_data (&c, cmd.v_dependent[k]->fv);
1211               gsl_vector_set (Y, row, val->f);
1212               row++;
1213             }
1214         }
1215       /*
1216          Now that we know the number of coefficients, allocate space
1217          and store pointers to the variables that correspond to the
1218          coefficients.
1219        */
1220       pspp_linreg_coeff_init (models[k], X);
1221
1222       /*
1223          Find the least-squares estimates and other statistics.
1224        */
1225       pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
1226       subcommand_statistics (cmd.a_statistics, models[k]);
1227       subcommand_export (cmd.sbc_export, models[k]);
1228
1229       gsl_vector_free (Y);
1230       design_matrix_destroy (X);
1231       free (indep_vars);
1232       free (lopts.get_indep_mean_std);
1233       casereader_destroy (r);
1234     }
1235
1236   free (is_missing_case);
1237
1238   return true;
1239 }
1240
1241 /*
1242   Local Variables:
1243   mode: c
1244   End:
1245 */