d0fff3fd0cb4505ecdaa801bf248fbde829c6415
[pspp] / src / language / stats / t-test.q
1 /* PSPP - computes sample statistics. -*-c-*-
2
3    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
4    Written by John Williams <johnr.williams@stonebow.otago.ac.nz>.
5    Almost completly re-written by John Darrington 2004
6
7    This program is free software; you can redistribute it and/or
8    modify it under the terms of the GNU General Public License as
9    published by the Free Software Foundation; either version 2 of the
10    License, or (at your option) any later version.
11
12    This program is distributed in the hope that it will be useful, but
13    WITHOUT ANY WARRANTY; without even the implied warranty of
14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15    General Public License for more details.
16
17    You should have received a copy of the GNU General Public License
18    along with this program; if not, write to the Free Software
19    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20    02110-1301, USA. */
21
22 #include <config.h>
23 #include <gsl/gsl_cdf.h>
24 #include "message.h"
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <math.h>
28 #include "alloc.h"
29 #include "case.h"
30 #include "casefile.h"
31 #include "command.h"
32 #include "dictionary.h"
33 #include "message.h"
34 #include "group-proc.h"
35 #include "hash.h"
36 #include "levene.h"
37 #include "lexer.h"
38 #include "magic.h"
39 #include "misc.h"
40 #include "size_max.h"
41 #include "manager.h"
42 #include "str.h"
43 #include "table.h"
44 #include "value-labels.h"
45 #include "variable.h"
46 #include "procedure.h"
47
48 #include "gettext.h"
49 #define _(msgid) gettext (msgid)
50
51 /* (headers) */
52
53 /* (specification)
54    "T-TEST" (tts_):
55      +groups=custom;
56      testval=double;
57      variables=varlist("PV_NO_SCRATCH | PV_NUMERIC");
58      pairs=custom;
59      +missing=miss:!analysis/listwise,
60              incl:include/!exclude;
61      format=fmt:!labels/nolabels;
62      criteria=:cin(d:criteria,"%s > 0. && %s < 1.").
63 */
64 /* (declarations) */
65 /* (functions) */
66
67
68
69
70 /* Function to use for testing for missing values */
71 static is_missing_func *value_is_missing;
72
73 /* Variable for the GROUPS subcommand, if given. */
74 static struct variable *indep_var;
75
76 enum comparison
77   {
78     CMP_LE = -2,
79     CMP_EQ = 0,
80   };
81
82 struct group_properties
83 {
84   /* The comparison criterion */
85   enum comparison criterion;
86
87   /* The width of the independent variable */
88   int indep_width ;  
89
90   union {
91     /* The value of the independent variable at which groups are determined to 
92        belong to one group or the other */
93     double critical_value;
94     
95
96     /* The values of the independent variable for each group */
97     union value g_value[2];
98   } v ;
99
100 };
101
102
103 static struct group_properties gp ;
104
105
106
107 /* PAIRS: Number of pairs to be compared ; each pair. */
108 static int n_pairs = 0 ;
109 struct pair 
110 {
111   /* The variables comprising the pair */
112   struct variable *v[2];
113
114   /* The number of valid variable pairs */
115   double n;
116
117   /* The sum of the members */
118   double sum[2];
119
120   /* sum of squares of the members */
121   double ssq[2];
122
123   /* Std deviation of the members */
124   double std_dev[2];
125
126
127   /* Sample Std deviation of the members */
128   double s_std_dev[2];
129
130   /* The means of the members */
131   double mean[2];
132
133   /* The correlation coefficient between the variables */
134   double correlation;
135
136   /* The sum of the differences */
137   double sum_of_diffs;
138
139   /* The sum of the products */
140   double sum_of_prod;
141
142   /* The mean of the differences */
143   double mean_diff;
144
145   /* The sum of the squares of the differences */
146   double ssq_diffs;
147
148   /* The std deviation of the differences */
149   double std_dev_diff;
150 };
151
152 static struct pair *pairs=0;
153
154 static int parse_value (union value * v, int type) ;
155
156 /* Structures and Functions for the Statistics Summary Box */
157 struct ssbox;
158 typedef void populate_ssbox_func(struct ssbox *ssb,
159                                             struct cmd_t_test *cmd);
160 typedef void finalize_ssbox_func(struct ssbox *ssb);
161
162 struct ssbox
163 {
164   struct tab_table *t;
165
166   populate_ssbox_func *populate;
167   finalize_ssbox_func *finalize;
168
169 };
170
171 /* Create a ssbox */
172 void ssbox_create(struct ssbox *ssb,   struct cmd_t_test *cmd, int mode);
173
174 /* Populate a ssbox according to cmd */
175 void ssbox_populate(struct ssbox *ssb, struct cmd_t_test *cmd);
176
177 /* Submit and destroy a ssbox */
178 void ssbox_finalize(struct ssbox *ssb);
179
180 /* A function to create, populate and submit the Paired Samples Correlation 
181    box */
182 void pscbox(void);
183
184
185 /* Structures and Functions for the Test Results Box */
186 struct trbox;
187
188 typedef void populate_trbox_func(struct trbox *trb,
189                                  struct cmd_t_test *cmd);
190 typedef void finalize_trbox_func(struct trbox *trb);
191
192 struct trbox {
193   struct tab_table *t;
194   populate_trbox_func *populate;
195   finalize_trbox_func *finalize;
196 };
197
198 /* Create a trbox */
199 void trbox_create(struct trbox *trb,   struct cmd_t_test *cmd, int mode);
200
201 /* Populate a ssbox according to cmd */
202 void trbox_populate(struct trbox *trb, struct cmd_t_test *cmd);
203
204 /* Submit and destroy a ssbox */
205 void trbox_finalize(struct trbox *trb);
206
207 /* Which mode was T-TEST invoked */
208 enum {
209   T_1_SAMPLE = 0 ,
210   T_IND_SAMPLES, 
211   T_PAIRED
212 };
213
214
215 static int common_calc (const struct ccase *, void *);
216 static void common_precalc (struct cmd_t_test *);
217 static void common_postcalc (struct cmd_t_test *);
218
219 static int one_sample_calc (const struct ccase *, void *);
220 static void one_sample_precalc (struct cmd_t_test *);
221 static void one_sample_postcalc (struct cmd_t_test *);
222
223 static int  paired_calc (const struct ccase *, void *);
224 static void paired_precalc (struct cmd_t_test *);
225 static void paired_postcalc (struct cmd_t_test *);
226
227 static void group_precalc (struct cmd_t_test *);
228 static int  group_calc (const struct ccase *, struct cmd_t_test *);
229 static void group_postcalc (struct cmd_t_test *);
230
231
232 static bool calculate(const struct casefile *cf, void *_mode);
233
234 static  int mode;
235
236 static struct cmd_t_test cmd;
237
238 static int bad_weight_warn;
239
240
241 static int compare_group_binary(const struct group_statistics *a, 
242                                 const struct group_statistics *b, 
243                                 const struct group_properties *p);
244
245
246 static unsigned  hash_group_binary(const struct group_statistics *g, 
247                                    const struct group_properties *p);
248
249
250
251 int
252 cmd_t_test(void)
253 {
254   bool ok;
255   
256   if ( !parse_t_test(&cmd) )
257     return CMD_FAILURE;
258
259   if (! cmd.sbc_criteria)
260     cmd.criteria=0.95;
261
262   {
263     int m=0;
264     if (cmd.sbc_testval) ++m;
265     if (cmd.sbc_groups) ++m;
266     if (cmd.sbc_pairs) ++m;
267
268     if ( m != 1)
269       {
270         msg(SE, 
271             _("TESTVAL, GROUPS and PAIRS subcommands are mutually exclusive.")
272             );
273         free_t_test(&cmd);
274         return CMD_FAILURE;
275       }
276   }
277
278   if (cmd.sbc_testval) 
279     mode=T_1_SAMPLE;
280   else if (cmd.sbc_groups)
281     mode=T_IND_SAMPLES;
282   else
283     mode=T_PAIRED;
284
285   if ( mode == T_PAIRED) 
286     {
287       if (cmd.sbc_variables) 
288         {
289           msg(SE, _("VARIABLES subcommand is not appropriate with PAIRS"));
290           free_t_test(&cmd);
291           return CMD_FAILURE;
292         }
293       else
294         {
295           /* Iterate through the pairs and put each variable that is a 
296              member of a pair into cmd.v_variables */
297
298           int i;
299           struct hsh_iterator hi;
300           struct hsh_table *hash;
301           struct variable *v;
302
303           hash = hsh_create (n_pairs, compare_var_names, hash_var_name, 0, 0);
304
305           for (i=0; i < n_pairs; ++i)
306             {
307               hsh_insert(hash,pairs[i].v[0]);
308               hsh_insert(hash,pairs[i].v[1]);
309             }
310
311           assert(cmd.n_variables == 0);
312           cmd.n_variables = hsh_count(hash);
313
314           cmd.v_variables = xnrealloc (cmd.v_variables, cmd.n_variables,
315                                        sizeof *cmd.v_variables);
316           /* Iterate through the hash */
317           for (i=0,v = (struct variable *) hsh_first(hash,&hi);
318                v != 0;
319                v=hsh_next(hash,&hi) ) 
320             cmd.v_variables[i++]=v;
321
322           hsh_destroy(hash);
323         }
324     }
325   else if ( !cmd.sbc_variables) 
326     {
327       msg(SE, _("One or more VARIABLES must be specified."));
328       free_t_test(&cmd);
329       return CMD_FAILURE;
330     }
331
332
333   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
334   if (cmd.incl == TTS_INCLUDE ) 
335     value_is_missing = mv_is_value_system_missing;
336   else
337     value_is_missing = mv_is_value_missing;
338
339   bad_weight_warn = 1;
340
341   ok = multipass_procedure_with_splits (calculate, &cmd);
342
343   n_pairs=0;
344   free(pairs);
345   pairs=0;
346
347   if ( mode == T_IND_SAMPLES) 
348     {
349       int v;
350       /* Destroy any group statistics we created */
351       for (v = 0 ; v < cmd.n_variables ; ++v ) 
352         {
353           struct group_proc *grpp = group_proc_get (cmd.v_variables[v]);
354           hsh_destroy (grpp->group_hash);
355         }
356     }
357     
358   free_t_test(&cmd);
359   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
360 }
361
362 static int
363 tts_custom_groups (struct cmd_t_test *cmd UNUSED)
364 {
365   int n_group_values=0;
366
367   lex_match('=');
368
369   indep_var = parse_variable ();
370   if (!indep_var)
371     {
372       lex_error ("expecting variable name in GROUPS subcommand");
373       return 0;
374     }
375
376   if (indep_var->type == T_STRING && indep_var->width > MAX_SHORT_STRING)
377     {
378       msg (SE, _("Long string variable %s is not valid here."),
379            indep_var->name);
380       return 0;
381     }
382
383   if (!lex_match ('('))
384     {
385       if (indep_var->type == NUMERIC)
386         {
387           gp.v.g_value[0].f = 1;
388           gp.v.g_value[1].f = 2;
389
390           gp.criterion = CMP_EQ;
391           
392           n_group_values = 2;
393
394           return 1;
395         }
396       else
397         {
398           msg (SE, _("When applying GROUPS to a string variable, two "
399                      "values must be specified."));
400           return 0;
401         }
402     }
403
404   if (!parse_value (&gp.v.g_value[0], indep_var->type))
405       return 0;
406
407   lex_match (',');
408   if (lex_match (')'))
409     {
410       if (indep_var->type != NUMERIC)
411         {
412
413           msg (SE, _("When applying GROUPS to a string variable, two "
414                      "values must be specified."));
415           return 0;
416         }
417       gp.criterion = CMP_LE;
418       gp.v.critical_value = gp.v.g_value[0].f;
419
420       n_group_values = 1;
421       return 1;
422     }
423
424   if (!parse_value (&gp.v.g_value[1], indep_var->type))
425     return 0;
426
427   n_group_values = 2;
428   if (!lex_force_match (')'))
429     return 0;
430
431   if ( n_group_values == 2 ) 
432     gp.criterion = CMP_EQ ;
433   else
434     gp.criterion = CMP_LE ;
435
436
437   return 1;
438 }
439
440
441 static int
442 tts_custom_pairs (struct cmd_t_test *cmd UNUSED)
443 {
444   struct variable **vars;
445   size_t n_vars;
446   size_t n_pairs_local;
447
448   size_t n_before_WITH;
449   size_t n_after_WITH = SIZE_MAX;
450   int paired ; /* Was the PAIRED keyword given ? */
451
452   lex_match('=');
453
454   n_vars=0;
455   if (!parse_variables (default_dict, &vars, &n_vars,
456                         PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH))
457     {
458       free (vars);
459       return 0;
460     }
461   assert (n_vars);
462
463   n_before_WITH = 0;
464   if (lex_match (T_WITH))
465     {
466       n_before_WITH = n_vars;
467       if (!parse_variables (default_dict, &vars, &n_vars,
468                             PV_DUPLICATE | PV_APPEND
469                             | PV_NUMERIC | PV_NO_SCRATCH))
470         {
471           free (vars);
472           return 0;
473         }
474       n_after_WITH = n_vars - n_before_WITH;
475     }
476
477   paired = (lex_match ('(') && lex_match_id ("PAIRED") && lex_match (')'));
478
479   /* Determine the number of pairs needed */
480   if (paired)
481     {
482       if (n_before_WITH != n_after_WITH)
483         {
484           free (vars);
485           msg (SE, _("PAIRED was specified but the number of variables "
486                      "preceding WITH (%d) did not match the number "
487                      "following (%d)."),
488                n_before_WITH, n_after_WITH );
489           return 0;
490         }
491       n_pairs_local = n_before_WITH;
492     }
493   else if (n_before_WITH > 0) /* WITH keyword given, but not PAIRED keyword */
494     {
495       n_pairs_local = n_before_WITH * n_after_WITH ;
496     }
497   else /* Neither WITH nor PAIRED keyword given */
498     {
499       if (n_vars < 2)
500         {
501           free (vars);
502           msg (SE, _("At least two variables must be specified "
503                      "on PAIRS."));
504           return 0;
505         }
506
507       /* how many ways can you pick 2 from n_vars ? */
508       n_pairs_local = n_vars * (n_vars - 1) / 2;
509     }
510
511
512   /* Allocate storage for the pairs */
513   pairs = xnrealloc (pairs, n_pairs + n_pairs_local, sizeof *pairs);
514
515   /* Populate the pairs with the appropriate variables */
516   if ( paired ) 
517     {
518       int i;
519
520       assert(n_pairs_local == n_vars / 2);
521       for (i = 0; i < n_pairs_local; ++i)
522         {
523           pairs[i].v[n_pairs] = vars[i];
524           pairs[i].v[n_pairs + 1] = vars[i + n_pairs_local];
525         }
526     }
527   else if (n_before_WITH > 0) /* WITH keyword given, but not PAIRED keyword */
528     {
529       int i,j;
530       size_t p = n_pairs;
531
532       for(i=0 ; i < n_before_WITH ; ++i ) 
533         {
534           for(j=0 ; j < n_after_WITH ; ++j)
535             {
536               pairs[p].v[0] = vars[i];
537               pairs[p].v[1] = vars[j+n_before_WITH];
538               ++p;
539             }
540         }
541     }
542   else /* Neither WITH nor PAIRED given */
543     {
544       size_t i,j;
545       size_t p=n_pairs;
546       
547       for(i=0 ; i < n_vars ; ++i ) 
548         {
549           for(j=i+1 ; j < n_vars ; ++j)
550             {
551               pairs[p].v[0] = vars[i];
552               pairs[p].v[1] = vars[j];
553               ++p;
554             }
555         }
556     }
557
558   n_pairs+=n_pairs_local;
559
560   free (vars);
561   return 1;
562 }
563
564 /* Parses the current token (numeric or string, depending on type)
565     value v and returns success. */
566 static int
567 parse_value (union value * v, int type )
568 {
569   if (type == NUMERIC)
570     {
571       if (!lex_force_num ())
572         return 0;
573       v->f = tokval;
574     }
575   else
576     {
577       if (!lex_force_string ())
578         return 0;
579       strncpy (v->s, ds_c_str (&tokstr), ds_length (&tokstr));
580     }
581
582   lex_get ();
583
584   return 1;
585 }
586
587
588 /* Implementation of the SSBOX object */
589
590 void ssbox_base_init(struct ssbox *this, int cols,int rows);
591
592 void ssbox_base_finalize(struct ssbox *ssb);
593
594 void ssbox_one_sample_init(struct ssbox *this, 
595                            struct cmd_t_test *cmd );
596
597 void ssbox_independent_samples_init(struct ssbox *this,
598                                     struct cmd_t_test *cmd);
599
600 void ssbox_paired_init(struct ssbox *this,
601                            struct cmd_t_test *cmd);
602
603
604 /* Factory to create an ssbox */
605 void 
606 ssbox_create(struct ssbox *ssb, struct cmd_t_test *cmd, int mode)
607 {
608     switch (mode) 
609       {
610       case T_1_SAMPLE:
611         ssbox_one_sample_init(ssb,cmd);
612         break;
613       case T_IND_SAMPLES:
614         ssbox_independent_samples_init(ssb,cmd);
615         break;
616       case T_PAIRED:
617         ssbox_paired_init(ssb,cmd);
618         break;
619       default:
620         assert(0);
621       }
622 }
623
624
625
626 /* Despatcher for the populate method */
627 void
628 ssbox_populate(struct ssbox *ssb,struct cmd_t_test *cmd)
629 {
630   ssb->populate(ssb,cmd);
631 }
632
633
634 /* Despatcher for finalize */
635 void
636 ssbox_finalize(struct ssbox *ssb)
637 {
638   ssb->finalize(ssb);
639 }
640
641
642 /* Submit the box and clear up */
643 void 
644 ssbox_base_finalize(struct ssbox *ssb)
645 {
646   tab_submit(ssb->t);
647 }
648
649
650
651 /* Initialize a ssbox struct */
652 void 
653 ssbox_base_init(struct ssbox *this, int cols,int rows)
654 {
655   this->finalize = ssbox_base_finalize;
656   this->t = tab_create (cols, rows, 0);
657
658   tab_columns (this->t, SOM_COL_DOWN, 1);
659   tab_headers (this->t,0,0,1,0); 
660   tab_box (this->t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols -1, rows -1 );
661   tab_hline(this->t, TAL_2,0,cols-1,1);
662   tab_dim (this->t, tab_natural_dimensions);
663 }
664
665 void  ssbox_one_sample_populate(struct ssbox *ssb,
666                               struct cmd_t_test *cmd);
667
668 /* Initialize the one_sample ssbox */
669 void 
670 ssbox_one_sample_init(struct ssbox *this, 
671                            struct cmd_t_test *cmd )
672 {
673   const int hsize=5;
674   const int vsize=cmd->n_variables+1;
675
676   this->populate = ssbox_one_sample_populate;
677
678   ssbox_base_init(this, hsize,vsize);
679   tab_title (this->t, 0, _("One-Sample Statistics"));
680   tab_vline(this->t, TAL_2, 1,0,vsize - 1);
681   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE, _("N"));
682   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
683   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
684   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("SE. Mean"));
685 }
686
687 void ssbox_independent_samples_populate(struct ssbox *ssb,
688                                         struct cmd_t_test *cmd);
689
690 /* Initialize the independent samples ssbox */
691 void 
692 ssbox_independent_samples_init(struct ssbox *this, 
693         struct cmd_t_test *cmd)
694 {
695   int hsize=6;
696   int vsize = cmd->n_variables*2 +1;
697
698   this->populate = ssbox_independent_samples_populate;
699
700   ssbox_base_init(this, hsize,vsize);
701   tab_title (this->t, 0, _("Group Statistics"));
702   tab_vline(this->t,0,1,0,vsize - 1);
703   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE, indep_var->name);
704   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("N"));
705   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
706   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
707   tab_text (this->t, 5, 0, TAB_CENTER | TAT_TITLE, _("SE. Mean"));
708 }
709
710
711 /* Populate the ssbox for independent samples */
712 void 
713 ssbox_independent_samples_populate(struct ssbox *ssb,
714                               struct cmd_t_test *cmd)
715 {
716   int i;
717
718   char *val_lab0=0;
719   char *val_lab1=0;
720   double indep_value[2];
721
722   char prefix[2][3]={"",""};
723
724   if ( indep_var->type == NUMERIC ) 
725     {
726       val_lab0 = val_labs_find( indep_var->val_labs,gp.v.g_value[0]); 
727       val_lab1 = val_labs_find( indep_var->val_labs,gp.v.g_value[1]);
728     }
729   else
730     {
731       val_lab0 = gp.v.g_value[0].s;
732       val_lab1 = gp.v.g_value[1].s;
733     }
734
735   if (gp.criterion == CMP_LE ) 
736     {
737       strcpy(prefix[0],"< ");
738       strcpy(prefix[1],">=");
739       indep_value[0] = gp.v.critical_value;
740       indep_value[1] = gp.v.critical_value;
741     }
742   else
743     {
744       indep_value[0] = gp.v.g_value[0].f;
745       indep_value[1] = gp.v.g_value[1].f;
746     }
747
748   assert(ssb->t);
749
750   for (i=0; i < cmd->n_variables; ++i)
751     {
752       struct variable *var = cmd->v_variables[i];
753       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
754       int count=0;
755
756       tab_text (ssb->t, 0, i*2+1, TAB_LEFT, cmd->v_variables[i]->name);
757
758       if (val_lab0)
759         tab_text (ssb->t, 1, i*2+1, TAB_LEFT | TAT_PRINTF, 
760                   "%s%s", prefix[0], val_lab0);
761       else
762           tab_text (ssb->t, 1, i*2+1, TAB_LEFT | TAT_PRINTF, 
763                     "%s%g", prefix[0], indep_value[0]);
764
765
766       if (val_lab1)
767         tab_text (ssb->t, 1, i*2+1+1, TAB_LEFT | TAT_PRINTF, 
768                   "%s%s", prefix[1], val_lab1);
769       else
770           tab_text (ssb->t, 1, i*2+1+1, TAB_LEFT | TAT_PRINTF, 
771                     "%s%g", prefix[1], indep_value[1]);
772
773
774       /* Fill in the group statistics */
775       for ( count = 0 ; count < 2 ; ++count ) 
776         {
777           union value search_val;
778
779           struct group_statistics *gs;
780
781           if ( gp.criterion == CMP_LE ) 
782             {
783               if ( count == 0 ) 
784                 {
785                   /*  less than ( < )  case */
786                   search_val.f = gp.v.critical_value - 1.0;
787                 }
788               else
789                 {
790                   /* >= case  */
791                   search_val.f = gp.v.critical_value + 1.0;
792                 }
793             }
794           else
795             {
796               search_val = gp.v.g_value[count];
797             }
798
799           gs = hsh_find(grp_hash, (void *) &search_val);
800           assert(gs);
801
802           tab_float(ssb->t, 2 ,i*2+count+1, TAB_RIGHT, gs->n, 2, 0);
803           tab_float(ssb->t, 3 ,i*2+count+1, TAB_RIGHT, gs->mean, 8, 2);
804           tab_float(ssb->t, 4 ,i*2+count+1, TAB_RIGHT, gs->std_dev, 8, 3);
805           tab_float(ssb->t, 5 ,i*2+count+1, TAB_RIGHT, gs->se_mean, 8, 3);
806         }
807     }
808 }
809
810
811 void ssbox_paired_populate(struct ssbox *ssb,
812                            struct cmd_t_test *cmd);
813
814 /* Initialize the paired values ssbox */
815 void 
816 ssbox_paired_init(struct ssbox *this, struct cmd_t_test *cmd UNUSED)
817 {
818   int hsize=6;
819
820   int vsize = n_pairs*2+1;
821
822   this->populate = ssbox_paired_populate;
823
824   ssbox_base_init(this, hsize,vsize);
825   tab_title (this->t, 0, _("Paired Sample Statistics"));
826   tab_vline(this->t,TAL_0,1,0,vsize-1);
827   tab_vline(this->t,TAL_2,2,0,vsize-1);
828   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
829   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("N"));
830   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
831   tab_text (this->t, 5, 0, TAB_CENTER | TAT_TITLE, _("SE. Mean"));
832 }
833
834
835 /* Populate the ssbox for paired values */
836 void 
837 ssbox_paired_populate(struct ssbox *ssb,struct cmd_t_test *cmd UNUSED)
838 {
839   int i;
840
841   assert(ssb->t);
842
843   for (i=0; i < n_pairs; ++i)
844     {
845       int j;
846
847       tab_text (ssb->t, 0, i*2+1, TAB_LEFT | TAT_PRINTF , _("Pair %d"),i);
848
849       for (j=0 ; j < 2 ; ++j) 
850         {
851           struct group_statistics *gs;
852
853           gs = &group_proc_get (pairs[i].v[j])->ugs;
854
855           /* Titles */
856
857           tab_text (ssb->t, 1, i*2+j+1, TAB_LEFT, pairs[i].v[j]->name);
858
859           /* Values */
860           tab_float (ssb->t,2, i*2+j+1, TAB_RIGHT, pairs[i].mean[j], 8, 2);
861           tab_float (ssb->t,3, i*2+j+1, TAB_RIGHT, pairs[i].n, 2, 0);
862           tab_float (ssb->t,4, i*2+j+1, TAB_RIGHT, pairs[i].std_dev[j], 8, 3);
863           tab_float (ssb->t,5, i*2+j+1, TAB_RIGHT, pairs[i].std_dev[j]/sqrt(pairs[i].n), 8, 3);
864
865         }
866     }
867 }
868
869 /* Populate the one sample ssbox */
870 void 
871 ssbox_one_sample_populate(struct ssbox *ssb, struct cmd_t_test *cmd)
872 {
873   int i;
874
875   assert(ssb->t);
876
877   for (i=0; i < cmd->n_variables; ++i)
878     {
879       struct group_statistics *gs = &group_proc_get (cmd->v_variables[i])->ugs;
880
881       tab_text (ssb->t, 0, i+1, TAB_LEFT, cmd->v_variables[i]->name);
882       tab_float (ssb->t,1, i+1, TAB_RIGHT, gs->n, 2, 0);
883       tab_float (ssb->t,2, i+1, TAB_RIGHT, gs->mean, 8, 2);
884       tab_float (ssb->t,3, i+1, TAB_RIGHT, gs->std_dev, 8, 2);
885       tab_float (ssb->t,4, i+1, TAB_RIGHT, gs->se_mean, 8, 3);
886     }
887   
888 }
889
890
891
892 /* Implementation of the Test Results box struct */
893
894 void trbox_base_init(struct trbox *self,size_t n_vars, int cols);
895 void trbox_base_finalize(struct trbox *trb);
896
897 void trbox_independent_samples_init(struct trbox *trb,
898                                     struct cmd_t_test *cmd );
899
900 void trbox_independent_samples_populate(struct trbox *trb,
901                                         struct cmd_t_test *cmd);
902
903 void trbox_one_sample_init(struct trbox *self,
904                       struct cmd_t_test *cmd );
905
906 void trbox_one_sample_populate(struct trbox *trb,
907                                struct cmd_t_test *cmd);
908
909 void trbox_paired_init(struct trbox *self,
910                        struct cmd_t_test *cmd );
911
912 void trbox_paired_populate(struct trbox *trb,
913                       struct cmd_t_test *cmd);
914
915
916
917 /* Create a trbox according to mode*/
918 void 
919 trbox_create(struct trbox *trb,   
920              struct cmd_t_test *cmd, int mode)
921 {
922     switch (mode) 
923       {
924       case T_1_SAMPLE:
925         trbox_one_sample_init(trb,cmd);
926         break;
927       case T_IND_SAMPLES:
928         trbox_independent_samples_init(trb,cmd);
929         break;
930       case T_PAIRED:
931         trbox_paired_init(trb,cmd);
932         break;
933       default:
934         assert(0);
935       }
936 }
937
938 /* Populate a trbox according to cmd */
939 void 
940 trbox_populate(struct trbox *trb, struct cmd_t_test *cmd)
941 {
942   trb->populate(trb,cmd);
943 }
944
945 /* Submit and destroy a trbox */
946 void 
947 trbox_finalize(struct trbox *trb)
948 {
949   trb->finalize(trb);
950 }
951
952 /* Initialize the independent samples trbox */
953 void 
954 trbox_independent_samples_init(struct trbox *self,
955                            struct cmd_t_test *cmd UNUSED)
956 {
957   const int hsize=11;
958   const int vsize=cmd->n_variables*2+3;
959
960   assert(self);
961   self->populate = trbox_independent_samples_populate;
962
963   trbox_base_init(self,cmd->n_variables*2,hsize);
964   tab_title(self->t,0,_("Independent Samples Test"));
965   tab_hline(self->t,TAL_1,2,hsize-1,1);
966   tab_vline(self->t,TAL_2,2,0,vsize-1);
967   tab_vline(self->t,TAL_1,4,0,vsize-1);
968   tab_box(self->t,-1,-1,-1,TAL_1, 2,1,hsize-2,vsize-1);
969   tab_hline(self->t,TAL_1, hsize-2,hsize-1,2);
970   tab_box(self->t,-1,-1,-1,TAL_1, hsize-2,2,hsize-1,vsize-1);
971   tab_joint_text(self->t, 2, 0, 3, 0, 
972                  TAB_CENTER,_("Levene's Test for Equality of Variances"));
973   tab_joint_text(self->t, 4,0,hsize-1,0,
974                  TAB_CENTER,_("t-test for Equality of Means"));
975
976   tab_text(self->t,2,2, TAB_CENTER | TAT_TITLE,_("F"));
977   tab_text(self->t,3,2, TAB_CENTER | TAT_TITLE,_("Sig."));
978   tab_text(self->t,4,2, TAB_CENTER | TAT_TITLE,_("t"));
979   tab_text(self->t,5,2, TAB_CENTER | TAT_TITLE,_("df"));
980   tab_text(self->t,6,2, TAB_CENTER | TAT_TITLE,_("Sig. (2-tailed)"));
981   tab_text(self->t,7,2, TAB_CENTER | TAT_TITLE,_("Mean Difference"));
982   tab_text(self->t,8,2, TAB_CENTER | TAT_TITLE,_("Std. Error Difference"));
983   tab_text(self->t,9,2, TAB_CENTER | TAT_TITLE,_("Lower"));
984   tab_text(self->t,10,2, TAB_CENTER | TAT_TITLE,_("Upper"));
985
986   tab_joint_text(self->t, 9, 1, 10, 1, TAB_CENTER | TAT_PRINTF, 
987                  _("%g%% Confidence Interval of the Difference"),
988                  cmd->criteria*100.0);
989
990 }
991
992 /* Populate the independent samples trbox */
993 void 
994 trbox_independent_samples_populate(struct trbox *self,
995                                    struct cmd_t_test *cmd )
996 {
997   int i;
998
999   assert(self);
1000   for (i=0; i < cmd->n_variables; ++i)
1001     {
1002       double p,q;
1003
1004       double t;
1005       double df;
1006
1007       double df1, df2;
1008
1009       double pooled_variance;
1010       double std_err_diff;
1011       double mean_diff;
1012
1013       struct variable *var = cmd->v_variables[i];
1014       struct group_proc *grp_data = group_proc_get (var);
1015
1016       struct hsh_table *grp_hash = grp_data->group_hash;
1017
1018       struct group_statistics *gs0 ;
1019       struct group_statistics *gs1 ;
1020           
1021       union value search_val;
1022           
1023       if ( gp.criterion == CMP_LE ) 
1024         search_val.f = gp.v.critical_value - 1.0;
1025       else
1026         search_val = gp.v.g_value[0];
1027
1028       gs0 = hsh_find(grp_hash, (void *) &search_val);
1029       assert(gs0);
1030
1031       if ( gp.criterion == CMP_LE ) 
1032         search_val.f = gp.v.critical_value + 1.0;
1033       else
1034         search_val = gp.v.g_value[1];
1035
1036       gs1 = hsh_find(grp_hash, (void *) &search_val);
1037       assert(gs1);
1038
1039           
1040       tab_text (self->t, 0, i*2+3, TAB_LEFT, cmd->v_variables[i]->name);
1041
1042       tab_text (self->t, 1, i*2+3, TAB_LEFT, _("Equal variances assumed"));
1043
1044
1045       tab_float(self->t, 2, i*2+3, TAB_CENTER, grp_data->levene, 8,3);
1046
1047       /* Now work out the significance of the Levene test */
1048       df1 = 1; df2 = grp_data->ugs.n - 2;
1049       q = gsl_cdf_fdist_Q(grp_data->levene, df1, df2);
1050
1051       tab_float(self->t, 3, i*2+3, TAB_CENTER, q, 8,3 );
1052
1053       df = gs0->n + gs1->n - 2.0 ;
1054       tab_float (self->t, 5, i*2+3, TAB_RIGHT, df, 2, 0);
1055
1056       pooled_variance = ( (gs0->n )*pow2(gs0->s_std_dev)
1057                           + 
1058                           (gs1->n )*pow2(gs1->s_std_dev) 
1059                         ) / df  ;
1060
1061       t = (gs0->mean - gs1->mean) / sqrt(pooled_variance) ;
1062       t /= sqrt((gs0->n + gs1->n)/(gs0->n*gs1->n)); 
1063
1064       tab_float (self->t, 4, i*2+3, TAB_RIGHT, t, 8, 3);
1065
1066       p = gsl_cdf_tdist_P(t, df);
1067       q = gsl_cdf_tdist_Q(t, df);
1068
1069       tab_float(self->t, 6, i*2+3, TAB_RIGHT, 2.0*(t>0?q:p) , 8, 3);
1070
1071       mean_diff = gs0->mean - gs1->mean;
1072       tab_float(self->t, 7, i*2+3, TAB_RIGHT, mean_diff, 8, 3);
1073
1074
1075       std_err_diff = sqrt( pow2(gs0->se_mean) + pow2(gs1->se_mean));
1076       tab_float(self->t, 8, i*2+3, TAB_RIGHT, std_err_diff, 8, 3);
1077
1078
1079       /* Now work out the confidence interval */
1080       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1081
1082       t = gsl_cdf_tdist_Qinv(q,df);
1083       tab_float(self->t, 9, i*2+3, TAB_RIGHT, 
1084                 mean_diff - t * std_err_diff, 8, 3); 
1085
1086       tab_float(self->t, 10, i*2+3, TAB_RIGHT, 
1087                 mean_diff + t * std_err_diff, 8, 3); 
1088
1089
1090       {
1091         double se2;
1092       /* Now for the \sigma_1 != \sigma_2 case */
1093       tab_text (self->t, 1, i*2+3+1, 
1094                 TAB_LEFT, _("Equal variances not assumed"));
1095
1096
1097       se2 = (pow2(gs0->s_std_dev)/(gs0->n -1) ) +
1098         (pow2(gs1->s_std_dev)/(gs1->n -1) );
1099
1100       t = mean_diff / sqrt(se2) ;
1101       tab_float (self->t, 4, i*2+3+1, TAB_RIGHT, t, 8, 3);
1102                 
1103       df = pow2(se2) / ( 
1104                        (pow2(pow2(gs0->s_std_dev)/(gs0->n - 1 )) 
1105                         /(gs0->n -1 )
1106                         )
1107                        + 
1108                        (pow2(pow2(gs1->s_std_dev)/(gs1->n - 1 ))
1109                         /(gs1->n -1 )
1110                         )
1111                        ) ;
1112       tab_float (self->t, 5, i*2+3+1, TAB_RIGHT, df, 8, 3);
1113
1114       p = gsl_cdf_tdist_P(t, df);
1115       q = gsl_cdf_tdist_Q(t, df);
1116
1117       tab_float(self->t, 6, i*2+3+1, TAB_RIGHT, 2.0*(t>0?q:p) , 8, 3);
1118
1119       /* Now work out the confidence interval */
1120       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1121
1122       t = gsl_cdf_tdist_Qinv(q, df);
1123
1124       tab_float(self->t, 7, i*2+3+1, TAB_RIGHT, mean_diff, 8, 3);
1125
1126
1127       tab_float(self->t, 8, i*2+3+1, TAB_RIGHT, std_err_diff, 8, 3);
1128
1129
1130       tab_float(self->t, 9, i*2+3+1, TAB_RIGHT, 
1131                 mean_diff - t * std_err_diff, 8, 3); 
1132
1133       tab_float(self->t, 10, i*2+3+1, TAB_RIGHT, 
1134                 mean_diff + t * std_err_diff, 8, 3); 
1135
1136       }
1137     }
1138 }
1139
1140 /* Initialize the paired samples trbox */
1141 void 
1142 trbox_paired_init(struct trbox *self,
1143                            struct cmd_t_test *cmd UNUSED)
1144 {
1145
1146   const int hsize=10;
1147   const int vsize=n_pairs+3;
1148
1149   self->populate = trbox_paired_populate;
1150
1151   trbox_base_init(self,n_pairs,hsize);
1152   tab_title (self->t, 0, _("Paired Samples Test"));
1153   tab_hline(self->t,TAL_1,2,6,1);
1154   tab_vline(self->t,TAL_2,2,0,vsize - 1);
1155   tab_joint_text(self->t,2,0,6,0,TAB_CENTER,_("Paired Differences"));
1156   tab_box(self->t,-1,-1,-1,TAL_1, 2,1,6,vsize-1);
1157   tab_box(self->t,-1,-1,-1,TAL_1, 6,0,hsize-1,vsize-1);
1158   tab_hline(self->t,TAL_1,5,6, 2);
1159   tab_vline(self->t,TAL_0,6,0,1);
1160
1161   tab_joint_text(self->t, 5, 1, 6, 1, TAB_CENTER | TAT_PRINTF, 
1162                  _("%g%% Confidence Interval of the Difference"),
1163                  cmd->criteria*100.0);
1164
1165   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("Mean"));
1166   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
1167   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("Std. Error Mean"));
1168   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
1169   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
1170   tab_text (self->t, 7, 2, TAB_CENTER | TAT_TITLE, _("t"));
1171   tab_text (self->t, 8, 2, TAB_CENTER | TAT_TITLE, _("df"));
1172   tab_text (self->t, 9, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
1173 }
1174
1175 /* Populate the paired samples trbox */
1176 void 
1177 trbox_paired_populate(struct trbox *trb,
1178                               struct cmd_t_test *cmd UNUSED)
1179 {
1180   int i;
1181
1182   for (i=0; i < n_pairs; ++i)
1183     {
1184       double p,q;
1185       double se_mean;
1186
1187       double n = pairs[i].n;
1188       double t;
1189       double df = n - 1;
1190       
1191       tab_text (trb->t, 0, i+3, TAB_LEFT | TAT_PRINTF, _("Pair %d"),i); 
1192
1193       tab_text (trb->t, 1, i+3, TAB_LEFT | TAT_PRINTF, "%s - %s",
1194                 pairs[i].v[0]->name, pairs[i].v[1]->name);
1195
1196       tab_float(trb->t, 2, i+3, TAB_RIGHT, pairs[i].mean_diff, 8, 4);
1197
1198       tab_float(trb->t, 3, i+3, TAB_RIGHT, pairs[i].std_dev_diff, 8, 5);
1199
1200       /* SE Mean */
1201       se_mean = pairs[i].std_dev_diff / sqrt(n) ;
1202       tab_float(trb->t, 4, i+3, TAB_RIGHT, se_mean, 8,5 );
1203
1204       /* Now work out the confidence interval */
1205       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1206
1207       t = gsl_cdf_tdist_Qinv(q, df);
1208
1209       tab_float(trb->t, 5, i+3, TAB_RIGHT, 
1210                 pairs[i].mean_diff - t * se_mean , 8, 4); 
1211
1212       tab_float(trb->t, 6, i+3, TAB_RIGHT, 
1213                 pairs[i].mean_diff + t * se_mean , 8, 4); 
1214
1215       t = (pairs[i].mean[0] - pairs[i].mean[1])
1216         / sqrt (
1217                 ( pow2 (pairs[i].s_std_dev[0]) + pow2 (pairs[i].s_std_dev[1]) -
1218                   2 * pairs[i].correlation * 
1219                   pairs[i].s_std_dev[0] * pairs[i].s_std_dev[1] )
1220                 / (n - 1)
1221                 );
1222
1223       tab_float(trb->t, 7, i+3, TAB_RIGHT, t , 8,3 );
1224
1225       /* Degrees of freedom */
1226       tab_float(trb->t, 8, i+3, TAB_RIGHT, df , 2, 0 );
1227
1228       p = gsl_cdf_tdist_P(t,df);
1229       q = gsl_cdf_tdist_P(t,df);
1230
1231       tab_float(trb->t, 9, i+3, TAB_RIGHT, 2.0*(t>0?q:p) , 8, 3);
1232
1233     }
1234 }
1235
1236 /* Initialize the one sample trbox */
1237 void 
1238 trbox_one_sample_init(struct trbox *self, struct cmd_t_test *cmd )
1239 {
1240   const int hsize=7;
1241   const int vsize=cmd->n_variables+3;
1242
1243   self->populate = trbox_one_sample_populate;
1244
1245   trbox_base_init(self, cmd->n_variables,hsize);
1246   tab_title (self->t, 0, _("One-Sample Test"));
1247   tab_hline(self->t, TAL_1, 1, hsize - 1, 1);
1248   tab_vline(self->t, TAL_2, 1, 0, vsize - 1);
1249
1250   tab_joint_text(self->t, 1, 0, hsize-1,0, TAB_CENTER | TAT_PRINTF, 
1251                  _("Test Value = %f"), cmd->n_testval[0]);
1252
1253   tab_box(self->t, -1, -1, -1, TAL_1, 1,1,hsize-1,vsize-1);
1254
1255
1256   tab_joint_text(self->t,5,1,6,1,TAB_CENTER  | TAT_PRINTF, 
1257                  _("%g%% Confidence Interval of the Difference"),
1258                  cmd->criteria*100.0);
1259
1260   tab_vline(self->t,TAL_0,6,1,1);
1261   tab_hline(self->t,TAL_1,5,6,2);
1262   tab_text (self->t, 1, 2, TAB_CENTER | TAT_TITLE, _("t"));
1263   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("df"));
1264   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
1265   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
1266   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
1267   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
1268
1269 }
1270
1271
1272 /* Populate the one sample trbox */
1273 void 
1274 trbox_one_sample_populate(struct trbox *trb, struct cmd_t_test *cmd)
1275 {
1276   int i;
1277
1278   assert(trb->t);
1279
1280   for (i=0; i < cmd->n_variables; ++i)
1281     {
1282       double t;
1283       double p,q;
1284       double df;
1285       struct group_statistics *gs = &group_proc_get (cmd->v_variables[i])->ugs;
1286
1287
1288       tab_text (trb->t, 0, i+3, TAB_LEFT, cmd->v_variables[i]->name);
1289
1290       t = (gs->mean - cmd->n_testval[0] ) * sqrt(gs->n) / gs->std_dev ;
1291
1292       tab_float (trb->t, 1, i+3, TAB_RIGHT, t, 8,3);
1293
1294       /* degrees of freedom */
1295       df = gs->n - 1;
1296
1297       tab_float (trb->t, 2, i+3, TAB_RIGHT, df, 8,0);
1298
1299       p = gsl_cdf_tdist_P(t, df);
1300       q = gsl_cdf_tdist_Q(t, df);
1301
1302       /* Multiply by 2 to get 2-tailed significance, makeing sure we've got 
1303          the correct tail*/
1304       tab_float (trb->t, 3, i+3, TAB_RIGHT, 2.0*(t>0?q:p), 8,3);
1305
1306       tab_float (trb->t, 4, i+3, TAB_RIGHT, gs->mean_diff, 8,3);
1307
1308
1309       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1310       t = gsl_cdf_tdist_Qinv(q, df);
1311
1312       tab_float (trb->t, 5, i+3, TAB_RIGHT,
1313                  gs->mean_diff - t * gs->se_mean, 8,4);
1314
1315       tab_float (trb->t, 6, i+3, TAB_RIGHT,
1316                  gs->mean_diff + t * gs->se_mean, 8,4);
1317     }
1318 }
1319
1320 /* Base initializer for the generalized trbox */
1321 void 
1322 trbox_base_init(struct trbox *self, size_t data_rows, int cols)
1323 {
1324   const size_t rows = 3 + data_rows;
1325
1326   self->finalize = trbox_base_finalize;
1327   self->t = tab_create (cols, rows, 0);
1328   tab_headers (self->t,0,0,3,0); 
1329   tab_box (self->t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols -1, rows -1);
1330   tab_hline(self->t, TAL_2,0,cols-1,3);
1331   tab_dim (self->t, tab_natural_dimensions);
1332 }
1333
1334
1335 /* Base finalizer for the trbox */
1336 void 
1337 trbox_base_finalize(struct trbox *trb)
1338 {
1339   tab_submit(trb->t);
1340 }
1341
1342
1343 /* Create , populate and submit the Paired Samples Correlation box */
1344 void
1345 pscbox(void)
1346 {
1347   const int rows=1+n_pairs;
1348   const int cols=5;
1349   int i;
1350   
1351   struct tab_table *table;
1352   
1353   table = tab_create (cols,rows,0);
1354
1355   tab_columns (table, SOM_COL_DOWN, 1);
1356   tab_headers (table,0,0,1,0); 
1357   tab_box (table, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols -1, rows -1 );
1358   tab_hline(table, TAL_2, 0, cols - 1, 1);
1359   tab_vline(table, TAL_2, 2, 0, rows - 1);
1360   tab_dim(table, tab_natural_dimensions);
1361   tab_title(table, 0, _("Paired Samples Correlations"));
1362
1363   /* column headings */
1364   tab_text(table, 2,0, TAB_CENTER | TAT_TITLE, _("N"));
1365   tab_text(table, 3,0, TAB_CENTER | TAT_TITLE, _("Correlation"));
1366   tab_text(table, 4,0, TAB_CENTER | TAT_TITLE, _("Sig."));
1367
1368   for (i=0; i < n_pairs; ++i)
1369     {
1370       double p,q;
1371
1372       double df = pairs[i].n -2;
1373
1374       double correlation_t = 
1375         pairs[i].correlation * sqrt(df) /
1376         sqrt(1 - pow2(pairs[i].correlation));
1377
1378
1379       /* row headings */
1380       tab_text(table, 0,i+1, TAB_LEFT | TAT_TITLE | TAT_PRINTF, 
1381                _("Pair %d"), i);
1382       
1383       tab_text(table, 1,i+1, TAB_LEFT | TAT_TITLE | TAT_PRINTF, 
1384                _("%s & %s"), pairs[i].v[0]->name, pairs[i].v[1]->name);
1385
1386
1387       /* row data */
1388       tab_float(table, 2, i+1, TAB_RIGHT, pairs[i].n, 4, 0);
1389       tab_float(table, 3, i+1, TAB_RIGHT, pairs[i].correlation, 8, 3);
1390
1391       p = gsl_cdf_tdist_P(correlation_t, df);
1392       q = gsl_cdf_tdist_Q(correlation_t, df);
1393
1394       tab_float(table, 4, i+1, TAB_RIGHT, 2.0*(correlation_t>0?q:p), 8, 3);
1395     }
1396
1397   tab_submit(table);
1398 }
1399
1400
1401
1402
1403 /* Calculation Implementation */
1404
1405 /* Per case calculations common to all variants of the T test */
1406 static int 
1407 common_calc (const struct ccase *c, void *_cmd)
1408 {
1409   int i;
1410   struct cmd_t_test *cmd = (struct cmd_t_test *)_cmd;  
1411
1412   double weight = dict_get_case_weight(default_dict,c,&bad_weight_warn);
1413
1414
1415   /* Skip the entire case if /MISSING=LISTWISE is set */
1416   if ( cmd->miss == TTS_LISTWISE ) 
1417     {
1418       for(i=0; i< cmd->n_variables ; ++i) 
1419         {
1420           struct variable *v = cmd->v_variables[i];
1421           const union value *val = case_data (c, v->fv);
1422
1423           if (value_is_missing(&v->miss, val) )
1424             {
1425               return 0;
1426             }
1427         }
1428     }
1429
1430   /* Listwise has to be implicit if the independent variable is missing ?? */
1431   if ( cmd->sbc_groups )
1432     {
1433       const union value *gv = case_data (c, indep_var->fv);
1434       if ( value_is_missing(&indep_var->miss, gv) )
1435         {
1436           return 0;
1437         }
1438     }
1439
1440
1441   for(i=0; i< cmd->n_variables ; ++i) 
1442     {
1443       struct group_statistics *gs;
1444       struct variable *v = cmd->v_variables[i];
1445       const union value *val = case_data (c, v->fv);
1446
1447       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1448
1449       if (! value_is_missing(&v->miss, val) )
1450         {
1451           gs->n+=weight;
1452           gs->sum+=weight * val->f;
1453           gs->ssq+=weight * val->f * val->f;
1454         }
1455     }
1456   return 0;
1457 }
1458
1459 /* Pre calculations common to all variants of the T test */
1460 static void 
1461 common_precalc ( struct cmd_t_test *cmd )
1462 {
1463   int i=0;
1464
1465   for(i=0; i< cmd->n_variables ; ++i) 
1466     {
1467       struct group_statistics *gs;
1468       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1469       
1470       gs->sum=0;
1471       gs->n=0;
1472       gs->ssq=0;
1473       gs->sum_diff=0;
1474     }
1475 }
1476
1477 /* Post calculations common to all variants of the T test */
1478 void 
1479 common_postcalc (  struct cmd_t_test *cmd )
1480 {
1481   int i=0;
1482
1483
1484   for(i=0; i< cmd->n_variables ; ++i) 
1485     {
1486       struct group_statistics *gs;
1487       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1488       
1489       gs->mean=gs->sum / gs->n;
1490       gs->s_std_dev= sqrt(
1491                          ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1492                          ) ;
1493
1494       gs->std_dev= sqrt(
1495                          gs->n/(gs->n-1) *
1496                          ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1497                          ) ;
1498
1499       gs->se_mean = gs->std_dev / sqrt(gs->n);
1500       gs->mean_diff= gs->sum_diff / gs->n;
1501     }
1502 }
1503
1504 /* Per case calculations for one sample t test  */
1505 static int 
1506 one_sample_calc (const struct ccase *c, void *cmd_)
1507 {
1508   int i;
1509   struct cmd_t_test *cmd = (struct cmd_t_test *)cmd_;
1510
1511
1512   double weight = dict_get_case_weight(default_dict,c,&bad_weight_warn);
1513
1514   /* Skip the entire case if /MISSING=LISTWISE is set */
1515   if ( cmd->miss == TTS_LISTWISE ) 
1516     {
1517       for(i=0; i< cmd->n_variables ; ++i) 
1518         {
1519           struct variable *v = cmd->v_variables[i];
1520           const union value *val = case_data (c, v->fv);
1521
1522           if (value_is_missing(&v->miss, val) )
1523             {
1524               return 0;
1525             }
1526         }
1527     }
1528
1529   for(i=0; i< cmd->n_variables ; ++i) 
1530     {
1531       struct group_statistics *gs;
1532       struct variable *v = cmd->v_variables[i];
1533       const union value *val = case_data (c, v->fv);
1534
1535       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1536       
1537       if ( ! value_is_missing(&v->miss, val))
1538         gs->sum_diff += weight * (val->f - cmd->n_testval[0]);
1539     }
1540
1541   return 0;
1542 }
1543
1544 /* Pre calculations for one sample t test */
1545 static void 
1546 one_sample_precalc ( struct cmd_t_test *cmd )
1547 {
1548   int i=0; 
1549  
1550   for(i=0; i< cmd->n_variables ; ++i) 
1551     {
1552       struct group_statistics *gs;
1553       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1554       
1555       gs->sum_diff=0;
1556     }
1557 }
1558
1559 /* Post calculations for one sample t test */
1560 static void 
1561 one_sample_postcalc (struct cmd_t_test *cmd)
1562 {
1563   int i=0;
1564   
1565   for(i=0; i< cmd->n_variables ; ++i) 
1566     {
1567       struct group_statistics *gs;
1568       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1569
1570       gs->mean_diff = gs->sum_diff / gs->n ;
1571     }
1572 }
1573
1574
1575
1576 static void 
1577 paired_precalc (struct cmd_t_test *cmd UNUSED)
1578 {
1579   int i;
1580
1581   for(i=0; i < n_pairs ; ++i )
1582     {
1583       pairs[i].n = 0;
1584       pairs[i].sum[0] = 0;      pairs[i].sum[1] = 0;
1585       pairs[i].ssq[0] = 0;      pairs[i].ssq[1] = 0;
1586       pairs[i].sum_of_prod = 0;
1587       pairs[i].correlation = 0;
1588       pairs[i].sum_of_diffs = 0;
1589       pairs[i].ssq_diffs = 0;
1590     }
1591
1592 }
1593
1594
1595 static int  
1596 paired_calc (const struct ccase *c, void *cmd_)
1597 {
1598   int i;
1599
1600   struct cmd_t_test *cmd  = (struct cmd_t_test *) cmd_;
1601
1602   double weight = dict_get_case_weight(default_dict,c,&bad_weight_warn);
1603
1604   /* Skip the entire case if /MISSING=LISTWISE is set , 
1605    AND one member of a pair is missing */
1606   if ( cmd->miss == TTS_LISTWISE ) 
1607     {
1608       for(i=0; i < n_pairs ; ++i )
1609         {
1610           struct variable *v0 = pairs[i].v[0];
1611           struct variable *v1 = pairs[i].v[1];
1612
1613           const union value *val0 = case_data (c, v0->fv);
1614           const union value *val1 = case_data (c, v1->fv);
1615           
1616           if ( value_is_missing(&v0->miss, val0) ||
1617                value_is_missing(&v1->miss, val1) )
1618             {
1619               return 0;
1620             }
1621         }
1622     }
1623
1624   for(i=0; i < n_pairs ; ++i )
1625     {
1626       struct variable *v0 = pairs[i].v[0];
1627       struct variable *v1 = pairs[i].v[1];
1628
1629       const union value *val0 = case_data (c, v0->fv);
1630       const union value *val1 = case_data (c, v1->fv);
1631
1632       if ( ( !value_is_missing(&v0->miss, val0)
1633              && !value_is_missing(&v1->miss, val1) ) )
1634       {
1635         pairs[i].n += weight;
1636         pairs[i].sum[0] += weight * val0->f;
1637         pairs[i].sum[1] += weight * val1->f;
1638
1639         pairs[i].ssq[0] += weight * pow2(val0->f);
1640         pairs[i].ssq[1] += weight * pow2(val1->f);
1641
1642         pairs[i].sum_of_prod += weight * val0->f * val1->f ;
1643
1644         pairs[i].sum_of_diffs += weight * ( val0->f - val1->f ) ;
1645         pairs[i].ssq_diffs += weight * pow2(val0->f - val1->f);
1646       }
1647     }
1648
1649   return 0;
1650 }
1651
1652 static void 
1653 paired_postcalc (struct cmd_t_test *cmd UNUSED)
1654 {
1655   int i;
1656
1657   for(i=0; i < n_pairs ; ++i )
1658     {
1659       int j;
1660       const double n = pairs[i].n;
1661
1662       for (j=0; j < 2 ; ++j) 
1663         {
1664           pairs[i].mean[j] = pairs[i].sum[j] / n ;
1665           pairs[i].s_std_dev[j] = sqrt((pairs[i].ssq[j] / n - 
1666                                               pow2(pairs[i].mean[j]))
1667                                      );
1668
1669           pairs[i].std_dev[j] = sqrt(n/(n-1)*(pairs[i].ssq[j] / n - 
1670                                               pow2(pairs[i].mean[j]))
1671                                      );
1672         }
1673       
1674       pairs[i].correlation = pairs[i].sum_of_prod / pairs[i].n - 
1675         pairs[i].mean[0] * pairs[i].mean[1] ;
1676       /* correlation now actually contains the covariance */
1677       
1678       pairs[i].correlation /= pairs[i].std_dev[0] * pairs[i].std_dev[1];
1679       pairs[i].correlation *= pairs[i].n / ( pairs[i].n - 1 );
1680       
1681       pairs[i].mean_diff = pairs[i].sum_of_diffs / n ;
1682
1683       pairs[i].std_dev_diff = sqrt (  n / (n - 1) * (
1684                                     ( pairs[i].ssq_diffs / n )
1685                                     - 
1686                                     pow2(pairs[i].mean_diff )
1687                                     ) );
1688     }
1689 }
1690
1691 static void 
1692 group_precalc (struct cmd_t_test *cmd )
1693 {
1694   int i;
1695   int j;
1696
1697   for(i=0; i< cmd->n_variables ; ++i) 
1698     {
1699       struct group_proc *ttpr = group_proc_get (cmd->v_variables[i]);
1700
1701       /* There's always 2 groups for a T - TEST */
1702       ttpr->n_groups = 2;
1703
1704       gp.indep_width = indep_var->width;
1705       
1706       ttpr->group_hash = hsh_create(2, 
1707                                     (hsh_compare_func *) compare_group_binary,
1708                                     (hsh_hash_func *) hash_group_binary,
1709                                     (hsh_free_func *) free_group,
1710                                     (void *) &gp );
1711
1712       for (j=0 ; j < 2 ; ++j)
1713         {
1714
1715           struct group_statistics *gs = xmalloc (sizeof *gs);
1716
1717           gs->sum = 0;
1718           gs->n = 0;
1719           gs->ssq = 0;
1720         
1721           if ( gp.criterion == CMP_EQ ) 
1722             {
1723               gs->id = gp.v.g_value[j];
1724             }
1725           else
1726             {
1727               if ( j == 0 ) 
1728                 gs->id.f = gp.v.critical_value - 1.0 ;
1729               else
1730                 gs->id.f = gp.v.critical_value + 1.0 ;
1731             }
1732           
1733           hsh_insert ( ttpr->group_hash, (void *) gs );
1734
1735         }
1736     }
1737
1738 }
1739
1740 static int  
1741 group_calc (const struct ccase *c, struct cmd_t_test *cmd)
1742 {
1743   int i;
1744
1745   const union value *gv = case_data (c, indep_var->fv);
1746
1747   const double weight = dict_get_case_weight(default_dict,c,&bad_weight_warn);
1748
1749   if ( value_is_missing(&indep_var->miss, gv) )
1750     {
1751       return 0;
1752     }
1753
1754   if ( cmd->miss == TTS_LISTWISE ) 
1755     {
1756       for(i=0; i< cmd->n_variables ; ++i) 
1757         {
1758           struct variable *v = cmd->v_variables[i];
1759           const union value *val = case_data (c, v->fv);
1760
1761           if (value_is_missing(&v->miss, val) )
1762             {
1763               return 0;
1764             }
1765         }
1766     }
1767
1768   gv = case_data (c, indep_var->fv);
1769
1770   for(i=0; i< cmd->n_variables ; ++i) 
1771     {
1772       struct variable *var = cmd->v_variables[i];
1773       const union value *val = case_data (c, var->fv);
1774       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1775       struct group_statistics *gs;
1776
1777       gs = hsh_find(grp_hash, (void *) gv);
1778
1779       /* If the independent variable doesn't match either of the values 
1780          for this case then move on to the next case */
1781       if ( ! gs ) 
1782         return 0;
1783
1784       if ( !value_is_missing(&var->miss, val) )
1785         {
1786           gs->n+=weight;
1787           gs->sum+=weight * val->f;
1788           gs->ssq+=weight * pow2(val->f);
1789         }
1790     }
1791
1792   return 0;
1793 }
1794
1795
1796 static void 
1797 group_postcalc ( struct cmd_t_test *cmd )
1798 {
1799   int i;
1800
1801   for(i=0; i< cmd->n_variables ; ++i) 
1802     {
1803       struct variable *var = cmd->v_variables[i];
1804       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1805       struct hsh_iterator g;
1806       struct group_statistics *gs;
1807       int count=0;
1808
1809       for (gs =  hsh_first (grp_hash,&g); 
1810            gs != 0; 
1811            gs = hsh_next(grp_hash,&g))
1812         {
1813           gs->mean = gs->sum / gs->n;
1814           
1815           gs->s_std_dev= sqrt(
1816                               ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1817                               ) ;
1818
1819           gs->std_dev= sqrt(
1820                             gs->n/(gs->n-1) *
1821                             ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1822                             ) ;
1823           
1824           gs->se_mean = gs->std_dev / sqrt(gs->n);
1825           count ++;
1826         }
1827       assert(count == 2);
1828     }
1829 }
1830
1831
1832
1833 static bool
1834 calculate(const struct casefile *cf, void *cmd_)
1835 {
1836   struct ssbox stat_summary_box;
1837   struct trbox test_results_box;
1838
1839   struct casereader *r;
1840   struct ccase c;
1841
1842   struct cmd_t_test *cmd = (struct cmd_t_test *) cmd_;
1843
1844   common_precalc(cmd);
1845   for(r = casefile_get_reader (cf);
1846       casereader_read (r, &c) ;
1847       case_destroy (&c)) 
1848     {
1849       common_calc(&c,cmd);
1850     }
1851   casereader_destroy (r);
1852   common_postcalc(cmd);
1853
1854   switch(mode)
1855     {
1856     case T_1_SAMPLE:
1857       one_sample_precalc(cmd);
1858       for(r = casefile_get_reader (cf);
1859           casereader_read (r, &c) ;
1860           case_destroy (&c)) 
1861         {
1862           one_sample_calc(&c,cmd);
1863         }
1864       casereader_destroy (r);
1865       one_sample_postcalc(cmd);
1866
1867       break;
1868     case T_PAIRED:
1869       paired_precalc(cmd);
1870       for(r = casefile_get_reader (cf);
1871           casereader_read (r, &c) ;
1872           case_destroy (&c)) 
1873         {
1874           paired_calc(&c,cmd);
1875         }
1876       casereader_destroy (r);
1877       paired_postcalc(cmd);
1878
1879       break;
1880     case T_IND_SAMPLES:
1881
1882       group_precalc(cmd);
1883       for(r = casefile_get_reader (cf);
1884           casereader_read (r, &c) ;
1885           case_destroy (&c)) 
1886         {
1887           group_calc(&c,cmd);
1888         }
1889       casereader_destroy (r);
1890       group_postcalc(cmd);
1891
1892       levene(cf, indep_var, cmd->n_variables, cmd->v_variables,
1893              (cmd->miss == TTS_LISTWISE)?LEV_LISTWISE:LEV_ANALYSIS ,
1894              value_is_missing);
1895       break;
1896     }
1897
1898   ssbox_create(&stat_summary_box,cmd,mode);
1899   ssbox_populate(&stat_summary_box,cmd);
1900   ssbox_finalize(&stat_summary_box);
1901
1902   if ( mode == T_PAIRED) 
1903       pscbox();
1904
1905   trbox_create(&test_results_box,cmd,mode);
1906   trbox_populate(&test_results_box,cmd);
1907   trbox_finalize(&test_results_box);
1908
1909   return true;
1910 }
1911
1912 short which_group(const struct group_statistics *g,
1913                   const struct group_properties *p);
1914
1915 /* Return -1 if the id of a is less than b; +1 if greater than and 
1916    0 if equal */
1917 static int 
1918 compare_group_binary(const struct group_statistics *a, 
1919                      const struct group_statistics *b, 
1920                      const struct group_properties *p)
1921 {
1922   short flag_a;
1923   short flag_b;
1924   
1925   if ( p->criterion == CMP_LE ) 
1926     {
1927       /* less-than-or-equal comparision is not meaningfull for
1928          alpha variables, so we shouldn't ever arrive here */
1929       assert(p->indep_width == 0 ) ;
1930       
1931       flag_a = ( a->id.f < p->v.critical_value ) ;
1932       flag_b = ( b->id.f < p->v.critical_value ) ;
1933     }
1934   else
1935     {
1936       flag_a = which_group(a, p);
1937       flag_b = which_group(b, p);
1938     }
1939
1940   if (flag_a < flag_b ) 
1941     return -1;
1942
1943   return (flag_a > flag_b);
1944 }
1945
1946 /* This is a degenerate case of a hash, since it can only return three possible
1947    values.  It's really a comparison, being used as a hash function */
1948
1949 static unsigned 
1950 hash_group_binary(const struct group_statistics *g, 
1951                   const struct group_properties *p)
1952 {
1953   short flag = -1;
1954
1955   if ( p->criterion == CMP_LE ) 
1956     {
1957       /* Not meaningfull to do a less than compare for alpha values ? */
1958       assert(p->indep_width == 0 ) ;
1959       flag = ( g->id.f < p->v.critical_value ) ; 
1960     }
1961   else if ( p->criterion == CMP_EQ) 
1962     {
1963       flag = which_group(g,p);
1964     }
1965   else
1966     assert(0);
1967
1968   return flag;
1969 }
1970
1971 /* return 0 if G belongs to group 0, 
1972           1 if it belongs to group 1,
1973           2 if it belongs to neither group */
1974 short
1975 which_group(const struct group_statistics *g,
1976             const struct group_properties *p)
1977 {
1978  
1979   if ( 0 == compare_values (&g->id, &p->v.g_value[0], p->indep_width))
1980     return 0;
1981
1982   if ( 0 == compare_values (&g->id, &p->v.g_value[1], p->indep_width))
1983     return 1;
1984
1985   return 2;
1986 }
1987