Add aux data parameter to q2c parse_<command>() and custom parser
[pspp-builds.git] / src / language / stats / oneway.q
1 /* PSPP - One way ANOVA. -*-c-*-
2
3 Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
4 Author: John Darrington 2004
5
6 This program is free software; you can redistribute it and/or
7 modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 2 of the
9 License, or (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful, but
12 WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
19 02110-1301, USA. */
20
21 #include <config.h>
22
23 #include <gsl/gsl_cdf.h>
24 #include <math.h>
25 #include <stdio.h>
26 #include <stdlib.h>
27
28 #include <data/case.h>
29 #include <data/casefile.h>
30 #include <data/dictionary.h>
31 #include <data/procedure.h>
32 #include <data/value-labels.h>
33 #include <data/variable.h>
34 #include <language/command.h>
35 #include <language/dictionary/split-file.h>
36 #include <language/lexer/lexer.h>
37 #include <libpspp/alloc.h>
38 #include <libpspp/compiler.h>
39 #include <libpspp/hash.h>
40 #include <libpspp/magic.h>
41 #include <libpspp/message.h>
42 #include <libpspp/message.h>
43 #include <libpspp/misc.h>
44 #include <libpspp/str.h>
45 #include <math/group-proc.h>
46 #include <math/group.h>
47 #include <math/levene.h>
48 #include <output/manager.h>
49 #include <output/table.h>
50 #include "sort-criteria.h"
51
52 #include "gettext.h"
53 #define _(msgid) gettext (msgid)
54
55 /* (headers) */
56
57 /* (specification)
58    "ONEWAY" (oneway_):
59    *^variables=custom;
60    +missing=miss:!analysis/listwise,
61    incl:include/!exclude;
62    contrast= double list;
63    statistics[st_]=descriptives,homogeneity.
64 */
65 /* (declarations) */
66 /* (functions) */
67
68
69
70 static int bad_weight_warn = 1;
71
72
73 static struct cmd_oneway cmd;
74
75 /* The independent variable */
76 static struct variable *indep_var;
77
78 /* Number of dependent variables */
79 static size_t n_vars;
80
81 /* The dependent variables */
82 static struct variable **vars;
83
84
85 /* A  hash table containing all the distinct values of the independent
86    variables */
87 static struct hsh_table *global_group_hash ;
88
89 /* The number of distinct values of the independent variable, when all 
90    missing values are disregarded */
91 static int ostensible_number_of_groups=-1;
92
93
94 /* Function to use for testing for missing values */
95 static is_missing_func *value_is_missing;
96
97
98 static bool run_oneway(const struct ccase *first,
99                        const struct casefile *cf, void *_mode);
100
101
102 /* Routines to show the output tables */
103 static void show_anova_table(void);
104 static void show_descriptives(void);
105 static void show_homogeneity(void);
106
107 static void show_contrast_coeffs(short *);
108 static void show_contrast_tests(short *);
109
110
111 enum stat_table_t {STAT_DESC = 1, STAT_HOMO = 2};
112
113 static enum stat_table_t stat_tables ;
114
115 void output_oneway(void);
116
117
118 int
119 cmd_oneway(void)
120 {
121   int i;
122   bool ok;
123
124   if ( !parse_oneway(&cmd, NULL) )
125     return CMD_FAILURE;
126
127   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
128   if (cmd.incl == ONEWAY_INCLUDE ) 
129     value_is_missing = mv_is_value_system_missing;
130   else
131     value_is_missing = mv_is_value_missing;
132
133   /* What statistics were requested */
134   if ( cmd.sbc_statistics ) 
135     {
136
137       for (i = 0 ; i < ONEWAY_ST_count ; ++i ) 
138         {
139           if  ( ! cmd.a_statistics[i]  ) continue;
140
141           switch (i) {
142           case ONEWAY_ST_DESCRIPTIVES:
143             stat_tables |= STAT_DESC;
144             break;
145           case ONEWAY_ST_HOMOGENEITY:
146             stat_tables |= STAT_HOMO;
147             break;
148           }
149         }
150     }
151
152   ok = multipass_procedure_with_splits (run_oneway, &cmd);
153
154   free (vars);
155   free_oneway (&cmd);
156
157   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
158 }
159
160
161 void
162 output_oneway(void)
163 {
164   size_t i;
165   short *bad_contrast ; 
166
167   bad_contrast = xnmalloc (cmd.sbc_contrast, sizeof *bad_contrast);
168
169   /* Check the sanity of the given contrast values */
170   for (i = 0 ; i < cmd.sbc_contrast ; ++i ) 
171     {
172       int j;
173       double sum = 0;
174
175       bad_contrast[i] = 0;
176       if ( subc_list_double_count(&cmd.dl_contrast[i]) != 
177            ostensible_number_of_groups )
178         {
179           msg(SW, 
180               _("Number of contrast coefficients must equal the number of groups"));
181           bad_contrast[i] = 1;
182           continue;
183         }
184
185       for (j=0; j < ostensible_number_of_groups ; ++j )
186         sum += subc_list_double_at(&cmd.dl_contrast[i],j);
187
188       if ( sum != 0.0 ) 
189         msg(SW,_("Coefficients for contrast %d do not total zero"),i + 1);
190     }
191
192   if ( stat_tables & STAT_DESC ) 
193     show_descriptives();
194
195   if ( stat_tables & STAT_HOMO )
196     show_homogeneity();
197
198   show_anova_table();
199      
200   if (cmd.sbc_contrast )
201     {
202       show_contrast_coeffs(bad_contrast);
203       show_contrast_tests(bad_contrast);
204     }
205
206
207   free(bad_contrast);
208
209   /* Clean up */
210   for (i = 0 ; i < n_vars ; ++i ) 
211     {
212       struct hsh_table *group_hash = group_proc_get (vars[i])->group_hash;
213
214       hsh_destroy(group_hash);
215     }
216
217   hsh_destroy(global_group_hash);
218
219 }
220
221
222
223
224 /* Parser for the variables sub command */
225 static int
226 oneway_custom_variables(struct cmd_oneway *cmd UNUSED, void *aux UNUSED)
227 {
228
229   lex_match('=');
230
231   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
232       && token != T_ALL)
233     return 2;
234   
235
236   if (!parse_variables (default_dict, &vars, &n_vars,
237                         PV_DUPLICATE 
238                         | PV_NUMERIC | PV_NO_SCRATCH) )
239     {
240       free (vars);
241       return 0;
242     }
243
244   assert(n_vars);
245
246   if ( ! lex_match(T_BY))
247     return 2;
248
249
250   indep_var = parse_variable();
251
252   if ( !indep_var ) 
253     {
254       msg(SE,_("`%s' is not a variable name"),tokid);
255       return 0;
256     }
257
258
259   return 1;
260 }
261
262
263 /* Show the ANOVA table */
264 static void  
265 show_anova_table(void)
266 {
267   size_t i;
268   int n_cols =7;
269   size_t n_rows = n_vars * 3 + 1;
270
271   struct tab_table *t;
272
273
274   t = tab_create (n_cols,n_rows,0);
275   tab_headers (t, 2, 0, 1, 0);
276   tab_dim (t, tab_natural_dimensions);
277
278
279   tab_box (t, 
280            TAL_2, TAL_2,
281            -1, TAL_1,
282            0, 0,
283            n_cols - 1, n_rows - 1);
284
285   tab_hline (t, TAL_2, 0, n_cols - 1, 1 );
286   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
287   tab_vline (t, TAL_0, 1, 0, 0);
288   
289   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
290   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
291   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
292   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
293   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
294
295
296   for ( i=0 ; i < n_vars ; ++i ) 
297     {
298       struct group_statistics *totals = &group_proc_get (vars[i])->ugs;
299       struct hsh_table *group_hash = group_proc_get (vars[i])->group_hash;
300       struct hsh_iterator g;
301       struct group_statistics *gs;
302       double ssa=0;
303       const char *s = var_to_string(vars[i]);
304
305       for (gs =  hsh_first (group_hash,&g); 
306            gs != 0; 
307            gs = hsh_next(group_hash,&g))
308         {
309           ssa += (gs->sum * gs->sum)/gs->n;
310         }
311       
312       ssa -= ( totals->sum * totals->sum ) / totals->n ;
313
314       tab_text (t, 0, i * 3 + 1, TAB_LEFT | TAT_TITLE, s);
315       tab_text (t, 1, i * 3 + 1, TAB_LEFT | TAT_TITLE, _("Between Groups"));
316       tab_text (t, 1, i * 3 + 2, TAB_LEFT | TAT_TITLE, _("Within Groups"));
317       tab_text (t, 1, i * 3 + 3, TAB_LEFT | TAT_TITLE, _("Total"));
318       
319       if (i > 0)
320         tab_hline(t, TAL_1, 0, n_cols - 1 , i * 3 + 1);
321
322       {
323         struct group_proc *gp = group_proc_get (vars[i]);
324         const double sst = totals->ssq - ( totals->sum * totals->sum) / totals->n ;
325         const double df1 = gp->n_groups - 1;
326         const double df2 = totals->n - gp->n_groups ;
327         const double msa = ssa / df1;
328         
329         gp->mse  = (sst - ssa) / df2;
330         
331         
332         /* Sums of Squares */
333         tab_float (t, 2, i * 3 + 1, 0, ssa, 10, 2);
334         tab_float (t, 2, i * 3 + 3, 0, sst, 10, 2);
335         tab_float (t, 2, i * 3 + 2, 0, sst - ssa, 10, 2);
336
337
338         /* Degrees of freedom */
339         tab_float (t, 3, i * 3 + 1, 0, df1, 4, 0);
340         tab_float (t, 3, i * 3 + 2, 0, df2, 4, 0);
341         tab_float (t, 3, i * 3 + 3, 0, totals->n - 1, 4, 0);
342
343         /* Mean Squares */
344         tab_float (t, 4, i * 3 + 1, TAB_RIGHT, msa, 8, 3);
345         tab_float (t, 4, i * 3 + 2, TAB_RIGHT, gp->mse, 8, 3);
346         
347
348         { 
349           const double F = msa/gp->mse ;
350
351           /* The F value */
352           tab_float (t, 5, i * 3 + 1, 0,  F, 8, 3);
353         
354           /* The significance */
355           tab_float (t, 6, i * 3 + 1, 0, gsl_cdf_fdist_Q(F,df1,df2), 8, 3);
356         }
357
358       }
359
360     }
361
362
363   tab_title (t, _("ANOVA"));
364   tab_submit (t);
365
366
367 }
368
369 /* Show the descriptives table */
370 static void  
371 show_descriptives(void)
372 {
373   size_t v;
374   int n_cols =10;
375   struct tab_table *t;
376   int row;
377
378   const double confidence=0.95;
379   const double q = (1.0 - confidence) / 2.0;
380
381   
382   int n_rows = 2 ; 
383
384   for ( v = 0 ; v < n_vars ; ++v ) 
385     n_rows += group_proc_get (vars[v])->n_groups + 1;
386
387   t = tab_create (n_cols,n_rows,0);
388   tab_headers (t, 2, 0, 2, 0);
389   tab_dim (t, tab_natural_dimensions);
390
391
392   /* Put a frame around the entire box, and vertical lines inside */
393   tab_box (t, 
394            TAL_2, TAL_2,
395            -1, TAL_1,
396            0, 0,
397            n_cols - 1, n_rows - 1);
398
399   /* Underline headers */
400   tab_hline (t, TAL_2, 0, n_cols - 1, 2 );
401   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
402   
403   tab_text (t, 2, 1, TAB_CENTER | TAT_TITLE, _("N"));
404   tab_text (t, 3, 1, TAB_CENTER | TAT_TITLE, _("Mean"));
405   tab_text (t, 4, 1, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
406   tab_text (t, 5, 1, TAB_CENTER | TAT_TITLE, _("Std. Error"));
407
408
409   tab_vline(t, TAL_0, 7, 0, 0);
410   tab_hline(t, TAL_1, 6, 7, 1);
411   tab_joint_text (t, 6, 0, 7, 0, TAB_CENTER | TAT_TITLE | TAT_PRINTF, _("%g%% Confidence Interval for Mean"),confidence*100.0);
412
413   tab_text (t, 6, 1, TAB_CENTER | TAT_TITLE, _("Lower Bound"));
414   tab_text (t, 7, 1, TAB_CENTER | TAT_TITLE, _("Upper Bound"));
415
416   tab_text (t, 8, 1, TAB_CENTER | TAT_TITLE, _("Minimum"));
417   tab_text (t, 9, 1, TAB_CENTER | TAT_TITLE, _("Maximum"));
418
419
420   tab_title (t, _("Descriptives"));
421
422
423   row = 2;
424   for ( v=0 ; v < n_vars ; ++v ) 
425     {
426       double T;
427       double std_error;
428       
429       struct group_proc *gp = group_proc_get (vars[v]);
430
431       struct group_statistics *gs;
432       struct group_statistics *totals = &gp->ugs; 
433
434       const char *s = var_to_string(vars[v]);
435
436       struct group_statistics *const *gs_array =
437         (struct group_statistics *const *) hsh_sort(gp->group_hash);
438       int count = 0;
439
440       tab_text (t, 0, row, TAB_LEFT | TAT_TITLE, s);
441       if ( v > 0) 
442         tab_hline(t, TAL_1, 0, n_cols - 1 , row);
443
444       for (count = 0 ; count < hsh_count(gp->group_hash) ; ++count)
445         {
446           gs = gs_array[count];
447
448           tab_text (t, 1, row + count, 
449                     TAB_LEFT | TAT_TITLE ,value_to_string(&gs->id,indep_var));
450
451           /* Now fill in the numbers ... */
452
453           tab_float (t, 2, row + count, 0, gs->n, 8,0);
454
455           tab_float (t, 3, row + count, 0, gs->mean,8,2);
456           
457           tab_float (t, 4, row + count, 0, gs->std_dev,8,2);
458
459           std_error = gs->std_dev/sqrt(gs->n) ;
460           tab_float (t, 5, row + count, 0, 
461                      std_error, 8,2);
462
463           /* Now the confidence interval */
464       
465           T = gsl_cdf_tdist_Qinv(q,gs->n - 1);
466
467           tab_float(t, 6, row + count, 0,
468                     gs->mean - T * std_error, 8, 2); 
469
470           tab_float(t, 7, row + count, 0,
471                     gs->mean + T * std_error, 8, 2); 
472
473           /* Min and Max */
474
475           tab_float(t, 8, row + count, 0,  gs->minimum, 8, 2); 
476           tab_float(t, 9, row + count, 0,  gs->maximum, 8, 2); 
477
478         }
479
480       tab_text (t, 1, row + count, 
481                 TAB_LEFT | TAT_TITLE ,_("Total"));
482
483       tab_float (t, 2, row + count, 0, totals->n, 8,0);
484
485       tab_float (t, 3, row + count, 0, totals->mean, 8,2);
486
487       tab_float (t, 4, row + count, 0, totals->std_dev,8,2);
488
489       std_error = totals->std_dev/sqrt(totals->n) ;
490
491       tab_float (t, 5, row + count, 0, std_error, 8,2);
492
493       /* Now the confidence interval */
494       
495       T = gsl_cdf_tdist_Qinv(q,totals->n - 1);
496
497       tab_float(t, 6, row + count, 0,
498                 totals->mean - T * std_error, 8, 2); 
499
500       tab_float(t, 7, row + count, 0,
501                 totals->mean + T * std_error, 8, 2); 
502
503       /* Min and Max */
504
505       tab_float(t, 8, row + count, 0,  totals->minimum, 8, 2); 
506       tab_float(t, 9, row + count, 0,  totals->maximum, 8, 2); 
507
508       row += gp->n_groups + 1;
509     }
510
511
512   tab_submit (t);
513
514
515 }
516
517 /* Show the homogeneity table */
518 static void 
519 show_homogeneity(void)
520 {
521   size_t v;
522   int n_cols = 5;
523   size_t n_rows = n_vars + 1;
524
525   struct tab_table *t;
526
527
528   t = tab_create (n_cols,n_rows,0);
529   tab_headers (t, 1, 0, 1, 0);
530   tab_dim (t, tab_natural_dimensions);
531
532   /* Put a frame around the entire box, and vertical lines inside */
533   tab_box (t, 
534            TAL_2, TAL_2,
535            -1, TAL_1,
536            0, 0,
537            n_cols - 1, n_rows - 1);
538
539
540   tab_hline(t, TAL_2, 0, n_cols - 1, 1);
541   tab_vline(t, TAL_2, 1, 0, n_rows - 1);
542
543
544   tab_text (t,  1, 0, TAB_CENTER | TAT_TITLE, _("Levene Statistic"));
545   tab_text (t,  2, 0, TAB_CENTER | TAT_TITLE, _("df1"));
546   tab_text (t,  3, 0, TAB_CENTER | TAT_TITLE, _("df2"));
547   tab_text (t,  4, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
548   
549
550   tab_title (t, _("Test of Homogeneity of Variances"));
551
552   for ( v=0 ; v < n_vars ; ++v ) 
553     {
554       double F;
555       const struct variable *var = vars[v];
556       const struct group_proc *gp = group_proc_get (vars[v]);
557       const char *s = var_to_string(var);
558       const struct group_statistics *totals = &gp->ugs;
559
560       const double df1 = gp->n_groups - 1;
561       const double df2 = totals->n - gp->n_groups ;
562
563       tab_text (t, 0, v + 1, TAB_LEFT | TAT_TITLE, s);
564
565       F = gp->levene;
566       tab_float (t, 1, v + 1, TAB_RIGHT, F, 8,3);
567       tab_float (t, 2, v + 1, TAB_RIGHT, df1 ,8,0);
568       tab_float (t, 3, v + 1, TAB_RIGHT, df2 ,8,0);
569
570       /* Now the significance */
571       tab_float (t, 4, v + 1, TAB_RIGHT,gsl_cdf_fdist_Q(F,df1,df2), 8, 3);
572     }
573
574   tab_submit (t);
575
576
577 }
578
579
580 /* Show the contrast coefficients table */
581 static void 
582 show_contrast_coeffs(short *bad_contrast)
583 {
584   int n_cols = 2 + ostensible_number_of_groups;
585   int n_rows = 2 + cmd.sbc_contrast;
586   union value *group_value;
587   int count = 0 ;      
588   void *const *group_values ;
589
590   struct tab_table *t;
591
592   t = tab_create (n_cols,n_rows,0);
593   tab_headers (t, 2, 0, 2, 0);
594   tab_dim (t, tab_natural_dimensions);
595
596   /* Put a frame around the entire box, and vertical lines inside */
597   tab_box (t, 
598            TAL_2, TAL_2,
599            -1, TAL_1,
600            0, 0,
601            n_cols - 1, n_rows - 1);
602
603   tab_box (t, 
604            -1,-1,
605            TAL_0, TAL_0,
606            2, 0,
607            n_cols - 1, 0);
608
609   tab_box (t,
610            -1,-1,
611            TAL_0, TAL_0,
612            0,0,
613            1,1);
614
615   tab_hline(t, TAL_1, 2, n_cols - 1, 1);
616   tab_hline(t, TAL_2, 0, n_cols - 1, 2);
617
618   tab_vline(t, TAL_2, 2, 0, n_rows - 1);
619
620   tab_title (t, _("Contrast Coefficients"));
621
622   tab_text (t,  0, 2, TAB_LEFT | TAT_TITLE, _("Contrast"));
623
624
625   tab_joint_text (t, 2, 0, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, 
626                   var_to_string(indep_var));
627
628   group_values = hsh_sort(global_group_hash);
629   for (count = 0 ; 
630        count < hsh_count(global_group_hash) ; 
631        ++count)
632     {
633       int i;
634       group_value = group_values[count];
635
636       tab_text (t, count + 2, 1, TAB_CENTER | TAT_TITLE, 
637                 value_to_string(group_value, indep_var));
638
639       for (i = 0 ; i < cmd.sbc_contrast ; ++i ) 
640         {
641           tab_text(t, 1, i + 2, TAB_CENTER | TAT_PRINTF, "%d", i + 1);
642
643           if ( bad_contrast[i] ) 
644             tab_text(t, count + 2, i + 2, TAB_RIGHT, "?" );
645           else
646             tab_text(t, count + 2, i + 2, TAB_RIGHT | TAT_PRINTF, "%g", 
647                      subc_list_double_at(&cmd.dl_contrast[i], count)
648                      );
649         }
650     }
651   
652   tab_submit (t);
653 }
654
655
656 /* Show the results of the contrast tests */
657 static void 
658 show_contrast_tests(short *bad_contrast)
659 {
660   size_t v;
661   int n_cols = 8;
662   size_t n_rows = 1 + n_vars * 2 * cmd.sbc_contrast;
663
664   struct tab_table *t;
665
666   t = tab_create (n_cols,n_rows,0);
667   tab_headers (t, 3, 0, 1, 0);
668   tab_dim (t, tab_natural_dimensions);
669
670   /* Put a frame around the entire box, and vertical lines inside */
671   tab_box (t, 
672            TAL_2, TAL_2,
673            -1, TAL_1,
674            0, 0,
675            n_cols - 1, n_rows - 1);
676
677   tab_box (t, 
678            -1,-1,
679            TAL_0, TAL_0,
680            0, 0,
681            2, 0);
682
683   tab_hline(t, TAL_2, 0, n_cols - 1, 1);
684   tab_vline(t, TAL_2, 3, 0, n_rows - 1);
685
686
687   tab_title (t, _("Contrast Tests"));
688
689   tab_text (t,  2, 0, TAB_CENTER | TAT_TITLE, _("Contrast"));
690   tab_text (t,  3, 0, TAB_CENTER | TAT_TITLE, _("Value of Contrast"));
691   tab_text (t,  4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
692   tab_text (t,  5, 0, TAB_CENTER | TAT_TITLE, _("t"));
693   tab_text (t,  6, 0, TAB_CENTER | TAT_TITLE, _("df"));
694   tab_text (t,  7, 0, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
695
696   for ( v = 0 ; v < n_vars ; ++v ) 
697     {
698       int i;
699       int lines_per_variable = 2 * cmd.sbc_contrast;
700
701
702       tab_text (t,  0, (v * lines_per_variable) + 1, TAB_LEFT | TAT_TITLE,
703                 var_to_string(vars[v]));
704
705       for ( i = 0 ; i < cmd.sbc_contrast ; ++i ) 
706         {
707           int ci;
708           double contrast_value = 0.0;
709           double coef_msq = 0.0;
710           struct group_proc *grp_data = group_proc_get (vars[v]);
711           struct hsh_table *group_hash = grp_data->group_hash;
712
713           void *const *group_stat_array;
714
715           double T;
716           double std_error_contrast ;
717           double df;
718           double sec_vneq=0.0;
719
720
721           /* Note: The calculation of the degrees of freedom in the 
722              "variances not equal" case is painfull!!
723              The following formula may help to understand it:
724              \frac{\left(\sum_{i=1}^k{c_i^2\frac{s_i^2}{n_i}}\right)^2}
725              {
726              \sum_{i=1}^k\left(
727              \frac{\left(c_i^2\frac{s_i^2}{n_i}\right)^2}  {n_i-1}
728              \right)
729              }
730           */
731
732           double df_denominator = 0.0;
733           double df_numerator = 0.0;
734           if ( i == 0 ) 
735             {
736               tab_text (t,  1, (v * lines_per_variable) + i + 1, 
737                         TAB_LEFT | TAT_TITLE,
738                         _("Assume equal variances"));
739
740               tab_text (t,  1, (v * lines_per_variable) + i + 1 + cmd.sbc_contrast, 
741                         TAB_LEFT | TAT_TITLE, 
742                         _("Does not assume equal"));
743             }
744
745           tab_text (t,  2, (v * lines_per_variable) + i + 1, 
746                     TAB_CENTER | TAT_TITLE | TAT_PRINTF, "%d",i+1);
747
748
749           tab_text (t,  2, (v * lines_per_variable) + i + 1 + cmd.sbc_contrast,
750                     TAB_CENTER | TAT_TITLE | TAT_PRINTF, "%d",i+1);
751
752
753           if ( bad_contrast[i]) 
754             continue;
755
756           group_stat_array = hsh_sort(group_hash);
757           
758           for (ci = 0 ; ci < hsh_count(group_hash) ;  ++ci)
759             {
760               const double coef = subc_list_double_at(&cmd.dl_contrast[i], ci);
761               struct group_statistics *gs = group_stat_array[ci];
762
763               const double winv = (gs->std_dev * gs->std_dev) / gs->n;
764
765               contrast_value += coef * gs->mean;
766
767               coef_msq += (coef * coef) / gs->n ; 
768
769               sec_vneq += (coef * coef) * (gs->std_dev * gs->std_dev ) /gs->n ;
770
771               df_numerator += (coef * coef) * winv;
772               df_denominator += pow2((coef * coef) * winv) / (gs->n - 1);
773             }
774           sec_vneq = sqrt(sec_vneq);
775
776           df_numerator = pow2(df_numerator);
777
778           tab_float (t,  3, (v * lines_per_variable) + i + 1, 
779                      TAB_RIGHT, contrast_value, 8,2);
780
781           tab_float (t,  3, (v * lines_per_variable) + i + 1 + 
782                      cmd.sbc_contrast,
783                      TAB_RIGHT, contrast_value, 8,2);
784
785           std_error_contrast = sqrt(grp_data->mse * coef_msq);
786
787           /* Std. Error */
788           tab_float (t,  4, (v * lines_per_variable) + i + 1, 
789                      TAB_RIGHT, std_error_contrast,
790                      8,3);
791
792           T = fabs(contrast_value / std_error_contrast) ;
793
794           /* T Statistic */
795
796           tab_float (t,  5, (v * lines_per_variable) + i + 1, 
797                      TAB_RIGHT, T,
798                      8,3);
799
800           df = grp_data->ugs.n - grp_data->n_groups;
801
802           /* Degrees of Freedom */
803           tab_float (t,  6, (v * lines_per_variable) + i + 1, 
804                      TAB_RIGHT,  df,
805                      8,0);
806
807
808           /* Significance TWO TAILED !!*/
809           tab_float (t,  7, (v * lines_per_variable) + i + 1, 
810                      TAB_RIGHT,  2 * gsl_cdf_tdist_Q(T,df),
811                      8,3);
812
813
814           /* Now for the Variances NOT Equal case */
815
816           /* Std. Error */
817           tab_float (t,  4, 
818                      (v * lines_per_variable) + i + 1 + cmd.sbc_contrast, 
819                      TAB_RIGHT, sec_vneq,
820                      8,3);
821
822
823           T = contrast_value / sec_vneq;
824           tab_float (t,  5, 
825                      (v * lines_per_variable) + i + 1 + cmd.sbc_contrast, 
826                      TAB_RIGHT, T,
827                      8,3);
828
829
830           df = df_numerator / df_denominator;
831
832           tab_float (t,  6, 
833                      (v * lines_per_variable) + i + 1 + cmd.sbc_contrast, 
834                      TAB_RIGHT, df,
835                      8,3);
836
837           /* The Significance */
838
839           tab_float (t, 7, (v * lines_per_variable) + i + 1 + cmd.sbc_contrast,
840                      TAB_RIGHT,  2 * gsl_cdf_tdist_Q(T,df),
841                      8,3);
842
843
844         }
845
846       if ( v > 0 ) 
847         tab_hline(t, TAL_1, 0, n_cols - 1, (v * lines_per_variable) + 1);
848     }
849
850   tab_submit (t);
851
852 }
853
854
855 /* ONEWAY ANOVA Calculations */
856
857 static void  postcalc (  struct cmd_oneway *cmd UNUSED );
858
859 static void  precalc ( struct cmd_oneway *cmd UNUSED );
860
861
862
863 /* Pre calculations */
864 static void 
865 precalc ( struct cmd_oneway *cmd UNUSED )
866 {
867   size_t i=0;
868
869   for(i=0; i< n_vars ; ++i) 
870     {
871       struct group_proc *gp = group_proc_get (vars[i]);
872       struct group_statistics *totals = &gp->ugs;
873       
874       /* Create a hash for each of the dependent variables.
875          The hash contains a group_statistics structure, 
876          and is keyed by value of the independent variable */
877
878       gp->group_hash = 
879         hsh_create(4, 
880                    (hsh_compare_func *) compare_group,
881                    (hsh_hash_func *) hash_group,
882                    (hsh_free_func *) free_group,
883                    (void *) indep_var->width );
884
885
886       totals->sum=0;
887       totals->n=0;
888       totals->ssq=0;
889       totals->sum_diff=0;
890       totals->maximum = - DBL_MAX;
891       totals->minimum = DBL_MAX;
892     }
893 }
894
895
896 static bool
897 run_oneway(const struct ccase *first, const struct casefile *cf, void *cmd_)
898 {
899   struct casereader *r;
900   struct ccase c;
901
902   struct cmd_oneway *cmd = (struct cmd_oneway *) cmd_;
903
904   output_split_file_values (first);
905
906   global_group_hash = hsh_create(4, 
907                                  (hsh_compare_func *) compare_values,
908                                  (hsh_hash_func *) hash_value,
909                                  0,
910                                  (void *) indep_var->width );
911   precalc(cmd);
912
913   for(r = casefile_get_reader (cf);
914       casereader_read (r, &c) ;
915       case_destroy (&c)) 
916     {
917       size_t i;
918
919       const double weight = 
920         dict_get_case_weight(default_dict,&c,&bad_weight_warn);
921       
922       const union value *indep_val = case_data (&c, indep_var->fv);
923
924       /* Deal with missing values */
925       if ( value_is_missing(&indep_var->miss, indep_val) )
926         continue;
927
928       /* Skip the entire case if /MISSING=LISTWISE is set */
929       if ( cmd->miss == ONEWAY_LISTWISE ) 
930         {
931           for(i = 0; i < n_vars ; ++i) 
932             {
933               const struct variable *v = vars[i];
934               const union value *val = case_data (&c, v->fv);
935
936               if (value_is_missing(&v->miss, val) )
937                 break;
938             }
939           if ( i != n_vars ) 
940             continue;
941
942         }
943       
944           
945       hsh_insert ( global_group_hash, (void *) indep_val );
946
947       for ( i = 0 ; i < n_vars ; ++i ) 
948         {
949           const struct variable *v = vars[i];
950
951           const union value *val = case_data (&c, v->fv);
952
953           struct group_proc *gp = group_proc_get (vars[i]);
954           struct hsh_table *group_hash = gp->group_hash;
955
956           struct group_statistics *gs;
957
958           gs = hsh_find(group_hash, (void *) indep_val );
959
960           if ( ! gs ) 
961             {
962               gs = xmalloc (sizeof *gs);
963               gs->id = *indep_val;
964               gs->sum=0;
965               gs->n=0;
966               gs->ssq=0;
967               gs->sum_diff=0;
968               gs->minimum = DBL_MAX;
969               gs->maximum = -DBL_MAX;
970
971               hsh_insert ( group_hash, (void *) gs );
972             }
973           
974           if (! value_is_missing(&v->miss, val) )
975             {
976               struct group_statistics *totals = &gp->ugs;
977
978               totals->n+=weight;
979               totals->sum+=weight * val->f;
980               totals->ssq+=weight * val->f * val->f;
981
982               if ( val->f * weight  < totals->minimum ) 
983                 totals->minimum = val->f * weight;
984
985               if ( val->f * weight  > totals->maximum ) 
986                 totals->maximum = val->f * weight;
987
988               gs->n+=weight;
989               gs->sum+=weight * val->f;
990               gs->ssq+=weight * val->f * val->f;
991
992               if ( val->f * weight  < gs->minimum ) 
993                 gs->minimum = val->f * weight;
994
995               if ( val->f * weight  > gs->maximum ) 
996                 gs->maximum = val->f * weight;
997             }
998
999           gp->n_groups = hsh_count ( group_hash );
1000         }
1001   
1002     }
1003   casereader_destroy (r);
1004
1005   postcalc(cmd);
1006
1007   
1008   if ( stat_tables & STAT_HOMO ) 
1009     levene(cf, indep_var, n_vars, vars, 
1010            (cmd->miss == ONEWAY_LISTWISE) ? LEV_LISTWISE : LEV_ANALYSIS ,
1011            value_is_missing);
1012
1013   ostensible_number_of_groups = hsh_count (global_group_hash);
1014
1015
1016   output_oneway();
1017
1018   return true;
1019 }
1020
1021
1022 /* Post calculations for the ONEWAY command */
1023 void 
1024 postcalc (  struct cmd_oneway *cmd UNUSED )
1025 {
1026   size_t i=0;
1027
1028
1029   for(i = 0; i < n_vars ; ++i) 
1030     {
1031       struct group_proc *gp = group_proc_get (vars[i]);
1032       struct hsh_table *group_hash = gp->group_hash;
1033       struct group_statistics *totals = &gp->ugs;
1034
1035       struct hsh_iterator g;
1036       struct group_statistics *gs;
1037
1038       for (gs =  hsh_first (group_hash,&g); 
1039            gs != 0; 
1040            gs = hsh_next(group_hash,&g))
1041         {
1042           gs->mean=gs->sum / gs->n;
1043           gs->s_std_dev= sqrt(
1044                               ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1045                               ) ;
1046
1047           gs->std_dev= sqrt(
1048                             gs->n/(gs->n-1) *
1049                             ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1050                             ) ;
1051
1052           gs->se_mean = gs->std_dev / sqrt(gs->n);
1053           gs->mean_diff= gs->sum_diff / gs->n;
1054
1055         }
1056
1057
1058
1059       totals->mean = totals->sum / totals->n;
1060       totals->std_dev= sqrt(
1061                             totals->n/(totals->n-1) *
1062                             ( (totals->ssq / totals->n ) - totals->mean * totals->mean )
1063                             ) ;
1064
1065       totals->se_mean = totals->std_dev / sqrt(totals->n);
1066         
1067     }
1068 }