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, "#include <string.h>\n\n");
 488       fprintf (fp, "/*\n   Estimate the mean of Y, the dependent variable for\n");
 489       fprintf (fp, "   the linear model of the form \n\n");
 490       fprintf (fp, "      Y = b0 + b1 * X1 + b2 * X2 + ... + bk * X2 + error\n\n");
 491       fprintf (fp, "   where X1, ..., Xk are the independent variables\n");
 492       fprintf (fp, "   whose values are stored in var_vals and whose names, \n");
 493       fprintf (fp, "   as known by PSPP, are stored in var_names. The estimated \n");
 494       fprintf (fp, "   regression coefficients (i.e., the estimates of b0,...,bk) \n");
 495       fprintf (fp, "   are stored in model_coeffs.\n*/\n");
 496       fprintf (fp, "double\npspp_reg_estimate (const double *var_vals, const char *var_names[])\n{\n\tchar *model_depvars[%d] = {", c->n_indeps);
 497       for (i = 1; i < c->n_indeps; i++)
 498         {
 499           coeff = c->coeff[i];
 500           fprintf (fp, "\"%s\",\n\t\t", coeff.v->name);
 501         }
 502       coeff = c->coeff[i];
 503       fprintf (fp, "\"%s\"};\n\t", coeff.v->name);
 504       fprintf (fp, "double model_coeffs[%d] = {", c->n_indeps);
 505       for (i = 1; i < c->n_indeps; i++)
 506         {
 507           coeff = c->coeff[i];
 508           fprintf (fp, "%.15e,\n\t\t", coeff.estimate);
 509         }
 510       coeff = c->coeff[i];
 511       fprintf (fp, "%.15e};\n\t", coeff.estimate);
 512       coeff = c->coeff[0];
 513       fprintf (fp, "double estimate = %.15e;\n\t", coeff.estimate);
 514       fprintf (fp, "int i;\n\tint j;\n\n\t");
 515       fprintf (fp, "for (i = 0; i < %d; i++)\n\t", c->n_indeps);
 516       fprintf (fp, "{\n\t\tfor (j = 0; j < %d; j++)\n\t\t", c->n_indeps);
 517       fprintf (fp, "{\n\t\t\tif (strcmp (var_names[i], model_depvars[j]) == 0)\n");
 518       fprintf (fp, "\t\t\t{\n\t\t\t\testimate += var_vals[i] * model_coeffs[j];\n");
 519       fprintf (fp, "\t\t\t}\n\t\t}\n\t}\n\treturn estimate;\n}\n");
 520       fclose (fp);
 521     }
 522 }
 523 static int
 524 regression_custom_export (struct cmd_regression *cmd)
 525 {
 526   /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
 527   if (!lex_force_match ('('))
 528     return 0;
 529
 530   if (lex_match ('*'))
 531     model_file = NULL;
 532   else
 533     {
 534       model_file = fh_parse ();
 535       if (model_file == NULL)
 536         return 0;
 537     }
 538
 539   if (!lex_force_match (')'))
 540     return 0;
 541
 542   return 1;
 543 }
 544
 545 int
 546 cmd_regression (void)
 547 {
 548   if (!parse_regression (&cmd))
 549     {
 550       return CMD_FAILURE;
 551     }
 552   multipass_procedure_with_splits (run_regression, &cmd);
 553
 554   return pspp_reg_rc;
 555 }
 556
 557 /*
 558   Is variable k one of the dependent variables?
 559  */
 560 static int
 561 is_depvar (size_t k)
 562 {
 563   size_t j = 0;
 564   for (j = 0; j < cmd.n_dependent; j++)
 565     {
 566       /*
 567          compare_var_names returns 0 if the variable
 568          names match.
 569        */
 570       if (!compare_var_names (cmd.v_dependent[j], cmd.v_variables[k], NULL))
 571         return 1;
 572     }
 573   return 0;
 574 }
 575
 576 static void
 577 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
 578 {
 579   size_t i;
 580   size_t n_data = 0;
 581   size_t row;
 582   size_t case_num;
 583   int n_indep;
 584   int j = 0;
 585   /*
 586      Keep track of the missing cases.
 587    */
 588   int *is_missing_case;
 589   const union value *val;
 590   struct casereader *r;
 591   struct casereader *r2;
 592   struct ccase c;
 593   struct variable *v;
 594   struct variable **indep_vars;
 595   struct design_matrix *X;
 596   gsl_vector *Y;
 597   pspp_linreg_cache *lcache;
 598   pspp_linreg_opts lopts;
 599
 600   n_data = casefile_get_case_cnt (cf);
 601
 602   is_missing_case = xnmalloc (n_data, sizeof (*is_missing_case));
 603   for (i = 0; i < n_data; i++)
 604     is_missing_case[i] = 0;
 605
 606   n_indep = cmd.n_variables - cmd.n_dependent;
 607   indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
 608
 609   lopts.get_depvar_mean_std = 1;
 610   lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
 611
 612
 613   /*
 614      Read from the active file. The first pass encodes categorical
 615      variables and drops cases with missing values.
 616    */
 617   j = 0;
 618   for (i = 0; i < cmd.n_variables; i++)
 619     {
 620       if (!is_depvar (i))
 621         {
 622           v = cmd.v_variables[i];
 623           indep_vars[j] = v;
 624           j++;
 625           if (v->type == ALPHA)
 626             {
 627               /* Make a place to hold the binary vectors
 628                  corresponding to this variable's values. */
 629               cat_stored_values_create (v);
 630             }
 631           for (r = casefile_get_reader (cf);
 632                casereader_read (r, &c); case_destroy (&c))
 633             {
 634               row = casereader_cnum (r) - 1;
 635
 636               val = case_data (&c, v->fv);
 637               cat_value_update (v, val);
 638               if (mv_is_value_missing (&v->miss, val))
 639                 {
 640                   if (!is_missing_case[row])
 641                     {
 642                       /* Now it is missing. */
 643                       n_data--;
 644                       is_missing_case[row] = 1;
 645                     }
 646                 }
 647             }
 648         }
 649     }
 650
 651   Y = gsl_vector_alloc (n_data);
 652   X =
 653     design_matrix_create (n_indep, (const struct variable **) indep_vars,
 654                           n_data);
 655   lcache = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
 656   lcache->indep_means = gsl_vector_alloc (X->m->size2);
 657   lcache->indep_std = gsl_vector_alloc (X->m->size2);
 658
 659   /*
 660      The second pass creates the design matrix.
 661    */
 662   row = 0;
 663   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
 664        case_destroy (&c))
 665     /* Iterate over the cases. */
 666     {
 667       case_num = casereader_cnum (r2) - 1;
 668       if (!is_missing_case[case_num])
 669         {
 670           for (i = 0; i < cmd.n_variables; ++i) /* Iterate over the variables
 671                                                    for the current case.
 672                                                  */
 673             {
 674               v = cmd.v_variables[i];
 675               val = case_data (&c, v->fv);
 676               /*
 677                  Independent/dependent variable separation. The
 678                  'variables' subcommand specifies a varlist which contains
 679                  both dependent and independent variables. The dependent
 680                  variables are specified with the 'dependent'
 681                  subcommand. We need to separate the two.
 682                */
 683               if (is_depvar (i))
 684                 {
 685                   if (v->type != NUMERIC)
 686                     {
 687                       msg (SE,
 688                            gettext ("Dependent variable must be numeric."));
 689                       pspp_reg_rc = CMD_FAILURE;
 690                       return;
 691                     }
 692                   lcache->depvar = (const struct variable *) v;
 693                   gsl_vector_set (Y, row, val->f);
 694                 }
 695               else
 696                 {
 697                   if (v->type == ALPHA)
 698                     {
 699                       design_matrix_set_categorical (X, row, v, val);
 700                     }
 701                   else if (v->type == NUMERIC)
 702                     {
 703                       design_matrix_set_numeric (X, row, v, val);
 704                     }
 705
 706                   lopts.get_indep_mean_std[i] = 1;
 707                 }
 708             }
 709           row++;
 710         }
 711     }
 712   /*
 713      Now that we know the number of coefficients, allocate space
 714      and store pointers to the variables that correspond to the
 715      coefficients.
 716    */
 717   lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
 718   for (i = 0; i < X->m->size2; i++)
 719     {
 720       j = i + 1;                /* The first coeff is the intercept. */
 721       lcache->coeff[j].v =
 722         (const struct variable *) design_matrix_col_to_var (X, i);
 723       assert (lcache->coeff[j].v != NULL);
 724     }
 725   /*
 726     For large data sets, use QR decomposition.
 727    */
 728   if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
 729     {
 730       lcache->method = PSPP_LINREG_SVD;
 731     }
 732   /*
 733      Find the least-squares estimates and other statistics.
 734    */
 735   pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
 736   subcommand_statistics (cmd.a_statistics, lcache);
 737   subcommand_export (cmd.sbc_export, lcache);
 738   gsl_vector_free (Y);
 739   design_matrix_destroy (X);
 740   pspp_linreg_cache_free (lcache);
 741   free (lopts.get_indep_mean_std);
 742   free (indep_vars);
 743   free (is_missing_case);
 744   casereader_destroy (r);
 745   return;
 746 }
 747
 748 /*
 749   Local Variables:
 750   mode: c
 751   End:
 752 */