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