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
  48 /* (specification)
  49    "REGRESSION" (regression_):
  50    *variables=varlist;
  51    statistics[st_]=r,
  52    coeff,
  53    anova,
  54    outs,
  55    zpp,
  56    label,
  57    sha,
  58    ci,
  59    bcov,
  60    ses,
  61    xtx,
  62    collin,
  63    tol,
  64    selection,
  65    f,
  66    defaults,
  67    all;
  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   Return value for the procedure.
  82  */
  83 int pspp_reg_rc = CMD_SUCCESS;
  84
  85 static void run_regression (const struct casefile *, void *);
  86 /*
  87    STATISTICS subcommand output functions.
  88  */
  89 static void reg_stats_r (pspp_linreg_cache *);
  90 static void reg_stats_coeff (pspp_linreg_cache *);
  91 static void reg_stats_anova (pspp_linreg_cache *);
  92 static void reg_stats_outs (pspp_linreg_cache *);
  93 static void reg_stats_zpp (pspp_linreg_cache *);
  94 static void reg_stats_label (pspp_linreg_cache *);
  95 static void reg_stats_sha (pspp_linreg_cache *);
  96 static void reg_stats_ci (pspp_linreg_cache *);
  97 static void reg_stats_f (pspp_linreg_cache *);
  98 static void reg_stats_bcov (pspp_linreg_cache *);
  99 static void reg_stats_ses (pspp_linreg_cache *);
 100 static void reg_stats_xtx (pspp_linreg_cache *);
 101 static void reg_stats_collin (pspp_linreg_cache *);
 102 static void reg_stats_tol (pspp_linreg_cache *);
 103 static void reg_stats_selection (pspp_linreg_cache *);
 104 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
 105                                        int, pspp_linreg_cache *);
 106
 107 static void
 108 reg_stats_r (pspp_linreg_cache * c)
 109 {
 110   struct tab_table *t;
 111   int n_rows = 2;
 112   int n_cols = 5;
 113   double rsq;
 114   double adjrsq;
 115   double std_error;
 116
 117   assert (c != NULL);
 118   rsq = c->ssm / c->sst;
 119   adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
 120   std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
 121   t = tab_create (n_cols, n_rows, 0);
 122   tab_dim (t, tab_natural_dimensions);
 123   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 124   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 125   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 126   tab_vline (t, TAL_0, 1, 0, 0);
 127
 128   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
 129   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
 130   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
 131   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
 132   tab_float (t, 1, 1, TAB_RIGHT, sqrt (rsq), 10, 2);
 133   tab_float (t, 2, 1, TAB_RIGHT, rsq, 10, 2);
 134   tab_float (t, 3, 1, TAB_RIGHT, adjrsq, 10, 2);
 135   tab_float (t, 4, 1, TAB_RIGHT, std_error, 10, 2);
 136   tab_title (t, 0, _("Model Summary"));
 137   tab_submit (t);
 138 }
 139
 140 /*
 141   Table showing estimated regression coefficients.
 142  */
 143 static void
 144 reg_stats_coeff (pspp_linreg_cache * c)
 145 {
 146   size_t i;
 147   size_t j;
 148   int n_cols = 7;
 149   int n_rows;
 150   double t_stat;
 151   double pval;
 152   double coeff;
 153   double std_err;
 154   double beta;
 155   const char *label;
 156   struct tab_table *t;
 157
 158   assert (c != NULL);
 159   n_rows = c->n_coeffs + 2;
 160
 161   t = tab_create (n_cols, n_rows, 0);
 162   tab_headers (t, 2, 0, 1, 0);
 163   tab_dim (t, tab_natural_dimensions);
 164   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 165   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 166   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 167   tab_vline (t, TAL_0, 1, 0, 0);
 168
 169   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
 170   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 171   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
 172   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
 173   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 174   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
 175   coeff = c->coeff[0].estimate;
 176   tab_float (t, 2, 1, 0, coeff, 10, 2);
 177   std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
 178   tab_float (t, 3, 1, 0, std_err, 10, 2);
 179   beta = coeff / c->depvar_std;
 180   tab_float (t, 4, 1, 0, beta, 10, 2);
 181   t_stat = coeff / std_err;
 182   tab_float (t, 5, 1, 0, t_stat, 10, 2);
 183   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 184   tab_float (t, 6, 1, 0, pval, 10, 2);
 185   for (j = 1; j <= c->n_indeps; j++)
 186     {
 187       i = indep_vars[j];
 188       label = var_to_string (c->coeff[j].v);
 189       tab_text (t, 1, j + 1, TAB_CENTER, label);
 190       /*
 191          Regression coefficients.
 192        */
 193       coeff = c->coeff[j].estimate;
 194       tab_float (t, 2, j + 1, 0, coeff, 10, 2);
 195       /*
 196          Standard error of the coefficients.
 197        */
 198       std_err = sqrt (gsl_matrix_get (c->cov, j, j));
 199       tab_float (t, 3, j + 1, 0, std_err, 10, 2);
 200       /*
 201          'Standardized' coefficient, i.e., regression coefficient
 202          if all variables had unit variance.
 203        */
 204       beta = gsl_vector_get (c->indep_std, j);
 205       beta *= coeff / c->depvar_std;
 206       tab_float (t, 4, j + 1, 0, beta, 10, 2);
 207
 208       /*
 209          Test statistic for H0: coefficient is 0.
 210        */
 211       t_stat = coeff / std_err;
 212       tab_float (t, 5, j + 1, 0, t_stat, 10, 2);
 213       /*
 214          P values for the test statistic above.
 215        */
 216       pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 217       tab_float (t, 6, j + 1, 0, pval, 10, 2);
 218     }
 219   tab_title (t, 0, _("Coefficients"));
 220   tab_submit (t);
 221 }
 222
 223 /*
 224   Display the ANOVA table.
 225  */
 226 static void
 227 reg_stats_anova (pspp_linreg_cache * c)
 228 {
 229   int n_cols = 7;
 230   int n_rows = 4;
 231   const double msm = c->ssm / c->dfm;
 232   const double mse = c->sse / c->dfe;
 233   const double F = msm / mse;
 234   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
 235
 236   struct tab_table *t;
 237
 238   assert (c != NULL);
 239   t = tab_create (n_cols, n_rows, 0);
 240   tab_headers (t, 2, 0, 1, 0);
 241   tab_dim (t, tab_natural_dimensions);
 242
 243   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 244
 245   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 246   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 247   tab_vline (t, TAL_0, 1, 0, 0);
 248
 249   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 250   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 251   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 252   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 253   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 254
 255   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
 256   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
 257   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
 258
 259   /* Sums of Squares */
 260   tab_float (t, 2, 1, 0, c->ssm, 10, 2);
 261   tab_float (t, 2, 3, 0, c->sst, 10, 2);
 262   tab_float (t, 2, 2, 0, c->sse, 10, 2);
 263
 264
 265   /* Degrees of freedom */
 266   tab_float (t, 3, 1, 0, c->dfm, 4, 0);
 267   tab_float (t, 3, 2, 0, c->dfe, 4, 0);
 268   tab_float (t, 3, 3, 0, c->dft, 4, 0);
 269
 270   /* Mean Squares */
 271
 272   tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
 273   tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
 274
 275   tab_float (t, 5, 1, 0, F, 8, 3);
 276
 277   tab_float (t, 6, 1, 0, pval, 8, 3);
 278
 279   tab_title (t, 0, _("ANOVA"));
 280   tab_submit (t);
 281 }
 282 static void
 283 reg_stats_outs (pspp_linreg_cache * c)
 284 {
 285   assert (c != NULL);
 286 }
 287 static void
 288 reg_stats_zpp (pspp_linreg_cache * c)
 289 {
 290   assert (c != NULL);
 291 }
 292 static void
 293 reg_stats_label (pspp_linreg_cache * c)
 294 {
 295   assert (c != NULL);
 296 }
 297 static void
 298 reg_stats_sha (pspp_linreg_cache * c)
 299 {
 300   assert (c != NULL);
 301 }
 302 static void
 303 reg_stats_ci (pspp_linreg_cache * c)
 304 {
 305   assert (c != NULL);
 306 }
 307 static void
 308 reg_stats_f (pspp_linreg_cache * c)
 309 {
 310   assert (c != NULL);
 311 }
 312 static void
 313 reg_stats_bcov (pspp_linreg_cache * c)
 314 {
 315   int n_cols;
 316   int n_rows;
 317   int i;
 318   int j;
 319   int k;
 320   int row;
 321   int col;
 322   const char *label;
 323   struct tab_table *t;
 324
 325   assert (c != NULL);
 326   n_cols = c->n_indeps + 1 + 2;
 327   n_rows = 2 * (c->n_indeps + 1);
 328   t = tab_create (n_cols, n_rows, 0);
 329   tab_headers (t, 2, 0, 1, 0);
 330   tab_dim (t, tab_natural_dimensions);
 331   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 332   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 333   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 334   tab_vline (t, TAL_0, 1, 0, 0);
 335   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 336   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 337   for (i = 1; i < c->n_indeps + 1; i++)
 338     {
 339       j = indep_vars[(i - 1)];
 340       struct variable *v = cmd.v_variables[j];
 341       label = var_to_string (v);
 342       tab_text (t, 2, i, TAB_CENTER, label);
 343       tab_text (t, i + 2, 0, TAB_CENTER, label);
 344       for (k = 1; k < c->n_indeps + 1; k++)
 345         {
 346           col = (i <= k) ? k : i;
 347           row = (i <= k) ? i : k;
 348           tab_float (t, k + 2, i, TAB_CENTER,
 349                      gsl_matrix_get (c->cov, row, col), 8, 3);
 350         }
 351     }
 352   tab_title (t, 0, _("Coefficient Correlations"));
 353   tab_submit (t);
 354 }
 355 static void
 356 reg_stats_ses (pspp_linreg_cache * c)
 357 {
 358   assert (c != NULL);
 359 }
 360 static void
 361 reg_stats_xtx (pspp_linreg_cache * c)
 362 {
 363   assert (c != NULL);
 364 }
 365 static void
 366 reg_stats_collin (pspp_linreg_cache * c)
 367 {
 368   assert (c != NULL);
 369 }
 370 static void
 371 reg_stats_tol (pspp_linreg_cache * c)
 372 {
 373   assert (c != NULL);
 374 }
 375 static void
 376 reg_stats_selection (pspp_linreg_cache * c)
 377 {
 378   assert (c != NULL);
 379 }
 380
 381 static void
 382 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 383                            int keyword, pspp_linreg_cache * c)
 384 {
 385   if (keyword)
 386     {
 387       (*function) (c);
 388     }
 389 }
 390
 391 static void
 392 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 393 {
 394   /*
 395      The order here must match the order in which the STATISTICS
 396      keywords appear in the specification section above.
 397    */
 398   enum
 399   { r,
 400     coeff,
 401     anova,
 402     outs,
 403     zpp,
 404     label,
 405     sha,
 406     ci,
 407     bcov,
 408     ses,
 409     xtx,
 410     collin,
 411     tol,
 412     selection,
 413     f,
 414     defaults,
 415     all
 416   };
 417   int i;
 418   int d = 1;
 419
 420   if (keywords[all])
 421     {
 422       /*
 423          Set everything but F.
 424        */
 425       for (i = 0; i < f; i++)
 426         {
 427           keywords[i] = 1;
 428         }
 429     }
 430   else
 431     {
 432       for (i = 0; i < all; i++)
 433         {
 434           if (keywords[i])
 435             {
 436               d = 0;
 437             }
 438         }
 439       /*
 440          Default output: ANOVA table, parameter estimates,
 441          and statistics for variables not entered into model,
 442          if appropriate.
 443        */
 444       if (keywords[defaults] | d)
 445         {
 446           keywords[anova] = 1;
 447           keywords[outs] = 1;
 448           keywords[coeff] = 1;
 449           keywords[r] = 1;
 450         }
 451     }
 452   statistics_keyword_output (reg_stats_r, keywords[r], c);
 453   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 454   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 455   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 456   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 457   statistics_keyword_output (reg_stats_label, keywords[label], c);
 458   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 459   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 460   statistics_keyword_output (reg_stats_f, keywords[f], c);
 461   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 462   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 463   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 464   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 465   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 466   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 467 }
 468
 469 int
 470 cmd_regression (void)
 471 {
 472   if (!parse_regression (&cmd))
 473     {
 474       return CMD_FAILURE;
 475     }
 476   multipass_procedure_with_splits (run_regression, &cmd);
 477
 478   return pspp_reg_rc;
 479 }
 480
 481 /*
 482   Is variable k one of the dependent variables?
 483  */
 484 static int
 485 is_depvar (size_t k)
 486 {
 487   size_t j = 0;
 488   for (j = 0; j < cmd.n_dependent; j++)
 489     {
 490       /*
 491          compare_var_names returns 0 if the variable
 492          names match.
 493        */
 494       if (!compare_var_names (cmd.v_dependent[j], cmd.v_variables[k], NULL))
 495         return 1;
 496     }
 497   return 0;
 498 }
 499
 500 static void
 501 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
 502 {
 503   size_t i;
 504   size_t n_data = 0;
 505   size_t row;
 506   size_t case_num;
 507   int n_indep;
 508   int j = 0;
 509   /*
 510      Keep track of the missing cases.
 511    */
 512   int *is_missing_case;
 513   const union value *val;
 514   struct casereader *r;
 515   struct casereader *r2;
 516   struct ccase c;
 517   struct variable *v;
 518   struct variable **indep_vars;
 519   struct design_matrix *X;
 520   gsl_vector *Y;
 521   pspp_linreg_cache *lcache;
 522   pspp_linreg_opts lopts;
 523
 524   n_data = casefile_get_case_cnt (cf);
 525
 526   is_missing_case = xnmalloc (n_data, sizeof (*is_missing_case));
 527   for (i = 0; i < n_data; i++)
 528     is_missing_case[i] = 0;
 529
 530   n_indep = cmd.n_variables - cmd.n_dependent;
 531   indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
 532
 533   lopts.get_depvar_mean_std = 1;
 534   lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
 535
 536
 537   /*
 538      Read from the active file. The first pass encodes categorical
 539      variables and drops cases with missing values.
 540    */
 541   j = 0;
 542   for (i = 0; i < cmd.n_variables; i++)
 543     {
 544       if (!is_depvar (i))
 545         {
 546           v = cmd.v_variables[i];
 547           indep_vars[j] = v;
 548           j++;
 549           if (v->type == ALPHA)
 550             {
 551               /* Make a place to hold the binary vectors
 552                  corresponding to this variable's values. */
 553               cat_stored_values_create (v);
 554             }
 555           for (r = casefile_get_reader (cf);
 556                casereader_read (r, &c); case_destroy (&c))
 557             {
 558               row = casereader_cnum (r) - 1;
 559
 560               val = case_data (&c, v->fv);
 561               cat_value_update (v, val);
 562               if (mv_is_value_missing (&v->miss, val))
 563                 {
 564                   if (!is_missing_case[row])
 565                     {
 566                       /* Now it is missing. */
 567                       n_data--;
 568                       is_missing_case[row] = 1;
 569                     }
 570                 }
 571             }
 572         }
 573     }
 574
 575   Y = gsl_vector_alloc (n_data);
 576   X =
 577     design_matrix_create (n_indep, (const struct variable **) indep_vars,
 578                           n_data);
 579   lcache = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
 580   lcache->indep_means = gsl_vector_alloc (X->m->size2);
 581   lcache->indep_std = gsl_vector_alloc (X->m->size2);
 582
 583   /*
 584      The second pass creates the design matrix.
 585    */
 586   row = 0;
 587   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
 588        case_destroy (&c))
 589     /* Iterate over the cases. */
 590     {
 591       case_num = casereader_cnum (r2) - 1;
 592       if (!is_missing_case[case_num])
 593         {
 594           for (i = 0; i < cmd.n_variables; ++i) /* Iterate over the variables
 595                                                    for the current case.
 596                                                  */
 597             {
 598               v = cmd.v_variables[i];
 599               val = case_data (&c, v->fv);
 600               /*
 601                  Independent/dependent variable separation. The
 602                  'variables' subcommand specifies a varlist which contains
 603                  both dependent and independent variables. The dependent
 604                  variables are specified with the 'dependent'
 605                  subcommand. We need to separate the two.
 606                */
 607               if (is_depvar (i))
 608                 {
 609                   if (v->type != NUMERIC)
 610                     {
 611                       msg (SE,
 612                            gettext ("Dependent variable must be numeric."));
 613                       pspp_reg_rc = CMD_FAILURE;
 614                       return;
 615                     }
 616                   lcache->depvar = (const struct variable *) v;
 617                   gsl_vector_set (Y, row, val->f);
 618                 }
 619               else
 620                 {
 621                   if (v->type == ALPHA)
 622                     {
 623                       design_matrix_set_categorical (X, row, v, val);
 624                     }
 625                   else if (v->type == NUMERIC)
 626                     {
 627                       design_matrix_set_numeric (X, row, v, val);
 628                     }
 629
 630                   lopts.get_indep_mean_std[i] = 1;
 631                 }
 632             }
 633           row++;
 634         }
 635     }
 636   /*
 637      Now that we know the number of coefficients, allocate space
 638      and store pointers to the variables that correspond to the
 639      coefficients.
 640    */
 641   lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
 642   for (i = 0; i < X->m->size2; i++)
 643     {
 644       j = i + 1;                /* The first coeff is the intercept. */
 645       lcache->coeff[j].v =
 646         (const struct variable *) design_matrix_col_to_var (X, i);
 647       assert (lcache->coeff[j].v != NULL);
 648     }
 649   /*
 650     For large data sets, use QR decomposition.
 651    */
 652   if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
 653     {
 654       lcache->method = PSPP_LINREG_SVD;
 655     }
 656   /*
 657      Find the least-squares estimates and other statistics.
 658    */
 659     pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
 660   subcommand_statistics (cmd.a_statistics, lcache);
 661   gsl_vector_free (Y);
 662   design_matrix_destroy (X);
 663   pspp_linreg_cache_free (lcache);
 664   free (lopts.get_indep_mean_std);
 665   free (indep_vars);
 666   free (is_missing_case);
 667   casereader_destroy (r);
 668   return;
 669 }
 670
 671 /*
 672   Local Variables:
 673   mode: c
 674   End:
 675 */