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