pspp: Add --no-output option to allow entirely disabling output.
[pspp] / src / language / stats / regression.q
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2005, 2009 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
16
17 #include <config.h>
18
19 #include <gsl/gsl_cdf.h>
20 #include <gsl/gsl_matrix.h>
21 #include <gsl/gsl_vector.h>
22 #include <math.h>
23 #include <stdlib.h>
24 #include <data/case.h>
25 #include <data/casegrouper.h>
26 #include <data/casereader.h>
27 #include <data/dictionary.h>
28 #include <data/missing-values.h>
29 #include <data/procedure.h>
30 #include <data/transformations.h>
31 #include <data/value-labels.h>
32 #include <data/variable.h>
33 #include <language/command.h>
34 #include <language/dictionary/split-file.h>
35 #include <language/data-io/file-handle.h>
36 #include <language/lexer/lexer.h>
37 #include <libpspp/compiler.h>
38 #include <libpspp/message.h>
39 #include <libpspp/taint.h>
40 #include <math/covariance.h>
41 #include <math/linreg.h>
42 #include <math/moments.h>
43 #include <output/tab.h>
44
45 #include "xalloc.h"
46
47 #include "gettext.h"
48 #define _(msgid) gettext (msgid)
49
50 #define REG_LARGE_DATA 1000
51
52 /* (headers) */
53
54 /* (specification)
55    "REGRESSION" (regression_):
56    *variables=custom;
57    +statistics[st_]=r,
58                     coeff,
59                     anova,
60                     outs,
61                     zpp,
62                     label,
63                     sha,
64                     ci,
65                     bcov,
66                     ses,
67                     xtx,
68                     collin,
69                     tol,
70                     selection,
71                     f,
72                     defaults,
73                     all;
74    ^dependent=varlist;
75    +save[sv_]=resid,pred;
76    +method=enter.
77 */
78 /* (declarations) */
79 /* (functions) */
80 static struct cmd_regression cmd;
81
82 /*
83   Moments for each of the variables used.
84  */
85 struct moments_var
86 {
87   struct moments1 *m;
88   const struct variable *v;
89 };
90
91 /*
92   Transformations for saving predicted values
93   and residuals, etc.
94  */
95 struct reg_trns
96 {
97   int n_trns;                   /* Number of transformations. */
98   int trns_id;                  /* Which trns is this one? */
99   linreg *c;            /* Linear model for this trns. */
100 };
101 /*
102   Variables used (both explanatory and response).
103  */
104 static const struct variable **v_variables;
105
106 /*
107   Number of variables.
108  */
109 static size_t n_variables;
110
111 static bool run_regression (struct casereader *, struct cmd_regression *,
112                             struct dataset *, linreg **);
113
114 /*
115    STATISTICS subcommand output functions.
116  */
117 static void reg_stats_r (linreg *, void *);
118 static void reg_stats_coeff (linreg *, void *);
119 static void reg_stats_anova (linreg *, void *);
120 static void reg_stats_outs (linreg *, void *);
121 static void reg_stats_zpp (linreg *, void *);
122 static void reg_stats_label (linreg *, void *);
123 static void reg_stats_sha (linreg *, void *);
124 static void reg_stats_ci (linreg *, void *);
125 static void reg_stats_f (linreg *, void *);
126 static void reg_stats_bcov (linreg *, void *);
127 static void reg_stats_ses (linreg *, void *);
128 static void reg_stats_xtx (linreg *, void *);
129 static void reg_stats_collin (linreg *, void *);
130 static void reg_stats_tol (linreg *, void *);
131 static void reg_stats_selection (linreg *, void *);
132 static void statistics_keyword_output (void (*)(linreg *, void *),
133                                        int, linreg *, void *);
134
135 static void
136 reg_stats_r (linreg *c, void *aux UNUSED)
137 {
138   struct tab_table *t;
139   int n_rows = 2;
140   int n_cols = 5;
141   double rsq;
142   double adjrsq;
143   double std_error;
144
145   assert (c != NULL);
146   rsq = linreg_ssreg (c) / linreg_sst (c);
147   adjrsq = 1.0 - (1.0 - rsq) * (linreg_n_obs (c) - 1.0) / (linreg_n_obs (c) - linreg_n_coeffs (c));
148   std_error = sqrt (linreg_mse (c));
149   t = tab_create (n_cols, n_rows);
150   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
151   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
152   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
153   tab_vline (t, TAL_0, 1, 0, 0);
154
155   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
156   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
157   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
158   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
159   tab_double (t, 1, 1, TAB_RIGHT, sqrt (rsq), NULL);
160   tab_double (t, 2, 1, TAB_RIGHT, rsq, NULL);
161   tab_double (t, 3, 1, TAB_RIGHT, adjrsq, NULL);
162   tab_double (t, 4, 1, TAB_RIGHT, std_error, NULL);
163   tab_title (t, _("Model Summary"));
164   tab_submit (t);
165 }
166
167 /*
168   Table showing estimated regression coefficients.
169  */
170 static void
171 reg_stats_coeff (linreg * c, void *aux_)
172 {
173   size_t j;
174   int n_cols = 7;
175   int n_rows;
176   int this_row;
177   double t_stat;
178   double pval;
179   double std_err;
180   double beta;
181   const char *label;
182
183   const struct variable *v;
184   struct tab_table *t;
185   gsl_matrix *cov = aux_;
186
187   assert (c != NULL);
188   n_rows = linreg_n_coeffs (c) + 3;
189
190   t = tab_create (n_cols, n_rows);
191   tab_headers (t, 2, 0, 1, 0);
192   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
193   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
194   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
195   tab_vline (t, TAL_0, 1, 0, 0);
196
197   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
198   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
199   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
200   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
201   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
202   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
203   tab_double (t, 2, 1, 0, linreg_intercept (c), NULL);
204   std_err = sqrt (gsl_matrix_get (linreg_cov (c), 0, 0));
205   tab_double (t, 3, 1, 0, std_err, NULL);
206   tab_double (t, 4, 1, 0, 0.0, NULL);
207   t_stat = linreg_intercept (c) / std_err;
208   tab_double (t, 5, 1, 0, t_stat, NULL);
209   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
210   tab_double (t, 6, 1, 0, pval, NULL);
211   for (j = 0; j < linreg_n_coeffs (c); j++)
212     {
213       struct string tstr;
214       ds_init_empty (&tstr);
215       this_row = j + 2;
216
217       v = linreg_indep_var (c, j);
218       label = var_to_string (v);
219       /* Do not overwrite the variable's name. */
220       ds_put_cstr (&tstr, label);
221       tab_text (t, 1, this_row, TAB_CENTER, ds_cstr (&tstr));
222       /*
223          Regression coefficients.
224        */
225       tab_double (t, 2, this_row, 0, linreg_coeff (c, j), NULL);
226       /*
227          Standard error of the coefficients.
228        */
229       std_err = sqrt (gsl_matrix_get (linreg_cov (c), j + 1, j + 1));
230       tab_double (t, 3, this_row, 0, std_err, NULL);
231       /*
232          Standardized coefficient, i.e., regression coefficient
233          if all variables had unit variance.
234        */
235       beta = sqrt (gsl_matrix_get (cov, j, j));
236       beta *= linreg_coeff (c, j) / 
237         sqrt (gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1));
238       tab_double (t, 4, this_row, 0, beta, NULL);
239
240       /*
241          Test statistic for H0: coefficient is 0.
242        */
243       t_stat = linreg_coeff (c, j) / std_err;
244       tab_double (t, 5, this_row, 0, t_stat, NULL);
245       /*
246          P values for the test statistic above.
247        */
248       pval =
249         2 * gsl_cdf_tdist_Q (fabs (t_stat),
250                              (double) (linreg_n_obs (c) - linreg_n_coeffs (c)));
251       tab_double (t, 6, this_row, 0, pval, NULL);
252       ds_destroy (&tstr);
253     }
254   tab_title (t, _("Coefficients"));
255   tab_submit (t);
256 }
257
258 /*
259   Display the ANOVA table.
260  */
261 static void
262 reg_stats_anova (linreg * c, void *aux UNUSED)
263 {
264   int n_cols = 7;
265   int n_rows = 4;
266   const double msm = linreg_ssreg (c) / linreg_dfmodel (c);
267   const double mse = linreg_mse (c);
268   const double F = msm / mse;
269   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
270
271   struct tab_table *t;
272
273   assert (c != NULL);
274   t = tab_create (n_cols, n_rows);
275   tab_headers (t, 2, 0, 1, 0);
276
277   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
278
279   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
280   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
281   tab_vline (t, TAL_0, 1, 0, 0);
282
283   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
284   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
285   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
286   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
287   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
288
289   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
290   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
291   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
292
293   /* Sums of Squares */
294   tab_double (t, 2, 1, 0, linreg_ssreg (c), NULL);
295   tab_double (t, 2, 3, 0, linreg_sst (c), NULL);
296   tab_double (t, 2, 2, 0, linreg_sse (c), NULL);
297
298
299   /* Degrees of freedom */
300   tab_text_format (t, 3, 1, TAB_RIGHT, "%g", c->dfm);
301   tab_text_format (t, 3, 2, TAB_RIGHT, "%g", c->dfe);
302   tab_text_format (t, 3, 3, TAB_RIGHT, "%g", c->dft);
303
304   /* Mean Squares */
305   tab_double (t, 4, 1, TAB_RIGHT, msm, NULL);
306   tab_double (t, 4, 2, TAB_RIGHT, mse, NULL);
307
308   tab_double (t, 5, 1, 0, F, NULL);
309
310   tab_double (t, 6, 1, 0, pval, NULL);
311
312   tab_title (t, _("ANOVA"));
313   tab_submit (t);
314 }
315
316 static void
317 reg_stats_outs (linreg * c, void *aux UNUSED)
318 {
319   assert (c != NULL);
320 }
321
322 static void
323 reg_stats_zpp (linreg * c, void *aux UNUSED)
324 {
325   assert (c != NULL);
326 }
327
328 static void
329 reg_stats_label (linreg * c, void *aux UNUSED)
330 {
331   assert (c != NULL);
332 }
333
334 static void
335 reg_stats_sha (linreg * c, void *aux UNUSED)
336 {
337   assert (c != NULL);
338 }
339 static void
340 reg_stats_ci (linreg * c, void *aux UNUSED)
341 {
342   assert (c != NULL);
343 }
344 static void
345 reg_stats_f (linreg * c, void *aux UNUSED)
346 {
347   assert (c != NULL);
348 }
349 static void
350 reg_stats_bcov (linreg * c, void *aux UNUSED)
351 {
352   int n_cols;
353   int n_rows;
354   int i;
355   int k;
356   int row;
357   int col;
358   const char *label;
359   struct tab_table *t;
360
361   assert (c != NULL);
362   n_cols = c->n_indeps + 1 + 2;
363   n_rows = 2 * (c->n_indeps + 1);
364   t = tab_create (n_cols, n_rows);
365   tab_headers (t, 2, 0, 1, 0);
366   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
367   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
368   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
369   tab_vline (t, TAL_0, 1, 0, 0);
370   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
371   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
372   for (i = 0; i < linreg_n_coeffs (c); i++)
373     {
374       const struct variable *v = linreg_indep_var (c, i);
375       label = var_to_string (v);
376       tab_text (t, 2, i, TAB_CENTER, label);
377       tab_text (t, i + 2, 0, TAB_CENTER, label);
378       for (k = 1; k < linreg_n_coeffs (c); k++)
379         {
380           col = (i <= k) ? k : i;
381           row = (i <= k) ? i : k;
382           tab_double (t, k + 2, i, TAB_CENTER,
383                      gsl_matrix_get (c->cov, row, col), NULL);
384         }
385     }
386   tab_title (t, _("Coefficient Correlations"));
387   tab_submit (t);
388 }
389 static void
390 reg_stats_ses (linreg * c, void *aux UNUSED)
391 {
392   assert (c != NULL);
393 }
394 static void
395 reg_stats_xtx (linreg * c, void *aux UNUSED)
396 {
397   assert (c != NULL);
398 }
399 static void
400 reg_stats_collin (linreg * c, void *aux UNUSED)
401 {
402   assert (c != NULL);
403 }
404 static void
405 reg_stats_tol (linreg * c, void *aux UNUSED)
406 {
407   assert (c != NULL);
408 }
409 static void
410 reg_stats_selection (linreg * c, void *aux UNUSED)
411 {
412   assert (c != NULL);
413 }
414
415 static void
416 statistics_keyword_output (void (*function) (linreg *, void *),
417                            int keyword, linreg * c, void *aux)
418 {
419   if (keyword)
420     {
421       (*function) (c, aux);
422     }
423 }
424
425 static void
426 subcommand_statistics (int *keywords, linreg * c, void *aux)
427 {
428   /*
429      The order here must match the order in which the STATISTICS
430      keywords appear in the specification section above.
431    */
432   enum
433   { r,
434     coeff,
435     anova,
436     outs,
437     zpp,
438     label,
439     sha,
440     ci,
441     bcov,
442     ses,
443     xtx,
444     collin,
445     tol,
446     selection,
447     f,
448     defaults,
449     all
450   };
451   int i;
452   int d = 1;
453
454   if (keywords[all])
455     {
456       /*
457          Set everything but F.
458        */
459       for (i = 0; i < f; i++)
460         {
461           keywords[i] = 1;
462         }
463     }
464   else
465     {
466       for (i = 0; i < all; i++)
467         {
468           if (keywords[i])
469             {
470               d = 0;
471             }
472         }
473       /*
474          Default output: ANOVA table, parameter estimates,
475          and statistics for variables not entered into model,
476          if appropriate.
477        */
478       if (keywords[defaults] | d)
479         {
480           keywords[anova] = 1;
481           keywords[outs] = 1;
482           keywords[coeff] = 1;
483           keywords[r] = 1;
484         }
485     }
486   statistics_keyword_output (reg_stats_r, keywords[r], c, aux);
487   statistics_keyword_output (reg_stats_anova, keywords[anova], c, aux);
488   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c, aux);
489   statistics_keyword_output (reg_stats_outs, keywords[outs], c, aux);
490   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c, aux);
491   statistics_keyword_output (reg_stats_label, keywords[label], c, aux);
492   statistics_keyword_output (reg_stats_sha, keywords[sha], c, aux);
493   statistics_keyword_output (reg_stats_ci, keywords[ci], c, aux);
494   statistics_keyword_output (reg_stats_f, keywords[f], c, aux);
495   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c, aux);
496   statistics_keyword_output (reg_stats_ses, keywords[ses], c, aux);
497   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c, aux);
498   statistics_keyword_output (reg_stats_collin, keywords[collin], c, aux);
499   statistics_keyword_output (reg_stats_tol, keywords[tol], c, aux);
500   statistics_keyword_output (reg_stats_selection, keywords[selection], c, aux);
501 }
502
503 /*
504   Free the transformation. Free its linear model if this
505   transformation is the last one.
506  */
507 static bool
508 regression_trns_free (void *t_)
509 {
510   bool result = true;
511   struct reg_trns *t = t_;
512
513   if (t->trns_id == t->n_trns)
514     {
515       result = linreg_free (t->c);
516     }
517   free (t);
518
519   return result;
520 }
521
522 /*
523   Gets the predicted values.
524  */
525 static int
526 regression_trns_pred_proc (void *t_, struct ccase **c,
527                            casenumber case_idx UNUSED)
528 {
529   size_t i;
530   size_t n_vals;
531   struct reg_trns *trns = t_;
532   linreg *model;
533   union value *output = NULL;
534   const union value *tmp;
535   double *vals;
536   const struct variable **vars = NULL;
537
538   assert (trns != NULL);
539   model = trns->c;
540   assert (model != NULL);
541   assert (model->depvar != NULL);
542   assert (model->pred != NULL);
543
544   vars = linreg_get_vars (model);
545   n_vals = linreg_n_coeffs (model);
546   vals = xnmalloc (n_vals, sizeof (*vals));
547   *c = case_unshare (*c);
548
549   output = case_data_rw (*c, model->pred);
550
551   for (i = 0; i < n_vals; i++)
552     {
553       tmp = case_data (*c, vars[i]);
554       vals[i] = tmp->f;
555     }
556   output->f = linreg_predict (model, vals, n_vals);
557   free (vals);
558   return TRNS_CONTINUE;
559 }
560
561 /*
562   Gets the residuals.
563  */
564 static int
565 regression_trns_resid_proc (void *t_, struct ccase **c,
566                             casenumber case_idx UNUSED)
567 {
568   size_t i;
569   size_t n_vals;
570   struct reg_trns *trns = t_;
571   linreg *model;
572   union value *output = NULL;
573   const union value *tmp;
574   double *vals = NULL;
575   double obs;
576   const struct variable **vars = NULL;
577
578   assert (trns != NULL);
579   model = trns->c;
580   assert (model != NULL);
581   assert (model->depvar != NULL);
582   assert (model->resid != NULL);
583
584   vars = linreg_get_vars (model);
585   n_vals = linreg_n_coeffs (model);
586
587   vals = xnmalloc (n_vals, sizeof (*vals));
588   *c = case_unshare (*c);
589   output = case_data_rw (*c, model->resid);
590   assert (output != NULL);
591
592   for (i = 0; i < n_vals; i++)
593     {
594       tmp = case_data (*c, vars[i]);
595       vals[i] = tmp->f;
596     }
597   tmp = case_data (*c, model->depvar);
598   obs = tmp->f;
599   output->f = linreg_residual (model, obs, vals, n_vals);
600   free (vals);
601
602   return TRNS_CONTINUE;
603 }
604
605 /*
606    Returns false if NAME is a duplicate of any existing variable name.
607 */
608 static bool
609 try_name (const struct dictionary *dict, const char *name)
610 {
611   if (dict_lookup_var (dict, name) != NULL)
612     return false;
613
614   return true;
615 }
616
617 static void
618 reg_get_name (const struct dictionary *dict, char name[VAR_NAME_LEN],
619               const char prefix[VAR_NAME_LEN])
620 {
621   int i = 1;
622
623   snprintf (name, VAR_NAME_LEN, "%s%d", prefix, i);
624   while (!try_name (dict, name))
625     {
626       i++;
627       snprintf (name, VAR_NAME_LEN, "%s%d", prefix, i);
628     }
629 }
630
631 static void
632 reg_save_var (struct dataset *ds, const char *prefix, trns_proc_func * f,
633               linreg * c, struct variable **v, int n_trns)
634 {
635   struct dictionary *dict = dataset_dict (ds);
636   static int trns_index = 1;
637   char name[VAR_NAME_LEN];
638   struct variable *new_var;
639   struct reg_trns *t = NULL;
640
641   t = xmalloc (sizeof (*t));
642   t->trns_id = trns_index;
643   t->n_trns = n_trns;
644   t->c = c;
645   reg_get_name (dict, name, prefix);
646   new_var = dict_create_var (dict, name, 0);
647   assert (new_var != NULL);
648   *v = new_var;
649   add_transformation (ds, f, regression_trns_free, t);
650   trns_index++;
651 }
652 static void
653 subcommand_save (struct dataset *ds, int save, linreg ** models)
654 {
655   linreg **lc;
656   int n_trns = 0;
657   int i;
658
659   if (save)
660     {
661       /* Count the number of transformations we will need. */
662       for (i = 0; i < REGRESSION_SV_count; i++)
663         {
664           if (cmd.a_save[i])
665             {
666               n_trns++;
667             }
668         }
669       n_trns *= cmd.n_dependent;
670
671       for (lc = models; lc < models + cmd.n_dependent; lc++)
672         {
673           if (*lc != NULL)
674             {
675               if ((*lc)->depvar != NULL)
676                 {
677                   if (cmd.a_save[REGRESSION_SV_RESID])
678                     {
679                       reg_save_var (ds, "RES", regression_trns_resid_proc, *lc,
680                                     &(*lc)->resid, n_trns);
681                     }
682                   if (cmd.a_save[REGRESSION_SV_PRED])
683                     {
684                       reg_save_var (ds, "PRED", regression_trns_pred_proc, *lc,
685                                     &(*lc)->pred, n_trns);
686                     }
687                 }
688             }
689         }
690     }
691   else
692     {
693       for (lc = models; lc < models + cmd.n_dependent; lc++)
694         {
695           if (*lc != NULL)
696             {
697               linreg_free (*lc);
698             }
699         }
700     }
701 }
702
703 int
704 cmd_regression (struct lexer *lexer, struct dataset *ds)
705 {
706   struct casegrouper *grouper;
707   struct casereader *group;
708   linreg **models;
709   bool ok;
710   size_t i;
711
712   if (!parse_regression (lexer, ds, &cmd, NULL))
713     {
714       return CMD_FAILURE;
715     }
716
717   models = xnmalloc (cmd.n_dependent, sizeof *models);
718   for (i = 0; i < cmd.n_dependent; i++)
719     {
720       models[i] = NULL;
721     }
722
723   /* Data pass. */
724   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
725   while (casegrouper_get_next_group (grouper, &group))
726     run_regression (group, &cmd, ds, models);
727   ok = casegrouper_destroy (grouper);
728   ok = proc_commit (ds) && ok;
729
730   subcommand_save (ds, cmd.sbc_save, models);
731   free (v_variables);
732   free (models);
733   free_regression (&cmd);
734
735   return ok ? CMD_SUCCESS : CMD_FAILURE;
736 }
737
738 /*
739   Is variable k the dependent variable?
740  */
741 static bool
742 is_depvar (size_t k, const struct variable *v)
743 {
744   return v == v_variables[k];
745 }
746
747 /* Parser for the variables sub command */
748 static int
749 regression_custom_variables (struct lexer *lexer, struct dataset *ds,
750                              struct cmd_regression *cmd UNUSED,
751                              void *aux UNUSED)
752 {
753   const struct dictionary *dict = dataset_dict (ds);
754
755   lex_match (lexer, '=');
756
757   if ((lex_token (lexer) != T_ID
758        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
759       && lex_token (lexer) != T_ALL)
760     return 2;
761
762
763   if (!parse_variables_const
764       (lexer, dict, &v_variables, &n_variables, PV_NONE))
765     {
766       free (v_variables);
767       return 0;
768     }
769   assert (n_variables);
770
771   return 1;
772 }
773
774 /* Identify the explanatory variables in v_variables.  Returns
775    the number of independent variables. */
776 static int
777 identify_indep_vars (const struct variable **indep_vars,
778                      const struct variable *depvar)
779 {
780   int n_indep_vars = 0;
781   int i;
782
783   for (i = 0; i < n_variables; i++)
784     if (!is_depvar (i, depvar))
785       indep_vars[n_indep_vars++] = v_variables[i];
786   if ((n_indep_vars < 1) && is_depvar (0, depvar))
787     {
788       /*
789         There is only one independent variable, and it is the same
790         as the dependent variable. Print a warning and continue.
791        */
792       msg (SE,
793            gettext ("The dependent variable is equal to the independent variable." 
794                     "The least squares line is therefore Y=X." 
795                     "Standard errors and related statistics may be meaningless."));
796       n_indep_vars = 1;
797       indep_vars[0] = v_variables[0];
798     }
799   return n_indep_vars;
800 }
801 static double
802 fill_covariance (gsl_matrix *cov, struct covariance *all_cov, 
803                  const struct variable **vars,
804                  size_t n_vars, const struct variable *dep_var, 
805                  const struct variable **all_vars, size_t n_all_vars,
806                  double *means)
807 {
808   size_t i;
809   size_t j;
810   size_t dep_subscript;
811   size_t *rows;
812   const gsl_matrix *ssizes;
813   const gsl_matrix *cm;
814   const gsl_matrix *mean_matrix;
815   const gsl_matrix *ssize_matrix;
816   double result = 0.0;
817   
818   cm = covariance_calculate_unnormalized (all_cov);
819   rows = xnmalloc (cov->size1 - 1, sizeof (*rows));
820   
821   for (i = 0; i < n_all_vars; i++)
822     {
823       for (j = 0; j < n_vars; j++)
824         {
825           if (vars[j] == all_vars[i])
826             {
827               rows[j] = i;
828             }
829         }
830       if (all_vars[i] == dep_var)
831         {
832           dep_subscript = i;
833         }
834     }
835   mean_matrix = covariance_moments (all_cov, MOMENT_MEAN);
836   ssize_matrix = covariance_moments (all_cov, MOMENT_NONE);
837   for (i = 0; i < cov->size1 - 1; i++)
838     {
839       means[i] = gsl_matrix_get (mean_matrix, rows[i], 0)
840         / gsl_matrix_get (ssize_matrix, rows[i], 0);
841       for (j = 0; j < cov->size2 - 1; j++)
842         {
843           gsl_matrix_set (cov, i, j, gsl_matrix_get (cm, rows[i], rows[j]));
844           gsl_matrix_set (cov, j, i, gsl_matrix_get (cm, rows[j], rows[i]));
845         }
846     }
847   means[cov->size1 - 1] = gsl_matrix_get (mean_matrix, dep_subscript, 0)
848     / gsl_matrix_get (ssize_matrix, dep_subscript, 0);
849   ssizes = covariance_moments (all_cov, MOMENT_NONE);
850   result = gsl_matrix_get (ssizes, dep_subscript, rows[0]);
851   for (i = 0; i < cov->size1 - 1; i++)
852     {
853       gsl_matrix_set (cov, i, cov->size1 - 1, 
854                       gsl_matrix_get (cm, rows[i], dep_subscript));
855       gsl_matrix_set (cov, cov->size1 - 1, i, 
856                       gsl_matrix_get (cm, rows[i], dep_subscript));
857       if (result > gsl_matrix_get (ssizes, rows[i], dep_subscript))
858         {
859           result = gsl_matrix_get (ssizes, rows[i], dep_subscript);
860         }
861     }
862   gsl_matrix_set (cov, cov->size1 - 1, cov->size1 - 1, 
863                   gsl_matrix_get (cm, dep_subscript, dep_subscript));
864   free (rows);
865   return result;
866 }
867
868 static bool
869 run_regression (struct casereader *input, struct cmd_regression *cmd,
870                 struct dataset *ds, linreg **models)
871 {
872   size_t i;
873   int n_indep = 0;
874   int k;
875   double n_data;
876   double *means;
877   struct ccase *c;
878   struct covariance *cov;
879   const struct variable **vars;
880   const struct variable *dep_var;
881   struct casereader *reader;
882   const struct dictionary *dict;
883   gsl_matrix *this_cm;
884
885   assert (models != NULL);
886
887   for (i = 0; i < n_variables; i++)
888     {
889       if (!var_is_numeric (v_variables[i]))
890         {
891           msg (SE, _("REGRESSION requires numeric variables."));
892           return false;
893         }
894     }
895
896   c = casereader_peek (input, 0);
897   if (c == NULL)
898     {
899       casereader_destroy (input);
900       return true;
901     }
902   output_split_file_values (ds, c);
903   case_unref (c);
904
905   dict = dataset_dict (ds);
906   if (!v_variables)
907     {
908       dict_get_vars (dict, &v_variables, &n_variables, 0);
909     }
910   vars = xnmalloc (n_variables, sizeof (*vars));
911   means  = xnmalloc (n_variables, sizeof (*means));
912   cov = covariance_1pass_create (n_variables, v_variables,
913                                  dict_get_weight (dict), MV_ANY);
914
915   reader = casereader_clone (input);
916   reader = casereader_create_filter_missing (reader, v_variables, n_variables,
917                                              MV_ANY, NULL, NULL);
918   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
919     {
920       covariance_accumulate (cov, c);
921     }
922   
923   for (k = 0; k < cmd->n_dependent; k++)
924     {
925       dep_var = cmd->v_dependent[k];
926       n_indep = identify_indep_vars (vars, dep_var);
927       
928       this_cm = gsl_matrix_alloc (n_indep + 1, n_indep + 1);
929       n_data = fill_covariance (this_cm, cov, vars, n_indep, 
930                                 dep_var, v_variables, n_variables, means);
931       models[k] = linreg_alloc (dep_var, (const struct variable **) vars,
932                                 n_data, n_indep);
933       models[k]->depvar = dep_var;
934       for (i = 0; i < n_indep; i++)
935         {
936           linreg_set_indep_variable_mean (models[k], i, means[i]);
937         }
938       linreg_set_depvar_mean (models[k], means[i]);
939       /*
940         For large data sets, use QR decomposition.
941       */
942       if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
943         {
944           models[k]->method = LINREG_QR;
945         }
946       
947       if (n_data > 0)
948         {
949           /*
950             Find the least-squares estimates and other statistics.
951           */
952           linreg_fit (this_cm, models[k]);
953           
954           if (!taint_has_tainted_successor (casereader_get_taint (input)))
955             {
956               subcommand_statistics (cmd->a_statistics, models[k], this_cm);
957             }
958         }
959       else
960         {
961           msg (SE,
962                gettext ("No valid data found. This command was skipped."));
963           linreg_free (models[k]);
964           models[k] = NULL;
965         }
966     }
967   
968   casereader_destroy (reader);
969   free (vars);
970   free (means);
971   casereader_destroy (input);
972   covariance_destroy (cov);
973   
974   return true;
975 }
976
977 /*
978   Local Variables:
979   mode: c
980   End:
981 */