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