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