pintos-os.org Git - pspp/blob - src/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 "alloc.h"
  26 #include "case.h"
  27 #include "casefile.h"
  28 #include "cat.h"
  29 #include "cat-routines.h"
  30 #include "command.h"
  31 #include "design-matrix.h"
  32 #include "dictionary.h"
  33 #include "error.h"
  34 #include "file-handle.h"
  35 #include "gettext.h"
  36 #include "lexer.h"
  37 #include <linreg/pspp_linreg.h>
  38 #include "missing-values.h"
  39 #include "tab.h"
  40 #include "var.h"
  41 #include "vfm.h"
  42
  43 #define REG_LARGE_DATA 1000
  44
  45 /* (headers) */
  46
  47 /* (specification)
  48    "REGRESSION" (regression_):
  49    *variables=varlist;
  50    statistics[st_]=r,
  51    coeff,
  52    anova,
  53    outs,
  54    zpp,
  55    label,
  56    sha,
  57    ci,
  58    bcov,
  59    ses,
  60    xtx,
  61    collin,
  62    tol,
  63    selection,
  64    f,
  65    defaults,
  66    all;
  67    export=custom;
  68    ^dependent=varlist;
  69    ^method=enter.
  70 */
  71 /* (declarations) */
  72 /* (functions) */
  73 static struct cmd_regression cmd;
  74
  75 /*
  76   Array holding the subscripts of the independent variables.
  77  */
  78 size_t *indep_vars;
  79
  80 /*
  81   File where the model will be saved if the EXPORT subcommand
  82   is given.
  83  */
  84 struct file_handle *model_file;
  85
  86 /*
  87   Return value for the procedure.
  88  */
  89 int pspp_reg_rc = CMD_SUCCESS;
  90
  91 static void run_regression (const struct casefile *, void *);
  92 /*
  93    STATISTICS subcommand output functions.
  94  */
  95 static void reg_stats_r (pspp_linreg_cache *);
  96 static void reg_stats_coeff (pspp_linreg_cache *);
  97 static void reg_stats_anova (pspp_linreg_cache *);
  98 static void reg_stats_outs (pspp_linreg_cache *);
  99 static void reg_stats_zpp (pspp_linreg_cache *);
 100 static void reg_stats_label (pspp_linreg_cache *);
 101 static void reg_stats_sha (pspp_linreg_cache *);
 102 static void reg_stats_ci (pspp_linreg_cache *);
 103 static void reg_stats_f (pspp_linreg_cache *);
 104 static void reg_stats_bcov (pspp_linreg_cache *);
 105 static void reg_stats_ses (pspp_linreg_cache *);
 106 static void reg_stats_xtx (pspp_linreg_cache *);
 107 static void reg_stats_collin (pspp_linreg_cache *);
 108 static void reg_stats_tol (pspp_linreg_cache *);
 109 static void reg_stats_selection (pspp_linreg_cache *);
 110 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
 111                                        int, pspp_linreg_cache *);
 112
 113 static void
 114 reg_stats_r (pspp_linreg_cache * c)
 115 {
 116   struct tab_table *t;
 117   int n_rows = 2;
 118   int n_cols = 5;
 119   double rsq;
 120   double adjrsq;
 121   double std_error;
 122
 123   assert (c != NULL);
 124   rsq = c->ssm / c->sst;
 125   adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
 126   std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
 127   t = tab_create (n_cols, n_rows, 0);
 128   tab_dim (t, tab_natural_dimensions);
 129   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 130   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 131   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 132   tab_vline (t, TAL_0, 1, 0, 0);
 133
 134   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
 135   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
 136   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
 137   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
 138   tab_float (t, 1, 1, TAB_RIGHT, sqrt (rsq), 10, 2);
 139   tab_float (t, 2, 1, TAB_RIGHT, rsq, 10, 2);
 140   tab_float (t, 3, 1, TAB_RIGHT, adjrsq, 10, 2);
 141   tab_float (t, 4, 1, TAB_RIGHT, std_error, 10, 2);
 142   tab_title (t, 0, _("Model Summary"));
 143   tab_submit (t);
 144 }
 145
 146 /*
 147   Table showing estimated regression coefficients.
 148  */
 149 static void
 150 reg_stats_coeff (pspp_linreg_cache * c)
 151 {
 152   size_t i;
 153   size_t j;
 154   int n_cols = 7;
 155   int n_rows;
 156   double t_stat;
 157   double pval;
 158   double coeff;
 159   double std_err;
 160   double beta;
 161   const char *label;
 162   struct tab_table *t;
 163
 164   assert (c != NULL);
 165   n_rows = c->n_coeffs + 2;
 166
 167   t = tab_create (n_cols, n_rows, 0);
 168   tab_headers (t, 2, 0, 1, 0);
 169   tab_dim (t, tab_natural_dimensions);
 170   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 171   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 172   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 173   tab_vline (t, TAL_0, 1, 0, 0);
 174
 175   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
 176   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 177   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
 178   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
 179   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 180   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
 181   coeff = c->coeff[0].estimate;
 182   tab_float (t, 2, 1, 0, coeff, 10, 2);
 183   std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
 184   tab_float (t, 3, 1, 0, std_err, 10, 2);
 185   beta = coeff / c->depvar_std;
 186   tab_float (t, 4, 1, 0, beta, 10, 2);
 187   t_stat = coeff / std_err;
 188   tab_float (t, 5, 1, 0, t_stat, 10, 2);
 189   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 190   tab_float (t, 6, 1, 0, pval, 10, 2);
 191   for (j = 1; j <= c->n_indeps; j++)
 192     {
 193       i = indep_vars[j];
 194       label = var_to_string (c->coeff[j].v);
 195       tab_text (t, 1, j + 1, TAB_CENTER, label);
 196       /*
 197          Regression coefficients.
 198        */
 199       coeff = c->coeff[j].estimate;
 200       tab_float (t, 2, j + 1, 0, coeff, 10, 2);
 201       /*
 202          Standard error of the coefficients.
 203        */
 204       std_err = sqrt (gsl_matrix_get (c->cov, j, j));
 205       tab_float (t, 3, j + 1, 0, std_err, 10, 2);
 206       /*
 207          'Standardized' coefficient, i.e., regression coefficient
 208          if all variables had unit variance.
 209        */
 210       beta = gsl_vector_get (c->indep_std, j);
 211       beta *= coeff / c->depvar_std;
 212       tab_float (t, 4, j + 1, 0, beta, 10, 2);
 213
 214       /*
 215          Test statistic for H0: coefficient is 0.
 216        */
 217       t_stat = coeff / std_err;
 218       tab_float (t, 5, j + 1, 0, t_stat, 10, 2);
 219       /*
 220          P values for the test statistic above.
 221        */
 222       pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 223       tab_float (t, 6, j + 1, 0, pval, 10, 2);
 224     }
 225   tab_title (t, 0, _("Coefficients"));
 226   tab_submit (t);
 227 }
 228
 229 /*
 230   Display the ANOVA table.
 231  */
 232 static void
 233 reg_stats_anova (pspp_linreg_cache * c)
 234 {
 235   int n_cols = 7;
 236   int n_rows = 4;
 237   const double msm = c->ssm / c->dfm;
 238   const double mse = c->sse / c->dfe;
 239   const double F = msm / mse;
 240   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
 241
 242   struct tab_table *t;
 243
 244   assert (c != NULL);
 245   t = tab_create (n_cols, n_rows, 0);
 246   tab_headers (t, 2, 0, 1, 0);
 247   tab_dim (t, tab_natural_dimensions);
 248
 249   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 250
 251   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 252   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 253   tab_vline (t, TAL_0, 1, 0, 0);
 254
 255   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 256   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 257   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 258   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 259   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 260
 261   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
 262   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
 263   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
 264
 265   /* Sums of Squares */
 266   tab_float (t, 2, 1, 0, c->ssm, 10, 2);
 267   tab_float (t, 2, 3, 0, c->sst, 10, 2);
 268   tab_float (t, 2, 2, 0, c->sse, 10, 2);
 269
 270
 271   /* Degrees of freedom */
 272   tab_float (t, 3, 1, 0, c->dfm, 4, 0);
 273   tab_float (t, 3, 2, 0, c->dfe, 4, 0);
 274   tab_float (t, 3, 3, 0, c->dft, 4, 0);
 275
 276   /* Mean Squares */
 277
 278   tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
 279   tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
 280
 281   tab_float (t, 5, 1, 0, F, 8, 3);
 282
 283   tab_float (t, 6, 1, 0, pval, 8, 3);
 284
 285   tab_title (t, 0, _("ANOVA"));
 286   tab_submit (t);
 287 }
 288 static void
 289 reg_stats_outs (pspp_linreg_cache * c)
 290 {
 291   assert (c != NULL);
 292 }
 293 static void
 294 reg_stats_zpp (pspp_linreg_cache * c)
 295 {
 296   assert (c != NULL);
 297 }
 298 static void
 299 reg_stats_label (pspp_linreg_cache * c)
 300 {
 301   assert (c != NULL);
 302 }
 303 static void
 304 reg_stats_sha (pspp_linreg_cache * c)
 305 {
 306   assert (c != NULL);
 307 }
 308 static void
 309 reg_stats_ci (pspp_linreg_cache * c)
 310 {
 311   assert (c != NULL);
 312 }
 313 static void
 314 reg_stats_f (pspp_linreg_cache * c)
 315 {
 316   assert (c != NULL);
 317 }
 318 static void
 319 reg_stats_bcov (pspp_linreg_cache * c)
 320 {
 321   int n_cols;
 322   int n_rows;
 323   int i;
 324   int j;
 325   int k;
 326   int row;
 327   int col;
 328   const char *label;
 329   struct tab_table *t;
 330
 331   assert (c != NULL);
 332   n_cols = c->n_indeps + 1 + 2;
 333   n_rows = 2 * (c->n_indeps + 1);
 334   t = tab_create (n_cols, n_rows, 0);
 335   tab_headers (t, 2, 0, 1, 0);
 336   tab_dim (t, tab_natural_dimensions);
 337   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 338   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 339   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 340   tab_vline (t, TAL_0, 1, 0, 0);
 341   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 342   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 343   for (i = 1; i < c->n_indeps + 1; i++)
 344     {
 345       j = indep_vars[(i - 1)];
 346       struct variable *v = cmd.v_variables[j];
 347       label = var_to_string (v);
 348       tab_text (t, 2, i, TAB_CENTER, label);
 349       tab_text (t, i + 2, 0, TAB_CENTER, label);
 350       for (k = 1; k < c->n_indeps + 1; k++)
 351         {
 352           col = (i <= k) ? k : i;
 353           row = (i <= k) ? i : k;
 354           tab_float (t, k + 2, i, TAB_CENTER,
 355                      gsl_matrix_get (c->cov, row, col), 8, 3);
 356         }
 357     }
 358   tab_title (t, 0, _("Coefficient Correlations"));
 359   tab_submit (t);
 360 }
 361 static void
 362 reg_stats_ses (pspp_linreg_cache * c)
 363 {
 364   assert (c != NULL);
 365 }
 366 static void
 367 reg_stats_xtx (pspp_linreg_cache * c)
 368 {
 369   assert (c != NULL);
 370 }
 371 static void
 372 reg_stats_collin (pspp_linreg_cache * c)
 373 {
 374   assert (c != NULL);
 375 }
 376 static void
 377 reg_stats_tol (pspp_linreg_cache * c)
 378 {
 379   assert (c != NULL);
 380 }
 381 static void
 382 reg_stats_selection (pspp_linreg_cache * c)
 383 {
 384   assert (c != NULL);
 385 }
 386
 387 static void
 388 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 389                            int keyword, pspp_linreg_cache * c)
 390 {
 391   if (keyword)
 392     {
 393       (*function) (c);
 394     }
 395 }
 396
 397 static void
 398 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 399 {
 400   /*
 401      The order here must match the order in which the STATISTICS
 402      keywords appear in the specification section above.
 403    */
 404   enum
 405   { r,
 406     coeff,
 407     anova,
 408     outs,
 409     zpp,
 410     label,
 411     sha,
 412     ci,
 413     bcov,
 414     ses,
 415     xtx,
 416     collin,
 417     tol,
 418     selection,
 419     f,
 420     defaults,
 421     all
 422   };
 423   int i;
 424   int d = 1;
 425
 426   if (keywords[all])
 427     {
 428       /*
 429          Set everything but F.
 430        */
 431       for (i = 0; i < f; i++)
 432         {
 433           keywords[i] = 1;
 434         }
 435     }
 436   else
 437     {
 438       for (i = 0; i < all; i++)
 439         {
 440           if (keywords[i])
 441             {
 442               d = 0;
 443             }
 444         }
 445       /*
 446          Default output: ANOVA table, parameter estimates,
 447          and statistics for variables not entered into model,
 448          if appropriate.
 449        */
 450       if (keywords[defaults] | d)
 451         {
 452           keywords[anova] = 1;
 453           keywords[outs] = 1;
 454           keywords[coeff] = 1;
 455           keywords[r] = 1;
 456         }
 457     }
 458   statistics_keyword_output (reg_stats_r, keywords[r], c);
 459   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 460   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 461   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 462   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 463   statistics_keyword_output (reg_stats_label, keywords[label], c);
 464   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 465   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 466   statistics_keyword_output (reg_stats_f, keywords[f], c);
 467   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 468   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 469   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 470   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 471   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 472   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 473 }
 474 static void
 475 subcommand_export (int export, pspp_linreg_cache *c)
 476 {
 477   FILE *fp;
 478   size_t i;
 479   struct pspp_linreg_coeff coeff;
 480
 481   if (export)
 482     {
 483       assert (c != NULL);
 484       assert (model_file != NULL);
 485       assert (fp != NULL);
 486       fp = fopen (handle_get_filename (model_file), "w");
 487       fprintf (fp, "/* PSPP-generated linear regression model.\n   Copyright (C) 2005 Free Software Foundation, Inc.\n   Generated by the GNU PSPP regression procedure.\n\n   This program is free software; you can redistribute it and/or\n   modify it under the terms of the GNU General Public License as\n   published by the Free Software Foundation; either version 2 of the\n   License, or (at your option) any later version.\n\n   This program is distributed in the hope that it will be useful, but\n   WITHOUT ANY WARRANTY; without even the implied warranty of\n   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU\n   General Public License for more details.\n\n   You should have received a copy of the GNU General Public License\n   along with GNU PSPP; if not, write to the Free Software\n   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA\n   02110-1301, USA. */\n\n");
 488       fprintf (fp, "#include <string.h>\n\n");
 489       fprintf (fp, "double\npspp_reg_estimate (const double *var_vals, conts char *[] var_names)\n{\n\tchar *model_depvars[%d] = {", c->n_indeps);
 490       for (i = 1; i < c->n_indeps; i++)
 491         {
 492           coeff = c->coeff[i];
 493           fprintf (fp, "%s,\n\t\t", coeff.v->name);
 494         }
 495       coeff = c->coeff[i];
 496       fprintf (fp, "%s};\n\t", coeff.v->name);
 497       fprintf (fp, "double model_coeffs[%d] = {", c->n_indeps);
 498       for (i = 1; i < c->n_indeps; i++)
 499         {
 500           coeff = c->coeff[i];
 501           fprintf (fp, "%.15e,\n\t\t", coeff.estimate);
 502         }
 503       coeff = c->coeff[i];
 504       fprintf (fp, "%.15e};\n\t", coeff.estimate);
 505       coeff = c->coeff[0];
 506       fprintf (fp, "double estimate = %.15e\n\t", coeff.estimate);
 507       fprintf (fp, "int i;\n\tint j;\n\n\t");
 508       fprintf (fp, "for (i = 0; i < %d; i++)\n\t", c->n_indeps);
 509       fprintf (fp, "{\n\t\tfor (j = 0; j < %d; j++)\n\t\t", c->n_indeps);
 510       fprintf (fp, "{\n\t\t\tif (strcmp (var_names[i], model_names[j]) == 0)\n");
 511       fprintf (fp, "\t\t\t{\n\t\t\t\testimate += var_vals[i] * model_coeffs[j];\n");
 512       fprintf (fp, "\t\t\t}\n\t\t}\n\t}\n\treturn estimate;\n}\n");
 513       fclose (fp);
 514     }
 515 }
 516 static int
 517 regression_custom_export (struct cmd_regression *cmd)
 518 {
 519   /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
 520   if (!lex_force_match ('('))
 521     return 0;
 522
 523   if (lex_match ('*'))
 524     model_file = NULL;
 525   else
 526     {
 527       model_file = fh_parse ();
 528       if (model_file == NULL)
 529         return 0;
 530     }
 531
 532   if (!lex_force_match (')'))
 533     return 0;
 534
 535   return 1;
 536 }
 537
 538 int
 539 cmd_regression (void)
 540 {
 541   if (!parse_regression (&cmd))
 542     {
 543       return CMD_FAILURE;
 544     }
 545   multipass_procedure_with_splits (run_regression, &cmd);
 546
 547   return pspp_reg_rc;
 548 }
 549
 550 /*
 551   Is variable k one of the dependent variables?
 552  */
 553 static int
 554 is_depvar (size_t k)
 555 {
 556   size_t j = 0;
 557   for (j = 0; j < cmd.n_dependent; j++)
 558     {
 559       /*
 560          compare_var_names returns 0 if the variable
 561          names match.
 562        */
 563       if (!compare_var_names (cmd.v_dependent[j], cmd.v_variables[k], NULL))
 564         return 1;
 565     }
 566   return 0;
 567 }
 568
 569 static void
 570 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
 571 {
 572   size_t i;
 573   size_t n_data = 0;
 574   size_t row;
 575   size_t case_num;
 576   int n_indep;
 577   int j = 0;
 578   /*
 579      Keep track of the missing cases.
 580    */
 581   int *is_missing_case;
 582   const union value *val;
 583   struct casereader *r;
 584   struct casereader *r2;
 585   struct ccase c;
 586   struct variable *v;
 587   struct variable **indep_vars;
 588   struct design_matrix *X;
 589   gsl_vector *Y;
 590   pspp_linreg_cache *lcache;
 591   pspp_linreg_opts lopts;
 592
 593   n_data = casefile_get_case_cnt (cf);
 594
 595   is_missing_case = xnmalloc (n_data, sizeof (*is_missing_case));
 596   for (i = 0; i < n_data; i++)
 597     is_missing_case[i] = 0;
 598
 599   n_indep = cmd.n_variables - cmd.n_dependent;
 600   indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
 601
 602   lopts.get_depvar_mean_std = 1;
 603   lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
 604
 605
 606   /*
 607      Read from the active file. The first pass encodes categorical
 608      variables and drops cases with missing values.
 609    */
 610   j = 0;
 611   for (i = 0; i < cmd.n_variables; i++)
 612     {
 613       if (!is_depvar (i))
 614         {
 615           v = cmd.v_variables[i];
 616           indep_vars[j] = v;
 617           j++;
 618           if (v->type == ALPHA)
 619             {
 620               /* Make a place to hold the binary vectors
 621                  corresponding to this variable's values. */
 622               cat_stored_values_create (v);
 623             }
 624           for (r = casefile_get_reader (cf);
 625                casereader_read (r, &c); case_destroy (&c))
 626             {
 627               row = casereader_cnum (r) - 1;
 628
 629               val = case_data (&c, v->fv);
 630               cat_value_update (v, val);
 631               if (mv_is_value_missing (&v->miss, val))
 632                 {
 633                   if (!is_missing_case[row])
 634                     {
 635                       /* Now it is missing. */
 636                       n_data--;
 637                       is_missing_case[row] = 1;
 638                     }
 639                 }
 640             }
 641         }
 642     }
 643
 644   Y = gsl_vector_alloc (n_data);
 645   X =
 646     design_matrix_create (n_indep, (const struct variable **) indep_vars,
 647                           n_data);
 648   lcache = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
 649   lcache->indep_means = gsl_vector_alloc (X->m->size2);
 650   lcache->indep_std = gsl_vector_alloc (X->m->size2);
 651
 652   /*
 653      The second pass creates the design matrix.
 654    */
 655   row = 0;
 656   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
 657        case_destroy (&c))
 658     /* Iterate over the cases. */
 659     {
 660       case_num = casereader_cnum (r2) - 1;
 661       if (!is_missing_case[case_num])
 662         {
 663           for (i = 0; i < cmd.n_variables; ++i) /* Iterate over the variables
 664                                                    for the current case.
 665                                                  */
 666             {
 667               v = cmd.v_variables[i];
 668               val = case_data (&c, v->fv);
 669               /*
 670                  Independent/dependent variable separation. The
 671                  'variables' subcommand specifies a varlist which contains
 672                  both dependent and independent variables. The dependent
 673                  variables are specified with the 'dependent'
 674                  subcommand. We need to separate the two.
 675                */
 676               if (is_depvar (i))
 677                 {
 678                   if (v->type != NUMERIC)
 679                     {
 680                       msg (SE,
 681                            gettext ("Dependent variable must be numeric."));
 682                       pspp_reg_rc = CMD_FAILURE;
 683                       return;
 684                     }
 685                   lcache->depvar = (const struct variable *) v;
 686                   gsl_vector_set (Y, row, val->f);
 687                 }
 688               else
 689                 {
 690                   if (v->type == ALPHA)
 691                     {
 692                       design_matrix_set_categorical (X, row, v, val);
 693                     }
 694                   else if (v->type == NUMERIC)
 695                     {
 696                       design_matrix_set_numeric (X, row, v, val);
 697                     }
 698
 699                   lopts.get_indep_mean_std[i] = 1;
 700                 }
 701             }
 702           row++;
 703         }
 704     }
 705   /*
 706      Now that we know the number of coefficients, allocate space
 707      and store pointers to the variables that correspond to the
 708      coefficients.
 709    */
 710   lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
 711   for (i = 0; i < X->m->size2; i++)
 712     {
 713       j = i + 1;                /* The first coeff is the intercept. */
 714       lcache->coeff[j].v =
 715         (const struct variable *) design_matrix_col_to_var (X, i);
 716       assert (lcache->coeff[j].v != NULL);
 717     }
 718   /*
 719     For large data sets, use QR decomposition.
 720    */
 721   if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
 722     {
 723       lcache->method = PSPP_LINREG_SVD;
 724     }
 725   /*
 726      Find the least-squares estimates and other statistics.
 727    */
 728   pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
 729   subcommand_statistics (cmd.a_statistics, lcache);
 730   subcommand_export (cmd.sbc_export, lcache);
 731   gsl_vector_free (Y);
 732   design_matrix_destroy (X);
 733   pspp_linreg_cache_free (lcache);
 734   free (lopts.get_indep_mean_std);
 735   free (indep_vars);
 736   free (is_missing_case);
 737   casereader_destroy (r);
 738   return;
 739 }
 740
 741 /*
 742   Local Variables:
 743   mode: c
 744   End:
 745 */