src/language/stats/regression.q

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