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