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