session: Fix two memory leaks.
[pspp] / src / language / stats / glm.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2010, 2011, 2012 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
16
17 #include <config.h>
18
19 #include <gsl/gsl_cdf.h>
20 #include <gsl/gsl_matrix.h>
21 #include <gsl/gsl_combination.h>
22 #include <math.h>
23
24 #include "data/case.h"
25 #include "data/casegrouper.h"
26 #include "data/casereader.h"
27 #include "data/dataset.h"
28 #include "data/dictionary.h"
29 #include "data/format.h"
30 #include "data/value.h"
31 #include "language/command.h"
32 #include "language/dictionary/split-file.h"
33 #include "language/lexer/lexer.h"
34 #include "language/lexer/value-parser.h"
35 #include "language/lexer/variable-parser.h"
36 #include "libpspp/assertion.h"
37 #include "libpspp/ll.h"
38 #include "libpspp/message.h"
39 #include "libpspp/misc.h"
40 #include "libpspp/taint.h"
41 #include "linreg/sweep.h"
42 #include "math/categoricals.h"
43 #include "math/covariance.h"
44 #include "math/interaction.h"
45 #include "math/moments.h"
46 #include "output/tab.h"
47
48 #include "gettext.h"
49 #define _(msgid) gettext (msgid)
50
51 struct glm_spec
52 {
53   size_t n_dep_vars;
54   const struct variable **dep_vars;
55
56   size_t n_factor_vars;
57   const struct variable **factor_vars;
58
59   size_t n_interactions;
60   struct interaction **interactions;
61
62   enum mv_class exclude;
63
64   /* The weight variable */
65   const struct variable *wv;
66
67   const struct dictionary *dict;
68
69   int ss_type;
70   bool intercept;
71
72   double alpha;
73
74   bool dump_coding;
75 };
76
77 struct glm_workspace
78 {
79   double total_ssq;
80   struct moments *totals;
81
82   struct categoricals *cats;
83
84   /* 
85      Sums of squares due to different variables. Element 0 is the SSE
86      for the entire model. For i > 0, element i is the SS due to
87      variable i.
88    */
89   gsl_vector *ssq;
90 };
91
92
93 /* Default design: all possible interactions */
94 static void
95 design_full (struct glm_spec *glm)
96 {
97   int sz;
98   int i = 0;
99   glm->n_interactions = (1 << glm->n_factor_vars) - 1;
100
101   glm->interactions = xcalloc (glm->n_interactions, sizeof *glm->interactions);
102
103   /* All subsets, with exception of the empty set, of [0, glm->n_factor_vars) */
104   for (sz = 1; sz <= glm->n_factor_vars; ++sz)
105     {
106       gsl_combination *c = gsl_combination_calloc (glm->n_factor_vars, sz);
107
108       do
109         {
110           struct interaction *iact = interaction_create (NULL);
111           int e;
112           for (e = 0 ; e < gsl_combination_k (c); ++e)
113             interaction_add_variable (iact, glm->factor_vars [gsl_combination_get (c, e)]);
114
115           glm->interactions[i++] = iact;
116         }
117       while (gsl_combination_next (c) == GSL_SUCCESS);
118
119       gsl_combination_free (c);
120     }
121 }
122
123 static void output_glm (const struct glm_spec *,
124                         const struct glm_workspace *ws);
125 static void run_glm (struct glm_spec *cmd, struct casereader *input,
126                      const struct dataset *ds);
127
128
129 static bool parse_design_spec (struct lexer *lexer, struct glm_spec *glm);
130
131
132 int
133 cmd_glm (struct lexer *lexer, struct dataset *ds)
134 {
135   int i;
136   struct const_var_set *factors = NULL;
137   struct glm_spec glm;
138   bool design = false;
139   glm.dict = dataset_dict (ds);
140   glm.n_dep_vars = 0;
141   glm.n_factor_vars = 0;
142   glm.n_interactions = 0;
143   glm.interactions = NULL;
144   glm.dep_vars = NULL;
145   glm.factor_vars = NULL;
146   glm.exclude = MV_ANY;
147   glm.intercept = true;
148   glm.wv = dict_get_weight (glm.dict);
149   glm.alpha = 0.05;
150   glm.dump_coding = false;
151   glm.ss_type = 3;
152
153   if (!parse_variables_const (lexer, glm.dict,
154                               &glm.dep_vars, &glm.n_dep_vars,
155                               PV_NO_DUPLICATE | PV_NUMERIC))
156     goto error;
157
158   lex_force_match (lexer, T_BY);
159
160   if (!parse_variables_const (lexer, glm.dict,
161                               &glm.factor_vars, &glm.n_factor_vars,
162                               PV_NO_DUPLICATE | PV_NUMERIC))
163     goto error;
164
165   if (glm.n_dep_vars > 1)
166     {
167       msg (ME, _("Multivariate analysis is not yet implemented"));
168       return CMD_FAILURE;
169     }
170
171   factors =
172     const_var_set_create_from_array (glm.factor_vars, glm.n_factor_vars);
173
174   while (lex_token (lexer) != T_ENDCMD)
175     {
176       lex_match (lexer, T_SLASH);
177
178       if (lex_match_id (lexer, "MISSING"))
179         {
180           lex_match (lexer, T_EQUALS);
181           while (lex_token (lexer) != T_ENDCMD
182                  && lex_token (lexer) != T_SLASH)
183             {
184               if (lex_match_id (lexer, "INCLUDE"))
185                 {
186                   glm.exclude = MV_SYSTEM;
187                 }
188               else if (lex_match_id (lexer, "EXCLUDE"))
189                 {
190                   glm.exclude = MV_ANY;
191                 }
192               else
193                 {
194                   lex_error (lexer, NULL);
195                   goto error;
196                 }
197             }
198         }
199       else if (lex_match_id (lexer, "INTERCEPT"))
200         {
201           lex_match (lexer, T_EQUALS);
202           while (lex_token (lexer) != T_ENDCMD
203                  && lex_token (lexer) != T_SLASH)
204             {
205               if (lex_match_id (lexer, "INCLUDE"))
206                 {
207                   glm.intercept = true;
208                 }
209               else if (lex_match_id (lexer, "EXCLUDE"))
210                 {
211                   glm.intercept = false;
212                 }
213               else
214                 {
215                   lex_error (lexer, NULL);
216                   goto error;
217                 }
218             }
219         }
220       else if (lex_match_id (lexer, "CRITERIA"))
221         {
222           lex_match (lexer, T_EQUALS);
223           if (lex_match_id (lexer, "ALPHA"))
224             {
225               if (lex_force_match (lexer, T_LPAREN))
226                 {
227                   if (! lex_force_num (lexer))
228                     {
229                       lex_error (lexer, NULL);
230                       goto error;
231                     }
232                   
233                   glm.alpha = lex_number (lexer);
234                   lex_get (lexer);
235                   if ( ! lex_force_match (lexer, T_RPAREN))
236                     {
237                       lex_error (lexer, NULL);
238                       goto error;
239                     }
240                 }
241             }
242           else
243             {
244               lex_error (lexer, NULL);
245               goto error;
246             }
247         }
248       else if (lex_match_id (lexer, "METHOD"))
249         {
250           lex_match (lexer, T_EQUALS);
251           if ( !lex_force_match_id (lexer, "SSTYPE"))
252             {
253               lex_error (lexer, NULL);
254               goto error;
255             }
256
257           if ( ! lex_force_match (lexer, T_LPAREN))
258             {
259               lex_error (lexer, NULL);
260               goto error;
261             }
262
263           if ( ! lex_force_int (lexer))
264             {
265               lex_error (lexer, NULL);
266               goto error;
267             }
268
269           glm.ss_type = lex_integer (lexer);
270           if (1 > glm.ss_type  && 3 < glm.ss_type )
271             {
272               msg (ME, _("Only types 1, 2 & 3 sums of squares are currently implemented"));
273               goto error;
274             }
275
276           lex_get (lexer);
277
278           if ( ! lex_force_match (lexer, T_RPAREN))
279             {
280               lex_error (lexer, NULL);
281               goto error;
282             }
283         }
284       else if (lex_match_id (lexer, "DESIGN"))
285         {
286           lex_match (lexer, T_EQUALS);
287
288           if (! parse_design_spec (lexer, &glm))
289             goto error;
290
291           if (glm.n_interactions > 0)
292             design = true;
293         }
294       else if (lex_match_id (lexer, "SHOWCODES"))
295         /* Undocumented debug option */
296         {
297           lex_match (lexer, T_EQUALS);
298
299           glm.dump_coding = true;
300         }
301       else
302         {
303           lex_error (lexer, NULL);
304           goto error;
305         }
306     }
307
308   if ( ! design )
309     {
310       design_full (&glm);
311     }
312
313   {
314     struct casegrouper *grouper;
315     struct casereader *group;
316     bool ok;
317
318     grouper = casegrouper_create_splits (proc_open (ds), glm.dict);
319     while (casegrouper_get_next_group (grouper, &group))
320       run_glm (&glm, group, ds);
321     ok = casegrouper_destroy (grouper);
322     ok = proc_commit (ds) && ok;
323   }
324
325   const_var_set_destroy (factors);
326   free (glm.factor_vars);
327   for (i = 0 ; i < glm.n_interactions; ++i)
328     interaction_destroy (glm.interactions[i]);
329
330   free (glm.interactions);
331   free (glm.dep_vars);
332
333
334   return CMD_SUCCESS;
335
336 error:
337
338   const_var_set_destroy (factors);
339   free (glm.factor_vars);
340   for (i = 0 ; i < glm.n_interactions; ++i)
341     interaction_destroy (glm.interactions[i]);
342
343   free (glm.interactions);
344   free (glm.dep_vars);
345
346   return CMD_FAILURE;
347 }
348
349 static inline bool
350 not_dropped (size_t j, const bool *ff)
351 {
352   return ! ff[j];
353 }
354
355 static void
356 fill_submatrix (const gsl_matrix * cov, gsl_matrix * submatrix, bool *dropped_f)
357 {
358   size_t i;
359   size_t j;
360   size_t n = 0;
361   size_t m = 0;
362   
363   for (i = 0; i < cov->size1; i++)
364     {
365       if (not_dropped (i, dropped_f))
366         {         
367           m = 0;
368           for (j = 0; j < cov->size2; j++)
369             {
370               if (not_dropped (j, dropped_f))
371                 {
372                   gsl_matrix_set (submatrix, n, m,
373                                   gsl_matrix_get (cov, i, j));
374                   m++;
375                 }       
376             }
377           n++;
378         }
379     }
380 }
381
382
383 /* 
384    Type 1 sums of squares.
385    Populate SSQ with the Type 1 sums of squares according to COV
386  */
387 static void
388 ssq_type1 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
389 {
390   gsl_matrix *cm = covariance_calculate_unnormalized (cov);
391   size_t i;
392   size_t k;
393   bool *model_dropped = xcalloc (covariance_dim (cov), sizeof (*model_dropped));
394   bool *submodel_dropped = xcalloc (covariance_dim (cov), sizeof (*submodel_dropped));
395   const struct categoricals *cats = covariance_get_categoricals (cov);
396
397   size_t n_dropped_model = 0;
398   size_t n_dropped_submodel = 0;
399
400   for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
401     {
402       n_dropped_model++;
403       n_dropped_submodel++;
404       model_dropped[i] = true;
405       submodel_dropped[i] = true;
406     }
407
408   for (k = 0; k < cmd->n_interactions; k++)
409     {
410       gsl_matrix *model_cov = NULL;
411       gsl_matrix *submodel_cov = NULL;
412       
413       n_dropped_submodel = n_dropped_model;
414       for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
415         {
416           submodel_dropped[i] = model_dropped[i];
417         }
418
419       for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
420         {
421           const struct interaction * x = 
422             categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars);
423
424           if ( x == cmd->interactions [k])
425             {
426               model_dropped[i] = false;
427               n_dropped_model--;
428             }
429         }
430
431       model_cov = gsl_matrix_alloc (cm->size1 - n_dropped_model, cm->size2 - n_dropped_model);
432       submodel_cov = gsl_matrix_alloc (cm->size1 - n_dropped_submodel, cm->size2 - n_dropped_submodel);
433
434       fill_submatrix (cm, model_cov,    model_dropped);
435       fill_submatrix (cm, submodel_cov, submodel_dropped);
436
437       reg_sweep (model_cov, 0);
438       reg_sweep (submodel_cov, 0);
439
440       gsl_vector_set (ssq, k + 1,
441                       gsl_matrix_get (submodel_cov, 0, 0) - gsl_matrix_get (model_cov, 0, 0)
442                       );
443
444       gsl_matrix_free (model_cov);
445       gsl_matrix_free (submodel_cov);
446     }
447
448   free (model_dropped);
449   free (submodel_dropped);
450   gsl_matrix_free (cm);
451 }
452
453 /* 
454    Type 2 sums of squares.
455    Populate SSQ with the Type 2 sums of squares according to COV
456  */
457 static void
458 ssq_type2 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
459 {
460   gsl_matrix *cm = covariance_calculate_unnormalized (cov);
461   size_t i;
462   size_t k;
463   bool *model_dropped = xcalloc (covariance_dim (cov), sizeof (*model_dropped));
464   bool *submodel_dropped = xcalloc (covariance_dim (cov), sizeof (*submodel_dropped));
465   const struct categoricals *cats = covariance_get_categoricals (cov);
466
467   for (k = 0; k < cmd->n_interactions; k++)
468     {
469       gsl_matrix *model_cov = NULL;
470       gsl_matrix *submodel_cov = NULL;
471       size_t n_dropped_model = 0;
472       size_t n_dropped_submodel = 0;
473       for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
474         {
475           const struct interaction * x = 
476             categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars);
477
478           model_dropped[i] = false;
479           submodel_dropped[i] = false;
480           if (interaction_is_subset (cmd->interactions [k], x))
481             {
482               assert (n_dropped_submodel < covariance_dim (cov));
483               n_dropped_submodel++;
484               submodel_dropped[i] = true;
485
486               if ( cmd->interactions [k]->n_vars < x->n_vars)
487                 {
488                   assert (n_dropped_model < covariance_dim (cov));
489                   n_dropped_model++;
490                   model_dropped[i] = true;
491                 }
492             }
493         }
494
495       model_cov = gsl_matrix_alloc (cm->size1 - n_dropped_model, cm->size2 - n_dropped_model);
496       submodel_cov = gsl_matrix_alloc (cm->size1 - n_dropped_submodel, cm->size2 - n_dropped_submodel);
497
498       fill_submatrix (cm, model_cov,    model_dropped);
499       fill_submatrix (cm, submodel_cov, submodel_dropped);
500
501       reg_sweep (model_cov, 0);
502       reg_sweep (submodel_cov, 0);
503
504       gsl_vector_set (ssq, k + 1,
505                       gsl_matrix_get (submodel_cov, 0, 0) - gsl_matrix_get (model_cov, 0, 0)
506                       );
507
508       gsl_matrix_free (model_cov);
509       gsl_matrix_free (submodel_cov);
510     }
511
512   free (model_dropped);
513   free (submodel_dropped);
514   gsl_matrix_free (cm);
515 }
516
517 /* 
518    Type 3 sums of squares.
519    Populate SSQ with the Type 2 sums of squares according to COV
520  */
521 static void
522 ssq_type3 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
523 {
524   gsl_matrix *cm = covariance_calculate_unnormalized (cov);
525   size_t i;
526   size_t k;
527   bool *model_dropped = xcalloc (covariance_dim (cov), sizeof (*model_dropped));
528   bool *submodel_dropped = xcalloc (covariance_dim (cov), sizeof (*submodel_dropped));
529   const struct categoricals *cats = covariance_get_categoricals (cov);
530
531   double ss0;
532   gsl_matrix *submodel_cov = gsl_matrix_alloc (cm->size1, cm->size2);
533   fill_submatrix (cm, submodel_cov, submodel_dropped);
534   reg_sweep (submodel_cov, 0);
535   ss0 = gsl_matrix_get (submodel_cov, 0, 0);
536   gsl_matrix_free (submodel_cov);
537   free (submodel_dropped);
538
539   for (k = 0; k < cmd->n_interactions; k++)
540     {
541       gsl_matrix *model_cov = NULL;
542       size_t n_dropped_model = 0;
543
544       for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
545         {
546           const struct interaction * x = 
547             categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars);
548
549           model_dropped[i] = false;
550
551           if ( cmd->interactions [k] == x)
552             {
553               assert (n_dropped_model < covariance_dim (cov));
554               n_dropped_model++;
555               model_dropped[i] = true;
556             }
557         }
558
559       model_cov = gsl_matrix_alloc (cm->size1 - n_dropped_model, cm->size2 - n_dropped_model);
560
561       fill_submatrix (cm, model_cov,    model_dropped);
562
563       reg_sweep (model_cov, 0);
564
565       gsl_vector_set (ssq, k + 1,
566                       gsl_matrix_get (model_cov, 0, 0) - ss0);
567
568       gsl_matrix_free (model_cov);
569     }
570   free (model_dropped);
571
572   gsl_matrix_free (cm);
573 }
574
575
576
577 //static  void dump_matrix (const gsl_matrix *m);
578
579 static void
580 run_glm (struct glm_spec *cmd, struct casereader *input,
581          const struct dataset *ds)
582 {
583   bool warn_bad_weight = true;
584   int v;
585   struct taint *taint;
586   struct dictionary *dict = dataset_dict (ds);
587   struct casereader *reader;
588   struct ccase *c;
589
590   struct glm_workspace ws;
591   struct covariance *cov;
592
593   ws.cats = categoricals_create (cmd->interactions, cmd->n_interactions,
594                                  cmd->wv, cmd->exclude, MV_ANY);
595
596   cov = covariance_2pass_create (cmd->n_dep_vars, cmd->dep_vars,
597                                  ws.cats, cmd->wv, cmd->exclude);
598
599
600   c = casereader_peek (input, 0);
601   if (c == NULL)
602     {
603       casereader_destroy (input);
604       return;
605     }
606   output_split_file_values (ds, c);
607   case_unref (c);
608
609   taint = taint_clone (casereader_get_taint (input));
610
611   ws.totals = moments_create (MOMENT_VARIANCE);
612
613   for (reader = casereader_clone (input);
614        (c = casereader_read (reader)) != NULL; case_unref (c))
615     {
616       double weight = dict_get_case_weight (dict, c, &warn_bad_weight);
617
618       for (v = 0; v < cmd->n_dep_vars; ++v)
619         moments_pass_one (ws.totals, case_data (c, cmd->dep_vars[v])->f,
620                           weight);
621
622       covariance_accumulate_pass1 (cov, c);
623     }
624   casereader_destroy (reader);
625
626   if (cmd->dump_coding)
627     reader = casereader_clone (input);
628   else
629     reader = input;
630
631   for (;
632        (c = casereader_read (reader)) != NULL; case_unref (c))
633     {
634       double weight = dict_get_case_weight (dict, c, &warn_bad_weight);
635
636       for (v = 0; v < cmd->n_dep_vars; ++v)
637         moments_pass_two (ws.totals, case_data (c, cmd->dep_vars[v])->f,
638                           weight);
639
640       covariance_accumulate_pass2 (cov, c);
641     }
642   casereader_destroy (reader);
643
644
645   if (cmd->dump_coding)
646     {
647       struct tab_table *t =
648         covariance_dump_enc_header (cov,
649                                     1 + casereader_count_cases (input));
650       for (reader = input;
651            (c = casereader_read (reader)) != NULL; case_unref (c))
652         {
653           covariance_dump_enc (cov, c, t);
654         }
655       casereader_destroy (reader);
656       tab_submit (t);
657     }
658
659   {
660     gsl_matrix *cm = covariance_calculate_unnormalized (cov);
661
662     //    dump_matrix (cm);
663
664     ws.total_ssq = gsl_matrix_get (cm, 0, 0);
665
666     reg_sweep (cm, 0);
667
668     /*
669       Store the overall SSE.
670     */
671     ws.ssq = gsl_vector_alloc (cm->size1);
672     gsl_vector_set (ws.ssq, 0, gsl_matrix_get (cm, 0, 0));
673     switch (cmd->ss_type)
674       {
675       case 1:
676         ssq_type1 (cov, ws.ssq, cmd);
677         break;
678       case 2:
679         ssq_type2 (cov, ws.ssq, cmd);
680         break;
681       case 3:
682         ssq_type3 (cov, ws.ssq, cmd);
683         break;
684       default:
685         NOT_REACHED ();
686         break;
687       }
688     //    dump_matrix (cm);
689
690     gsl_matrix_free (cm);
691   }
692
693   if (!taint_has_tainted_successor (taint))
694     output_glm (cmd, &ws);
695
696   gsl_vector_free (ws.ssq);
697
698   covariance_destroy (cov);
699   moments_destroy (ws.totals);
700
701   taint_destroy (taint);
702 }
703
704 static const char *roman[] = 
705   {
706     "", /* The Romans had no concept of zero */
707     "I",
708     "II",
709     "III",
710     "IV"
711   };
712
713 static void
714 output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
715 {
716   const struct fmt_spec *wfmt =
717     cmd->wv ? var_get_print_format (cmd->wv) : &F_8_0;
718
719   double intercept_ssq;
720   double ssq_effects;
721   double n_total, mean;
722   double df_corr = 1.0;
723   double mse = 0;
724
725   int f;
726   int r;
727   const int heading_columns = 1;
728   const int heading_rows = 1;
729   struct tab_table *t;
730
731   const int nc = 6;
732   int nr = heading_rows + 3 + cmd->n_interactions;
733   if (cmd->intercept)
734     nr += 2;
735
736   msg (MW, "GLM is experimental.  Do not rely on these results.");
737   t = tab_create (nc, nr);
738   tab_title (t, _("Tests of Between-Subjects Effects"));
739
740   tab_headers (t, heading_columns, 0, heading_rows, 0);
741
742   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
743
744   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
745   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
746
747   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
748
749   /* TRANSLATORS: The parameter is a roman numeral */
750   tab_text_format (t, 1, 0, TAB_CENTER | TAT_TITLE,
751                    _("Type %s Sum of Squares"), 
752                    roman[cmd->ss_type]);
753   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df"));
754   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
755   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("F"));
756   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
757
758   moments_calculate (ws->totals, &n_total, &mean, NULL, NULL, NULL);
759
760   df_corr += categoricals_df_total (ws->cats);
761
762   r = heading_rows;
763   if (cmd->intercept)
764     tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Model"));
765   else
766     tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Model"));
767
768   r++;
769
770   mse = gsl_vector_get (ws->ssq, 0) / (n_total - df_corr);
771
772   intercept_ssq = pow2 (mean * n_total) / n_total;
773
774   ssq_effects = 0.0;
775   if (cmd->intercept)
776     {
777       const double df = 1.0;
778       const double F = intercept_ssq / df / mse;
779       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Intercept"));
780       tab_double (t, 1, r, 0, intercept_ssq, NULL);
781       tab_double (t, 2, r, 0, 1.00, wfmt);
782       tab_double (t, 3, r, 0, intercept_ssq / df, NULL);
783       tab_double (t, 4, r, 0, F, NULL);
784       tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr),
785                   NULL);
786       r++;
787     }
788
789   for (f = 0; f < cmd->n_interactions; ++f)
790     {
791       struct string str = DS_EMPTY_INITIALIZER;
792       double df = categoricals_df (ws->cats, f);
793
794       double ssq = gsl_vector_get (ws->ssq, f + 1);
795       double F;
796
797       ssq_effects += ssq;
798
799       if (! cmd->intercept) 
800         {
801           df++;
802           ssq += intercept_ssq;
803         }
804
805       F = ssq / df / mse;
806       interaction_to_string (cmd->interactions[f], &str);
807       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, ds_cstr (&str));
808       ds_destroy (&str);
809
810       tab_double (t, 1, r, 0, ssq, NULL);
811       tab_double (t, 2, r, 0, df, wfmt);
812       tab_double (t, 3, r, 0, ssq / df, NULL);
813       tab_double (t, 4, r, 0, F, NULL);
814
815       tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr),
816                   NULL);
817       r++;
818     }
819
820   {
821     /* Model / Corrected Model */
822     double df = df_corr;
823     double ssq = ws->total_ssq - gsl_vector_get (ws->ssq, 0);
824     double F;
825
826     if ( cmd->intercept )
827       df --;
828     else
829       ssq += intercept_ssq;
830
831     F = ssq / df / mse;
832     tab_double (t, 1, heading_rows, 0, ssq, NULL);
833     tab_double (t, 2, heading_rows, 0, df, wfmt);
834     tab_double (t, 3, heading_rows, 0, ssq / df, NULL);
835     tab_double (t, 4, heading_rows, 0, F, NULL);
836
837     tab_double (t, 5, heading_rows, 0,
838                 gsl_cdf_fdist_Q (F, df, n_total - df_corr), NULL);
839   }
840
841   {
842     const double df = n_total - df_corr;
843     const double ssq = gsl_vector_get (ws->ssq, 0);
844     const double mse = ssq / df;
845     tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Error"));
846     tab_double (t, 1, r, 0, ssq, NULL);
847     tab_double (t, 2, r, 0, df, wfmt);
848     tab_double (t, 3, r++, 0, mse, NULL);
849   }
850
851   {
852     tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Total"));
853     tab_double (t, 1, r, 0, ws->total_ssq + intercept_ssq, NULL);
854     tab_double (t, 2, r, 0, n_total, wfmt);
855     
856     r++;
857   }
858
859   if (cmd->intercept)
860     {
861       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Total"));
862       tab_double (t, 1, r, 0, ws->total_ssq, NULL);
863       tab_double (t, 2, r, 0, n_total - 1.0, wfmt);
864     }
865
866   tab_submit (t);
867 }
868
869 #if 0
870 static void
871 dump_matrix (const gsl_matrix * m)
872 {
873   size_t i, j;
874   for (i = 0; i < m->size1; ++i)
875     {
876       for (j = 0; j < m->size2; ++j)
877         {
878           double x = gsl_matrix_get (m, i, j);
879           printf ("%.3f ", x);
880         }
881       printf ("\n");
882     }
883   printf ("\n");
884 }
885 #endif
886
887
888 \f
889 static bool
890 parse_nested_variable (struct lexer *lexer, struct glm_spec *glm)
891 {
892   const struct variable *v = NULL;
893   if ( ! lex_match_variable (lexer, glm->dict, &v))
894     return false;
895
896   if (lex_match (lexer, T_LPAREN))
897     {
898       if ( ! parse_nested_variable (lexer, glm))
899         return false;
900
901       if ( ! lex_force_match (lexer, T_RPAREN))
902         return false;
903     }
904
905   lex_error (lexer, "Nested variables are not yet implemented"); return false;  
906   return true;
907 }
908
909 /* A design term is an interaction OR a nested variable */
910 static bool
911 parse_design_term (struct lexer *lexer, struct glm_spec *glm)
912 {
913   struct interaction *iact = NULL;
914   if (parse_design_interaction (lexer, glm->dict, &iact))
915     {
916       /* Interaction parsing successful.  Add to list of interactions */
917       glm->interactions = xrealloc (glm->interactions, sizeof *glm->interactions * ++glm->n_interactions);
918       glm->interactions[glm->n_interactions - 1] = iact;
919       return true;
920     }
921
922   if ( parse_nested_variable (lexer, glm))
923     return true;
924
925   return false;
926 }
927
928
929
930 /* Parse a complete DESIGN specification.
931    A design spec is a design term, optionally followed by a comma,
932    and another design spec.
933 */
934 static bool
935 parse_design_spec (struct lexer *lexer, struct glm_spec *glm)
936 {
937   if  (lex_token (lexer) == T_ENDCMD || lex_token (lexer) == T_SLASH)
938     return true;
939
940   if ( ! parse_design_term (lexer, glm))
941     return false;
942
943   lex_match (lexer, T_COMMA);
944
945   return parse_design_spec (lexer, glm);
946 }
947