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