src/language/stats/t-test.q

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