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