0f8e8532a13ea2b55a997952073a8515aff737ae
[pspp-builds.git] / 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 k;
503   size_t n_data = 0;
504   size_t row;
505   size_t case_num;
506   int n_indep;
507   int j = 0;
508   /*
509      Keep track of the missing cases.
510    */
511   int *is_missing_case;
512   const union value *val;
513   struct casereader *r;
514   struct casereader *r2;
515   struct ccase c;
516   struct variable *v;
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   for (i = 0; i < cmd.n_variables; i++)
540     {
541       v = cmd.v_variables[i];
542       if (v->type == ALPHA)
543         {
544           /* Make a place to hold the binary vectors 
545              corresponding to this variable's values. */
546           cat_stored_values_create (v);
547         }
548       for (r = casefile_get_reader (cf);
549            casereader_read (r, &c); case_destroy (&c))
550         {
551           row = casereader_cnum (r) - 1;
552
553           val = case_data (&c, v->fv);
554           cat_value_update (v, val);
555           if (mv_is_value_missing (&v->miss, val))
556             {
557               if (!is_missing_case[row])
558                 {
559                   /* Now it is missing. */
560                   n_data--;
561                   is_missing_case[row] = 1;
562                 }
563             }
564         }
565     }
566
567   Y = gsl_vector_alloc (n_data);
568   X =
569     design_matrix_create (n_indep, (const struct variable **) cmd.v_variables,
570                           n_data);
571   lcache = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
572   lcache->indep_means = gsl_vector_alloc (X->m->size2);
573   lcache->indep_std = gsl_vector_alloc (X->m->size2);
574
575   /*
576      The second pass creates the design matrix.
577    */
578   row = 0;
579   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
580        case_destroy (&c))
581     /* Iterate over the cases. */
582     {
583       k = 0;
584       case_num = casereader_cnum (r2) - 1;
585       if (!is_missing_case[case_num])
586         {
587           for (i = 0; i < cmd.n_variables; ++i) /* Iterate over the variables
588                                                    for the current case. 
589                                                  */
590             {
591               v = cmd.v_variables[i];
592               val = case_data (&c, v->fv);
593               /*
594                  Independent/dependent variable separation. The
595                  'variables' subcommand specifies a varlist which contains
596                  both dependent and independent variables. The dependent
597                  variables are specified with the 'dependent'
598                  subcommand. We need to separate the two.
599                */
600               if (is_depvar (i))
601                 {
602                   if (v->type != NUMERIC)
603                     {
604                       msg (SE,
605                            gettext ("Dependent variable must be numeric."));
606                       pspp_reg_rc = CMD_FAILURE;
607                       return;
608                     }
609                   lcache->depvar = (const struct variable *) v;
610                   gsl_vector_set (Y, row, val->f);
611                 }
612               else
613                 {
614                   if (v->type == ALPHA)
615                     {
616                       design_matrix_set_categorical (X, row, v, val);
617                     }
618                   else if (v->type == NUMERIC)
619                     {
620                       design_matrix_set_numeric (X, row, v, val);
621                     }
622
623                   indep_vars[k] = i;
624                   k++;
625                   lopts.get_indep_mean_std[i] = 1;
626                 }
627             }
628           row++;
629         }
630     }
631   /*
632      Now that we know the number of coefficients, allocate space
633      and store pointers to the variables that correspond to the
634      coefficients.
635    */
636   lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
637   for (i = 0; i < X->m->size2; i++)
638     {
639       j = i + 1;                /* The first coeff is the intercept. */
640       lcache->coeff[j].v =
641         (const struct variable *) design_matrix_col_to_var (X, i);
642       assert (lcache->coeff[j].v != NULL);
643     }
644   /* 
645      Find the least-squares estimates and other statistics.
646    */
647   pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
648   subcommand_statistics (cmd.a_statistics, lcache);
649   gsl_vector_free (Y);
650   design_matrix_destroy (X);
651   pspp_linreg_cache_free (lcache);
652   free (lopts.get_indep_mean_std);
653   free (indep_vars);
654   free (is_missing_case);
655   casereader_destroy (r);
656   return;
657 }
658
659 /*
660   Local Variables:   
661   mode: c
662   End:
663 */