28289622aaded41975b27df35eda4830b4d8c2bf
[pspp] / 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 */