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