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