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