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