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