Fixed column/variable lookup
[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 "tab.h"
37 #include "var.h"
38 #include "vfm.h"
39
40 /* (headers) */
41
42
43 /* (specification)
44    "REGRESSION" (regression_):
45    *variables=varlist;
46    statistics[st_]=r,
47    coeff,
48    anova,
49    outs,
50    zpp,
51    label,
52    sha,
53    ci,
54    bcov,
55    ses,
56    xtx,
57    collin,
58    tol,
59    selection,
60    f,
61    defaults,
62    all;
63    ^dependent=varlist;
64    ^method=enter.
65 */
66 /* (declarations) */
67 /* (functions) */
68 static struct cmd_regression cmd;
69
70 /*
71   Array holding the subscripts of the independent variables.
72  */
73 size_t *indep_vars;
74
75 /*
76   Return value for the procedure.
77  */
78 int pspp_reg_rc = CMD_SUCCESS;
79
80 static void run_regression (const struct casefile *, void *);
81 /* 
82    STATISTICS subcommand output functions.
83  */
84 static void reg_stats_r (pspp_linreg_cache *);
85 static void reg_stats_coeff (pspp_linreg_cache *);
86 static void reg_stats_anova (pspp_linreg_cache *);
87 static void reg_stats_outs (pspp_linreg_cache *);
88 static void reg_stats_zpp (pspp_linreg_cache *);
89 static void reg_stats_label (pspp_linreg_cache *);
90 static void reg_stats_sha (pspp_linreg_cache *);
91 static void reg_stats_ci (pspp_linreg_cache *);
92 static void reg_stats_f (pspp_linreg_cache *);
93 static void reg_stats_bcov (pspp_linreg_cache *);
94 static void reg_stats_ses (pspp_linreg_cache *);
95 static void reg_stats_xtx (pspp_linreg_cache *);
96 static void reg_stats_collin (pspp_linreg_cache *);
97 static void reg_stats_tol (pspp_linreg_cache *);
98 static void reg_stats_selection (pspp_linreg_cache *);
99 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
100                                        int, pspp_linreg_cache *);
101
102 static void
103 reg_stats_r (pspp_linreg_cache * c)
104 {
105   struct tab_table *t;
106   int n_rows = 2;
107   int n_cols = 5;
108   double rsq;
109   double adjrsq;
110   double std_error;
111
112   assert (c != NULL);
113   rsq = c->ssm / c->sst;
114   adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
115   std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
116   t = tab_create (n_cols, n_rows, 0);
117   tab_dim (t, tab_natural_dimensions);
118   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
119   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
120   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
121   tab_vline (t, TAL_0, 1, 0, 0);
122
123   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
124   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
125   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
126   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
127   tab_float (t, 1, 1, TAB_RIGHT, sqrt (rsq), 10, 2);
128   tab_float (t, 2, 1, TAB_RIGHT, rsq, 10, 2);
129   tab_float (t, 3, 1, TAB_RIGHT, adjrsq, 10, 2);
130   tab_float (t, 4, 1, TAB_RIGHT, std_error, 10, 2);
131   tab_title (t, 0, _("Model Summary"));
132   tab_submit (t);
133 }
134
135 /*
136   Table showing estimated regression coefficients.
137  */
138 static void
139 reg_stats_coeff (pspp_linreg_cache * c)
140 {
141   size_t i;
142   size_t j;
143   int n_cols = 7;
144   int n_rows;
145   double t_stat;
146   double pval;
147   double coeff;
148   double std_err;
149   double beta;
150   const char *label;
151   struct tab_table *t;
152
153   assert (c != NULL);
154   n_rows = 2;
155
156   t = tab_create (n_cols, n_rows, 0);
157   tab_headers (t, 2, 0, 1, 0);
158   tab_dim (t, tab_natural_dimensions);
159   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
160   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
161   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
162   tab_vline (t, TAL_0, 1, 0, 0);
163
164   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
165   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
166   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
167   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
168   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
169   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
170   coeff = gsl_vector_get (c->param_estimates, 0);
171   tab_float (t, 2, 1, 0, coeff, 10, 2);
172   std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
173   tab_float (t, 3, 1, 0, std_err, 10, 2);
174   beta = coeff / c->depvar_std;
175   tab_float (t, 4, 1, 0, beta, 10, 2);
176   t_stat = coeff / std_err;
177   tab_float (t, 5, 1, 0, t_stat, 10, 2);
178   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
179   tab_float (t, 6, 1, 0, pval, 10, 2);
180   for (j = 0; j < c->n_indeps; j++)
181     {
182       i = indep_vars[j];
183       struct variable *v = cmd.v_variables[i];
184       label = var_to_string (v);
185       tab_text (t, 1, j + 2, TAB_CENTER, label);
186       /*
187          Regression coefficients.
188        */
189       coeff = gsl_vector_get (c->param_estimates, j + 1);
190       tab_float (t, 2, j + 2, 0, coeff, 10, 2);
191       /*
192          Standard error of the coefficients.
193        */
194       std_err = sqrt (gsl_matrix_get (c->cov, j + 1, j + 1));
195       tab_float (t, 3, j + 2, 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 + 1);
201       beta *= coeff / c->depvar_std;
202       tab_float (t, 4, j + 2, 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 + 2, 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 + 2, 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   int n_indep;
504   int j = 0;
505   const union value *val;
506   struct casereader *r;
507   struct casereader *r2;
508   struct ccase c;
509   const struct variable *v;
510   struct recoded_categorical_array *ca;
511   struct recoded_categorical *rc;
512   struct design_matrix *X;
513   gsl_vector *Y;
514   pspp_linreg_cache *lcache;
515   pspp_linreg_opts lopts;
516
517   n_data = casefile_get_case_cnt (cf);
518   n_indep = cmd.n_variables - cmd.n_dependent;
519   indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
520
521   Y = gsl_vector_alloc (n_data);
522   lopts.get_depvar_mean_std = 1;
523   lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
524
525   lcache = pspp_linreg_cache_alloc (n_data, n_indep);
526   lcache->indep_means = gsl_vector_alloc (n_indep);
527   lcache->indep_std = gsl_vector_alloc (n_indep);
528
529   /*
530      Read from the active file. The first pass encodes categorical
531      variables.
532    */
533   ca = cr_recoded_cat_ar_create (cmd.n_variables, cmd.v_variables);
534   for (r = casefile_get_reader (cf);
535        casereader_read (r, &c); case_destroy (&c))
536     {
537       for (i = 0; i < ca->n_vars; i++)
538         {
539           v = (*(ca->a + i))->v;
540           val = case_data (&c, v->fv);
541           cr_value_update (*(ca->a + i), val);
542         }
543     }
544   cr_create_value_matrices (ca);
545   X =
546     design_matrix_create (n_indep, (const struct variable **) cmd.v_variables,
547                           ca, n_data);
548
549   /*
550      The second pass creates the design matrix.
551    */
552   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
553        case_destroy (&c))
554     /* Iterate over the cases. */
555     {
556       k = 0;
557       row = casereader_cnum (r2) - 1;
558       for (i = 0; i < cmd.n_variables; ++i)     /* Iterate over the variables
559                                                    for the current case. 
560                                                  */
561         {
562           v = cmd.v_variables[i];
563           val = case_data (&c, v->fv);
564           /*
565              Independent/dependent variable separation. The
566              'variables' subcommand specifies a varlist which contains
567              both dependent and independent variables. The dependent
568              variables are specified with the 'dependent'
569              subcommand. We need to separate the two.
570            */
571           if (is_depvar (i))
572             {
573               if (v->type != NUMERIC)
574                 {
575                   msg (SE, gettext ("Dependent variable must be numeric."));
576                   pspp_reg_rc = CMD_FAILURE;
577                   return;
578                 }
579               lcache->depvar = (const struct var *) v;
580               gsl_vector_set (Y, row, val->f);
581             }
582           else
583             {
584               if (v->type == ALPHA)
585                 {
586                   rc = cr_var_to_recoded_categorical (v, ca);
587                   design_matrix_set_categorical (X, row, v, val, rc);
588                 }
589               else if (v->type == NUMERIC)
590                 {
591                   design_matrix_set_numeric (X, row, v, val);
592                 }
593
594               indep_vars[k] = i;
595               k++;
596               lopts.get_indep_mean_std[i] = 1;
597             }
598         }
599     }
600   /*
601      Now that we know the number of coefficients, allocate space
602      and store pointers to the variables that correspond to the
603      coefficients.
604    */
605   lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
606   for (i = 0; i < X->m->size2; i++)
607     {
608       j = i + 1;                /* The first coeff is the intercept. */
609       lcache->coeff[j].v =
610         (const struct variable *) design_matrix_col_to_var (X, i);
611     }
612   /* 
613      Find the least-squares estimates and other statistics.
614    */
615   pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
616   subcommand_statistics (cmd.a_statistics, lcache);
617   gsl_vector_free (Y);
618   design_matrix_destroy (X);
619   pspp_linreg_cache_free (lcache);
620   free (lopts.get_indep_mean_std);
621   free (indep_vars);
622   casereader_destroy (r);
623   return;
624 }
625
626 /*
627   Local Variables:   
628   mode: c
629   End:
630 */