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