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 <language/lexer/value-parser.h>
  36 #include <libpspp/array.h>
  37 #include <libpspp/assertion.h>
  38 #include <libpspp/compiler.h>
  39 #include <libpspp/hash.h>
  40 #include <libpspp/message.h>
  41 #include <libpspp/misc.h>
  42 #include <libpspp/str.h>
  43 #include <libpspp/taint.h>
  44 #include <math/group-proc.h>
  45 #include <math/levene.h>
  46 #include <output/manager.h>
  47 #include <output/table.h>
  48 #include <data/format.h>
  49
  50 #include "minmax.h"
  51 #include "xalloc.h"
  52 #include "xmemdup0.h"
  53
  54 #include "gettext.h"
  55 #define _(msgid) gettext (msgid)
  56
  57 /* (headers) */
  58
  59 /* (specification)
  60    "T-TEST" (tts_):
  61    +groups=custom;
  62    testval=double;
  63    +variables=varlist("PV_NO_SCRATCH | PV_NUMERIC");
  64    +pairs=custom;
  65    missing=miss:!analysis/listwise,
  66    incl:include/!exclude;
  67    +format=fmt:!labels/nolabels;
  68    criteria=:cin(d:criteria,"%s > 0. && %s < 1.").
  69 */
  70 /* (declarations) */
  71 /* (functions) */
  72
  73 enum comparison
  74   {
  75     CMP_LE,
  76     CMP_EQ,
  77   };
  78
  79 /* A pair of variables to be compared. */
  80 struct pair
  81   {
  82     const struct variable *v[2]; /* The paired variables. */
  83     double n;             /* The number of valid variable pairs */
  84     double sum[2];        /* The sum of the members */
  85     double ssq[2];        /* sum of squares of the members */
  86     double std_dev[2];    /* Std deviation of the members */
  87     double s_std_dev[2];  /* Sample Std deviation of the members */
  88     double mean[2];       /* The means of the members */
  89     double correlation;   /* Correlation coefficient between the variables. */
  90     double sum_of_diffs;  /* The sum of the differences */
  91     double sum_of_prod;   /* The sum of the products */
  92     double mean_diff;     /* The mean of the differences */
  93     double ssq_diffs;     /* The sum of the squares of the differences */
  94     double std_dev_diff;  /* The std deviation of the differences */
  95   };
  96
  97 /* Which mode was T-TEST invoked */
  98 enum t_test_mode {
  99   T_1_SAMPLE,                   /* One-sample tests. */
 100   T_IND_SAMPLES,                /* Independent-sample tests. */
 101   T_PAIRED                      /* Paired-sample tests. */
 102 };
 103
 104 /* Total state of a T-TEST procedure. */
 105 struct t_test_proc
 106   {
 107     enum t_test_mode mode;      /* Mode that T-TEST was invoked in. */
 108     double criteria;            /* Confidence interval in (0, 1). */
 109     enum mv_class exclude;      /* Classes of missing values to exclude. */
 110     bool listwise_missing;      /* Drop whole case if one missing var? */
 111     struct fmt_spec weight_format; /* Format of weight variable. */
 112
 113     /* Dependent variables. */
 114     const struct variable **vars;
 115     size_t n_vars;
 116
 117     /* For mode == T_1_SAMPLE. */
 118     double testval;
 119
 120     /* For mode == T_PAIRED only. */
 121     struct pair *pairs;
 122     size_t n_pairs;
 123
 124     /* For mode == T_IND_SAMPLES only. */
 125     struct variable *indep_var; /* Independent variable. */
 126     enum comparison criterion;  /* Type of comparison. */
 127     double critical_value;      /* CMP_LE only: Grouping threshold value. */
 128     union value g_value[2];     /* CMP_EQ only: Per-group indep var values. */
 129   };
 130
 131 /* Statistics Summary Box */
 132 struct ssbox
 133   {
 134     struct tab_table *t;
 135     void (*populate) (struct ssbox *, struct t_test_proc *);
 136     void (*finalize) (struct ssbox *);
 137   };
 138
 139 static void ssbox_create (struct ssbox *, struct t_test_proc *);
 140 static void ssbox_populate (struct ssbox *, struct t_test_proc *);
 141 static void ssbox_finalize (struct ssbox *);
 142
 143 /* Paired Samples Correlation box */
 144 static void pscbox (struct t_test_proc *);
 145
 146
 147 /* Test Results Box. */
 148 struct trbox {
 149   struct tab_table *t;
 150   void (*populate) (struct trbox *, struct t_test_proc *);
 151   void (*finalize) (struct trbox *);
 152   };
 153
 154 static void trbox_create (struct trbox *, struct t_test_proc *);
 155 static void trbox_populate (struct trbox *, struct t_test_proc *);
 156 static void trbox_finalize (struct trbox *);
 157
 158 static void calculate (struct t_test_proc *, struct casereader *,
 159                        const struct dataset *);
 160
 161 static int compare_group_binary (const struct group_statistics *a,
 162                                  const struct group_statistics *b,
 163                                  const struct t_test_proc *);
 164 static unsigned hash_group_binary (const struct group_statistics *g,
 165                                    const struct t_test_proc *p);
 166
 167 int
 168 cmd_t_test (struct lexer *lexer, struct dataset *ds)
 169 {
 170   struct cmd_t_test cmd;
 171   struct t_test_proc proc;
 172   struct casegrouper *grouper;
 173   struct casereader *group;
 174   struct variable *wv;
 175   bool ok = false;
 176
 177   proc.pairs = NULL;
 178   proc.n_pairs = 0;
 179   proc.vars = NULL;
 180   proc.indep_var = NULL;
 181   if (!parse_t_test (lexer, ds, &cmd, &proc))
 182     goto parse_failed;
 183
 184   wv = dict_get_weight (dataset_dict (ds));
 185   proc.weight_format = wv ? *var_get_print_format (wv) : F_8_0;
 186
 187   if ((cmd.sbc_testval != 0) + (cmd.sbc_groups != 0) + (cmd.sbc_pairs != 0)
 188       != 1)
 189     {
 190       msg (SE, _("Exactly one of TESTVAL, GROUPS and PAIRS subcommands "
 191                  "must be specified."));
 192       goto done;
 193     }
 194
 195   proc.mode = (cmd.sbc_testval ? T_1_SAMPLE
 196                : cmd.sbc_groups ? T_IND_SAMPLES
 197                : T_PAIRED);
 198   proc.criteria = cmd.sbc_criteria ? cmd.criteria : 0.95;
 199   proc.exclude = cmd.incl != TTS_INCLUDE ? MV_ANY : MV_SYSTEM;
 200   proc.listwise_missing = cmd.miss == TTS_LISTWISE;
 201
 202   if (proc.mode == T_1_SAMPLE)
 203     proc.testval = cmd.n_testval[0];
 204
 205   if (proc.mode == T_PAIRED)
 206     {
 207       size_t i, j;
 208
 209       if (cmd.sbc_variables)
 210         {
 211           msg (SE, _("VARIABLES subcommand may not be used with PAIRS."));
 212           goto done;
 213         }
 214
 215       /* Fill proc.vars with the unique variables from pairs. */
 216       proc.n_vars = proc.n_pairs * 2;
 217       proc.vars = xmalloc (sizeof *proc.vars * proc.n_vars);
 218       for (i = j = 0; i < proc.n_pairs; i++)
 219         {
 220           proc.vars[j++] = proc.pairs[i].v[0];
 221           proc.vars[j++] = proc.pairs[i].v[1];
 222         }
 223       proc.n_vars = sort_unique (proc.vars, proc.n_vars, sizeof *proc.vars,
 224                                  compare_var_ptrs_by_name, NULL);
 225     }
 226   else
 227     {
 228       if (!cmd.n_variables)
 229         {
 230           msg (SE, _("One or more VARIABLES must be specified."));
 231           goto done;
 232         }
 233       proc.n_vars = cmd.n_variables;
 234       proc.vars = cmd.v_variables;
 235       cmd.v_variables = NULL;
 236     }
 237
 238   /* Data pass. */
 239   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 240   while (casegrouper_get_next_group (grouper, &group))
 241     calculate (&proc, group, ds);
 242   ok = casegrouper_destroy (grouper);
 243   ok = proc_commit (ds) && ok;
 244
 245   if (proc.mode == T_IND_SAMPLES)
 246     {
 247       int v;
 248       /* Destroy any group statistics we created */
 249       for (v = 0; v < proc.n_vars; v++)
 250         {
 251           struct group_proc *grpp = group_proc_get (proc.vars[v]);
 252           hsh_destroy (grpp->group_hash);
 253         }
 254     }
 255
 256 done:
 257   free_t_test (&cmd);
 258 parse_failed:
 259   if (proc.indep_var != NULL)
 260     {
 261       int width = var_get_width (proc.indep_var);
 262       value_destroy (&proc.g_value[0], width);
 263       value_destroy (&proc.g_value[1], width);
 264     }
 265   free (proc.vars);
 266   free (proc.pairs);
 267   return ok ? CMD_SUCCESS : CMD_FAILURE;
 268 }
 269
 270 static int
 271 tts_custom_groups (struct lexer *lexer, struct dataset *ds,
 272                    struct cmd_t_test *cmd UNUSED, void *proc_)
 273 {
 274   struct t_test_proc *proc = proc_;
 275   int n_values;
 276   int width;
 277
 278   lex_match (lexer, '=');
 279
 280   proc->indep_var = parse_variable (lexer, dataset_dict (ds));
 281   if (proc->indep_var == NULL)
 282     {
 283       lex_error (lexer, "expecting variable name in GROUPS subcommand");
 284       return 0;
 285     }
 286   width = var_get_width (proc->indep_var);
 287   value_init (&proc->g_value[0], width);
 288   value_init (&proc->g_value[1], width);
 289
 290   if (!lex_match (lexer, '('))
 291     n_values = 0;
 292   else
 293     {
 294       if (!parse_value (lexer, &proc->g_value[0], width))
 295         return 0;
 296       lex_match (lexer, ',');
 297       if (lex_match (lexer, ')'))
 298         n_values = 1;
 299       else
 300         {
 301           if (!parse_value (lexer, &proc->g_value[1], width)
 302               || !lex_force_match (lexer, ')'))
 303             return 0;
 304           n_values = 2;
 305         }
 306     }
 307
 308   if (var_is_numeric (proc->indep_var))
 309     {
 310       proc->criterion = n_values == 1 ? CMP_LE : CMP_EQ;
 311       if (n_values == 1)
 312         proc->critical_value = proc->g_value[0].f;
 313       else if (n_values == 0)
 314         {
 315           proc->g_value[0].f = 1;
 316           proc->g_value[1].f = 2;
 317         }
 318     }
 319   else
 320     {
 321       proc->criterion = CMP_EQ;
 322       if (n_values != 2)
 323         {
 324           msg (SE, _("When applying GROUPS to a string variable, two "
 325                      "values must be specified."));
 326           return 0;
 327         }
 328     }
 329   return 1;
 330 }
 331
 332 static void
 333 add_pair (struct t_test_proc *proc,
 334           const struct variable *v0, const struct variable *v1)
 335 {
 336   struct pair *p = &proc->pairs[proc->n_pairs++];
 337   p->v[0] = v0;
 338   p->v[1] = v1;
 339 }
 340
 341 static int
 342 tts_custom_pairs (struct lexer *lexer, struct dataset *ds,
 343                   struct cmd_t_test *cmd UNUSED, void *proc_)
 344 {
 345   struct t_test_proc *proc = proc_;
 346
 347   const struct variable **vars1 = NULL;
 348   size_t n_vars1 = 0;
 349
 350   const struct variable **vars2 = NULL;
 351   size_t n_vars2 = 0;
 352
 353   bool paired = false;
 354
 355   size_t n_total_pairs;
 356   size_t i, j;
 357
 358   lex_match (lexer, '=');
 359
 360   if (!parse_variables_const (lexer, dataset_dict (ds), &vars1, &n_vars1,
 361                               PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH))
 362     return 0;
 363
 364   if (lex_match (lexer, T_WITH))
 365     {
 366       if (!parse_variables_const (lexer, dataset_dict (ds), &vars2, &n_vars2,
 367                                   PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH))
 368         {
 369           free (vars1);
 370           return 0;
 371         }
 372
 373       if (lex_match (lexer, '(')
 374           && lex_match_id (lexer, "PAIRED")
 375           && lex_match (lexer, ')'))
 376         {
 377           paired = true;
 378           if (n_vars1 != n_vars2)
 379             {
 380               msg (SE, _("PAIRED was specified but the number of variables "
 381                          "preceding WITH (%zu) did not match the number "
 382                          "following (%zu)."),
 383                    n_vars1, n_vars2);
 384               free (vars1);
 385               free (vars2);
 386               return 0;
 387             }
 388         }
 389     }
 390   else
 391     {
 392       if (n_vars1 < 2)
 393         {
 394           free (vars1);
 395           msg (SE, _("At least two variables must be specified on PAIRS."));
 396           return 0;
 397         }
 398     }
 399
 400   /* Allocate storage for the new pairs. */
 401   n_total_pairs = proc->n_pairs + (paired ? n_vars1
 402                                    : n_vars2 > 0 ? n_vars1 * n_vars2
 403                                    : n_vars1 * (n_vars1 - 1) / 2);
 404   proc->pairs = xnrealloc (proc->pairs, n_total_pairs, sizeof *proc->pairs);
 405
 406   /* Populate the pairs with the appropriate variables. */
 407   if (paired)
 408     for (i = 0; i < n_vars1; i++)
 409       add_pair (proc, vars1[i], vars2[i]);
 410   else if (n_vars2 > 0)
 411     for (i = 0; i < n_vars1; i++)
 412       for (j = 0; j < n_vars2; j++)
 413         add_pair (proc, vars1[i], vars2[j]);
 414   else
 415     for (i = 0; i < n_vars1; i++)
 416       for (j = i + 1; j < n_vars1; j++)
 417         add_pair (proc, vars1[i], vars1[j]);
 418   assert (proc->n_pairs == n_total_pairs);
 419
 420   free (vars1);
 421   free (vars2);
 422   return 1;
 423 }
 424 \f
 425 /* Implementation of the SSBOX object. */
 426
 427 static void ssbox_base_init (struct ssbox *, int cols, int rows);
 428 static void ssbox_base_finalize (struct ssbox *);
 429 static void ssbox_one_sample_init (struct ssbox *, struct t_test_proc *);
 430 static void ssbox_independent_samples_init (struct ssbox *, struct t_test_proc *);
 431 static void ssbox_paired_init (struct ssbox *, struct t_test_proc *);
 432
 433 /* Factory to create an ssbox. */
 434 static void
 435 ssbox_create (struct ssbox *ssb, struct t_test_proc *proc)
 436 {
 437   switch (proc->mode)
 438     {
 439     case T_1_SAMPLE:
 440       ssbox_one_sample_init (ssb, proc);
 441       break;
 442     case T_IND_SAMPLES:
 443       ssbox_independent_samples_init (ssb, proc);
 444       break;
 445     case T_PAIRED:
 446       ssbox_paired_init (ssb, proc);
 447       break;
 448     default:
 449       NOT_REACHED ();
 450     }
 451 }
 452
 453 /* Despatcher for the populate method */
 454 static void
 455 ssbox_populate (struct ssbox *ssb, struct t_test_proc *proc)
 456 {
 457   ssb->populate (ssb, proc);
 458 }
 459
 460 /* Despatcher for finalize */
 461 static void
 462 ssbox_finalize (struct ssbox *ssb)
 463 {
 464   ssb->finalize (ssb);
 465 }
 466
 467 /* Submit the box and clear up */
 468 static void
 469 ssbox_base_finalize (struct ssbox *ssb)
 470 {
 471   tab_submit (ssb->t);
 472 }
 473
 474 /* Initialize a ssbox struct */
 475 static void
 476 ssbox_base_init (struct ssbox *this, int cols, int rows)
 477 {
 478   this->finalize = ssbox_base_finalize;
 479   this->t = tab_create (cols, rows, 0);
 480
 481   tab_columns (this->t, SOM_COL_DOWN, 1);
 482   tab_headers (this->t, 0, 0, 1, 0);
 483   tab_box (this->t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols - 1, rows - 1);
 484   tab_hline (this->t, TAL_2, 0, cols- 1, 1);
 485   tab_dim (this->t, tab_natural_dimensions, NULL);
 486 }
 487 \f
 488 /* ssbox implementations. */
 489
 490 static void ssbox_one_sample_populate (struct ssbox *, struct t_test_proc *);
 491 static void ssbox_independent_samples_populate (struct ssbox *,
 492                                                 struct t_test_proc *);
 493 static void ssbox_paired_populate (struct ssbox *, struct t_test_proc *);
 494
 495 /* Initialize the one_sample ssbox */
 496 static void
 497 ssbox_one_sample_init (struct ssbox *this, struct t_test_proc *proc)
 498 {
 499   const int hsize = 5;
 500   const int vsize = proc->n_vars + 1;
 501
 502   this->populate = ssbox_one_sample_populate;
 503
 504   ssbox_base_init (this, hsize, vsize);
 505   tab_title (this->t, _("One-Sample Statistics"));
 506   tab_vline (this->t, TAL_2, 1, 0, vsize - 1);
 507   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE, _("N"));
 508   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
 509   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 510   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("S.E. Mean"));
 511 }
 512
 513 /* Initialize the independent samples ssbox */
 514 static void
 515 ssbox_independent_samples_init (struct ssbox *this, struct t_test_proc *proc)
 516 {
 517   int hsize=6;
 518   int vsize = proc->n_vars * 2 + 1;
 519
 520   this->populate = ssbox_independent_samples_populate;
 521
 522   ssbox_base_init (this, hsize, vsize);
 523   tab_vline (this->t, TAL_GAP, 1, 0, vsize - 1);
 524   tab_title (this->t, _("Group Statistics"));
 525   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE,
 526             var_get_name (proc->indep_var));
 527   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("N"));
 528   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
 529   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 530   tab_text (this->t, 5, 0, TAB_CENTER | TAT_TITLE, _("S.E. Mean"));
 531 }
 532
 533 /* Populate the ssbox for independent samples */
 534 static void
 535 ssbox_independent_samples_populate (struct ssbox *ssb,
 536                                     struct t_test_proc *proc)
 537 {
 538   int i;
 539
 540   char *val_lab[2];
 541   double indep_value[2];
 542
 543   char prefix[2][3];
 544
 545   for (i = 0; i < 2; i++)
 546     {
 547       union value *value = &proc->g_value[i];
 548       int width = var_get_width (proc->indep_var);
 549
 550       indep_value[i] = (proc->criterion == CMP_LE ? proc->critical_value
 551                         : value->f);
 552
 553       if (val_type_from_width (width) == VAL_NUMERIC)
 554         {
 555           const char *s = var_lookup_value_label (proc->indep_var, value);
 556           val_lab[i] = s ? xstrdup (s) : xasprintf ("%g", indep_value[i]);
 557         }
 558       else
 559         val_lab[i] = xmemdup0 (value_str (value, width), width);
 560     }
 561
 562   if (proc->criterion == CMP_LE)
 563     {
 564       strcpy (prefix[0], ">=");
 565       strcpy (prefix[1], "<");
 566     }
 567   else
 568     {
 569       strcpy (prefix[0], "");
 570       strcpy (prefix[1], "");
 571     }
 572
 573   for (i = 0; i < proc->n_vars; i++)
 574     {
 575       const struct variable *var = proc->vars[i];
 576       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
 577       int count=0;
 578
 579       tab_text (ssb->t, 0, i * 2 + 1, TAB_LEFT,
 580                 var_get_name (proc->vars[i]));
 581       tab_text_format (ssb->t, 1, i * 2 + 1, TAB_LEFT,
 582                        "%s%s", prefix[0], val_lab[0]);
 583       tab_text_format (ssb->t, 1, i * 2 + 1+ 1, TAB_LEFT,
 584                        "%s%s", prefix[1], val_lab[1]);
 585
 586       /* Fill in the group statistics */
 587       for (count = 0; count < 2; count++)
 588         {
 589           union value search_val;
 590           struct group_statistics *gs;
 591
 592           if (proc->criterion == CMP_LE)
 593             search_val.f = proc->critical_value + (count == 0 ? 1.0 : -1.0);
 594           else
 595             search_val = proc->g_value[count];
 596
 597           gs = hsh_find (grp_hash, &search_val);
 598           assert (gs);
 599
 600           tab_double (ssb->t, 2, i * 2 + count+ 1, TAB_RIGHT, gs->n,
 601                       &proc->weight_format);
 602           tab_double (ssb->t, 3, i * 2 + count+ 1, TAB_RIGHT, gs->mean, NULL);
 603           tab_double (ssb->t, 4, i * 2 + count+ 1, TAB_RIGHT, gs->std_dev,
 604                       NULL);
 605           tab_double (ssb->t, 5, i * 2 + count+ 1, TAB_RIGHT, gs->se_mean,
 606                       NULL);
 607         }
 608     }
 609   free (val_lab[0]);
 610   free (val_lab[1]);
 611 }
 612
 613 /* Initialize the paired values ssbox */
 614 static void
 615 ssbox_paired_init (struct ssbox *this, struct t_test_proc *proc)
 616 {
 617   int hsize = 6;
 618   int vsize = proc->n_pairs * 2 + 1;
 619
 620   this->populate = ssbox_paired_populate;
 621
 622   ssbox_base_init (this, hsize, vsize);
 623   tab_title (this->t, _("Paired Sample Statistics"));
 624   tab_vline (this->t, TAL_GAP, 1, 0, vsize - 1);
 625   tab_vline (this->t, TAL_2, 2, 0, vsize - 1);
 626   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
 627   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("N"));
 628   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 629   tab_text (this->t, 5, 0, TAB_CENTER | TAT_TITLE, _("S.E. Mean"));
 630 }
 631
 632 /* Populate the ssbox for paired values */
 633 static void
 634 ssbox_paired_populate (struct ssbox *ssb, struct t_test_proc *proc)
 635 {
 636   int i;
 637
 638   for (i = 0; i < proc->n_pairs; i++)
 639     {
 640       struct pair *p = &proc->pairs[i];
 641       int j;
 642
 643       tab_text_format (ssb->t, 0, i * 2 + 1, TAB_LEFT, _("Pair %d"), i);
 644       for (j=0; j < 2; j++)
 645         {
 646           /* Titles */
 647           tab_text (ssb->t, 1, i * 2 + j + 1, TAB_LEFT,
 648                     var_get_name (p->v[j]));
 649
 650           /* Values */
 651           tab_double (ssb->t, 2, i * 2 + j + 1, TAB_RIGHT, p->mean[j], NULL);
 652           tab_double (ssb->t, 3, i * 2 + j + 1, TAB_RIGHT, p->n,
 653                       &proc->weight_format);
 654           tab_double (ssb->t, 4, i * 2 + j + 1, TAB_RIGHT, p->std_dev[j],
 655                       NULL);
 656           tab_double (ssb->t, 5, i * 2 + j + 1, TAB_RIGHT,
 657                      p->std_dev[j] /sqrt (p->n), NULL);
 658         }
 659     }
 660 }
 661
 662 /* Populate the one sample ssbox */
 663 static void
 664 ssbox_one_sample_populate (struct ssbox *ssb, struct t_test_proc *proc)
 665 {
 666   int i;
 667
 668   for (i = 0; i < proc->n_vars; i++)
 669     {
 670       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
 671
 672       tab_text (ssb->t, 0, i + 1, TAB_LEFT, var_get_name (proc->vars[i]));
 673       tab_double (ssb->t, 1, i + 1, TAB_RIGHT, gs->n, &proc->weight_format);
 674       tab_double (ssb->t, 2, i + 1, TAB_RIGHT, gs->mean, NULL);
 675       tab_double (ssb->t, 3, i + 1, TAB_RIGHT, gs->std_dev, NULL);
 676       tab_double (ssb->t, 4, i + 1, TAB_RIGHT, gs->se_mean, NULL);
 677     }
 678 }
 679 \f
 680 /* Implementation of the Test Results box struct */
 681
 682 static void trbox_base_init (struct trbox *, size_t n_vars, int cols);
 683 static void trbox_base_finalize (struct trbox *);
 684 static void trbox_independent_samples_init (struct trbox *,
 685                                             struct t_test_proc *);
 686 static void trbox_independent_samples_populate (struct trbox *,
 687                                                 struct t_test_proc *);
 688 static void trbox_one_sample_init (struct trbox *, struct t_test_proc *);
 689 static void trbox_one_sample_populate (struct trbox *, struct t_test_proc *);
 690 static void trbox_paired_init (struct trbox *, struct t_test_proc *);
 691 static void trbox_paired_populate (struct trbox *, struct t_test_proc *);
 692
 693 /* Create a trbox according to mode*/
 694 static void
 695 trbox_create (struct trbox *trb, struct t_test_proc *proc)
 696 {
 697   switch (proc->mode)
 698     {
 699     case T_1_SAMPLE:
 700       trbox_one_sample_init (trb, proc);
 701       break;
 702     case T_IND_SAMPLES:
 703       trbox_independent_samples_init (trb, proc);
 704       break;
 705     case T_PAIRED:
 706       trbox_paired_init (trb, proc);
 707       break;
 708     default:
 709       NOT_REACHED ();
 710     }
 711 }
 712
 713 /* Populate a trbox according to proc */
 714 static void
 715 trbox_populate (struct trbox *trb, struct t_test_proc *proc)
 716 {
 717   trb->populate (trb, proc);
 718 }
 719
 720 /* Submit and destroy a trbox */
 721 static void
 722 trbox_finalize (struct trbox *trb)
 723 {
 724   trb->finalize (trb);
 725 }
 726
 727 /* Initialize the independent samples trbox */
 728 static void
 729 trbox_independent_samples_init (struct trbox *self,
 730                                 struct t_test_proc *proc)
 731 {
 732   const int hsize = 11;
 733   const int vsize = proc->n_vars * 2 + 3;
 734
 735   assert (self);
 736   self->populate = trbox_independent_samples_populate;
 737
 738   trbox_base_init (self, proc->n_vars * 2, hsize);
 739   tab_title (self->t, _("Independent Samples Test"));
 740   tab_hline (self->t, TAL_1, 2, hsize - 1, 1);
 741   tab_vline (self->t, TAL_2, 2, 0, vsize - 1);
 742   tab_vline (self->t, TAL_1, 4, 0, vsize - 1);
 743   tab_box (self->t, -1, -1, -1, TAL_1, 2, 1, hsize - 2, vsize - 1);
 744   tab_hline (self->t, TAL_1, hsize - 2, hsize - 1, 2);
 745   tab_box (self->t, -1, -1, -1, TAL_1, hsize - 2, 2, hsize - 1, vsize - 1);
 746   tab_joint_text (self->t, 2, 0, 3, 0,
 747                   TAB_CENTER, _("Levene's Test for Equality of Variances"));
 748   tab_joint_text (self->t, 4, 0, hsize- 1, 0,
 749                   TAB_CENTER, _("t-test for Equality of Means"));
 750
 751   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("F"));
 752   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Sig."));
 753   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("t"));
 754   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("df"));
 755   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
 756   tab_text (self->t, 7, 2, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
 757   tab_text (self->t, 8, 2, TAB_CENTER | TAT_TITLE, _("Std. Error Difference"));
 758   tab_text (self->t, 9, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
 759   tab_text (self->t, 10, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
 760
 761   tab_joint_text_format (self->t, 9, 1, 10, 1, TAB_CENTER,
 762                          _("%g%% Confidence Interval of the Difference"),
 763                          proc->criteria * 100.0);
 764 }
 765
 766 /* Populate the independent samples trbox */
 767 static void
 768 trbox_independent_samples_populate (struct trbox *self,
 769                                     struct t_test_proc *proc)
 770 {
 771   int i;
 772
 773   for (i = 0; i < proc->n_vars; i++)
 774     {
 775       double p, q;
 776
 777       double t;
 778       double df;
 779
 780       double df1, df2;
 781
 782       double pooled_variance;
 783       double std_err_diff;
 784       double mean_diff;
 785
 786       double se2;
 787
 788       const struct variable *var = proc->vars[i];
 789       struct group_proc *grp_data = group_proc_get (var);
 790
 791       struct hsh_table *grp_hash = grp_data->group_hash;
 792
 793       struct group_statistics *gs0;
 794       struct group_statistics *gs1;
 795
 796       union value search_val;
 797
 798       if (proc->criterion == CMP_LE)
 799         search_val.f = proc->critical_value - 1.0;
 800       else
 801         search_val = proc->g_value[0];
 802
 803       gs0 = hsh_find (grp_hash, &search_val);
 804       assert (gs0);
 805
 806       if (proc->criterion == CMP_LE)
 807         search_val.f = proc->critical_value + 1.0;
 808       else
 809         search_val = proc->g_value[1];
 810
 811       gs1 = hsh_find (grp_hash, &search_val);
 812       assert (gs1);
 813
 814
 815       tab_text (self->t, 0, i * 2 + 3, TAB_LEFT, var_get_name (proc->vars[i]));
 816       tab_text (self->t, 1, i * 2 + 3, TAB_LEFT, _("Equal variances assumed"));
 817       tab_double (self->t, 2, i * 2 + 3, TAB_CENTER, grp_data->levene, NULL);
 818
 819       /* Now work out the significance of the Levene test */
 820       df1 = 1;
 821       df2 = grp_data->ugs.n - 2;
 822       q = gsl_cdf_fdist_Q (grp_data->levene, df1, df2);
 823       tab_double (self->t, 3, i * 2 + 3, TAB_CENTER, q, NULL);
 824
 825       df = gs0->n + gs1->n - 2.0;
 826       tab_double (self->t, 5, i * 2 + 3, TAB_RIGHT, df, NULL);
 827
 828       pooled_variance = (gs0->n * pow2 (gs0->s_std_dev)
 829                          + gs1->n *pow2 (gs1->s_std_dev)) / df ;
 830
 831       t = (gs0->mean - gs1->mean) / sqrt (pooled_variance);
 832       t /= sqrt ((gs0->n + gs1->n) / (gs0->n * gs1->n));
 833
 834       tab_double (self->t, 4, i * 2 + 3, TAB_RIGHT, t, NULL);
 835
 836       p = gsl_cdf_tdist_P (t, df);
 837       q = gsl_cdf_tdist_Q (t, df);
 838
 839       tab_double (self->t, 6, i * 2 + 3, TAB_RIGHT, 2.0 * (t > 0 ? q : p),
 840                   NULL);
 841
 842       mean_diff = gs0->mean - gs1->mean;
 843       tab_double (self->t, 7, i * 2 + 3, TAB_RIGHT, mean_diff, NULL);
 844
 845
 846       std_err_diff = sqrt (pow2 (gs0->se_mean) + pow2 (gs1->se_mean));
 847       tab_double (self->t, 8, i * 2 + 3, TAB_RIGHT, std_err_diff, NULL);
 848
 849       /* Now work out the confidence interval */
 850       q = (1 - proc->criteria)/2.0;  /* 2-tailed test */
 851
 852       t = gsl_cdf_tdist_Qinv (q, df);
 853       tab_double (self->t, 9, i * 2 + 3, TAB_RIGHT,
 854                  mean_diff - t * std_err_diff, NULL);
 855
 856       tab_double (self->t, 10, i * 2 + 3, TAB_RIGHT,
 857                  mean_diff + t * std_err_diff, NULL);
 858
 859
 860       /* Now for the \sigma_1 != \sigma_2 case */
 861       tab_text (self->t, 1, i * 2 + 3 + 1,
 862                 TAB_LEFT, _("Equal variances not assumed"));
 863
 864       se2 = ((pow2 (gs0->s_std_dev) / (gs0->n - 1)) +
 865              (pow2 (gs1->s_std_dev) / (gs1->n - 1)));
 866
 867       t = mean_diff / sqrt (se2);
 868       tab_double (self->t, 4, i * 2 + 3 + 1, TAB_RIGHT, t, NULL);
 869
 870       df = pow2 (se2) / ((pow2 (pow2 (gs0->s_std_dev) / (gs0->n - 1))
 871                           / (gs0->n - 1))
 872                          + (pow2 (pow2 (gs1->s_std_dev) / (gs1->n - 1))
 873                             / (gs1->n - 1)));
 874       tab_double (self->t, 5, i * 2 + 3 + 1, TAB_RIGHT, df, NULL);
 875
 876       p = gsl_cdf_tdist_P (t, df);
 877       q = gsl_cdf_tdist_Q (t, df);
 878
 879       tab_double (self->t, 6, i * 2 + 3 + 1, TAB_RIGHT, 2.0 * (t > 0 ? q : p),
 880                   NULL);
 881
 882       /* Now work out the confidence interval */
 883       q = (1 - proc->criteria) / 2.0;  /* 2-tailed test */
 884
 885       t = gsl_cdf_tdist_Qinv (q, df);
 886
 887       tab_double (self->t, 7, i * 2 + 3 + 1, TAB_RIGHT, mean_diff, NULL);
 888       tab_double (self->t, 8, i * 2 + 3 + 1, TAB_RIGHT, std_err_diff, NULL);
 889       tab_double (self->t, 9, i * 2 + 3 + 1, TAB_RIGHT,
 890                  mean_diff - t * std_err_diff, NULL);
 891       tab_double (self->t, 10, i * 2 + 3 + 1, TAB_RIGHT,
 892                  mean_diff + t * std_err_diff, NULL);
 893     }
 894 }
 895
 896 /* Initialize the paired samples trbox */
 897 static void
 898 trbox_paired_init (struct trbox *self, struct t_test_proc *proc)
 899 {
 900   const int hsize=10;
 901   const int vsize=proc->n_pairs+ 3;
 902
 903   self->populate = trbox_paired_populate;
 904
 905   trbox_base_init (self, proc->n_pairs, hsize);
 906   tab_title (self->t, _("Paired Samples Test"));
 907   tab_hline (self->t, TAL_1, 2, 6, 1);
 908   tab_vline (self->t, TAL_2, 2, 0, vsize - 1);
 909   tab_joint_text (self->t, 2, 0, 6, 0, TAB_CENTER, _("Paired Differences"));
 910   tab_box (self->t, -1, -1, -1, TAL_1, 2, 1, 6, vsize - 1);
 911   tab_box (self->t, -1, -1, -1, TAL_1, 6, 0, hsize - 1, vsize - 1);
 912   tab_hline (self->t, TAL_1, 5, 6, 2);
 913   tab_vline (self->t, TAL_GAP, 6, 0, 1);
 914
 915   tab_joint_text_format (self->t, 5, 1, 6, 1, TAB_CENTER,
 916                          _("%g%% Confidence Interval of the Difference"),
 917                          proc->criteria*100.0);
 918
 919   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("Mean"));
 920   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 921   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("Std. Error Mean"));
 922   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
 923   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
 924   tab_text (self->t, 7, 2, TAB_CENTER | TAT_TITLE, _("t"));
 925   tab_text (self->t, 8, 2, TAB_CENTER | TAT_TITLE, _("df"));
 926   tab_text (self->t, 9, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
 927 }
 928
 929 /* Populate the paired samples trbox */
 930 static void
 931 trbox_paired_populate (struct trbox *trb,
 932                        struct t_test_proc *proc)
 933 {
 934   int i;
 935
 936   for (i = 0; i < proc->n_pairs; i++)
 937     {
 938       struct pair *pair = &proc->pairs[i];
 939       double p, q;
 940       double se_mean;
 941
 942       double n = pair->n;
 943       double t;
 944       double df = n - 1;
 945
 946       tab_text_format (trb->t, 0, i + 3, TAB_LEFT, _("Pair %d"), i);
 947       tab_text_format (trb->t, 1, i + 3, TAB_LEFT, "%s - %s",
 948                        var_get_name (pair->v[0]),
 949                        var_get_name (pair->v[1]));
 950       tab_double (trb->t, 2, i + 3, TAB_RIGHT, pair->mean_diff, NULL);
 951       tab_double (trb->t, 3, i + 3, TAB_RIGHT, pair->std_dev_diff, NULL);
 952
 953       /* SE Mean */
 954       se_mean = pair->std_dev_diff / sqrt (n);
 955       tab_double (trb->t, 4, i + 3, TAB_RIGHT, se_mean, NULL);
 956
 957       /* Now work out the confidence interval */
 958       q = (1 - proc->criteria) / 2.0;  /* 2-tailed test */
 959
 960       t = gsl_cdf_tdist_Qinv (q, df);
 961
 962       tab_double (trb->t, 5, i + 3, TAB_RIGHT,
 963                  pair->mean_diff - t * se_mean, NULL);
 964       tab_double (trb->t, 6, i + 3, TAB_RIGHT,
 965                  pair->mean_diff + t * se_mean, NULL);
 966
 967       t = ((pair->mean[0] - pair->mean[1])
 968            / sqrt ((pow2 (pair->s_std_dev[0]) + pow2 (pair->s_std_dev[1])
 969                     - (2 * pair->correlation
 970                        * pair->s_std_dev[0] * pair->s_std_dev[1]))
 971                    / (n - 1)));
 972
 973       tab_double (trb->t, 7, i + 3, TAB_RIGHT, t, NULL);
 974
 975       /* Degrees of freedom */
 976       tab_double (trb->t, 8, i + 3, TAB_RIGHT, df, &proc->weight_format);
 977
 978       p = gsl_cdf_tdist_P (t,df);
 979       q = gsl_cdf_tdist_Q (t,df);
 980
 981       tab_double (trb->t, 9, i + 3, TAB_RIGHT, 2.0 * (t > 0 ? q : p), NULL);
 982     }
 983 }
 984
 985 /* Initialize the one sample trbox */
 986 static void
 987 trbox_one_sample_init (struct trbox *self, struct t_test_proc *proc)
 988 {
 989   const int hsize = 7;
 990   const int vsize = proc->n_vars + 3;
 991
 992   self->populate = trbox_one_sample_populate;
 993
 994   trbox_base_init (self, proc->n_vars, hsize);
 995   tab_title (self->t, _("One-Sample Test"));
 996   tab_hline (self->t, TAL_1, 1, hsize - 1, 1);
 997   tab_vline (self->t, TAL_2, 1, 0, vsize - 1);
 998
 999   tab_joint_text_format (self->t, 1, 0, hsize - 1, 0, TAB_CENTER,
1000                          _("Test Value = %f"), proc->testval);
1001
1002   tab_box (self->t, -1, -1, -1, TAL_1, 1, 1, hsize - 1, vsize - 1);
1003
1004
1005   tab_joint_text_format (self->t, 5, 1, 6, 1, TAB_CENTER,
1006                          _("%g%% Confidence Interval of the Difference"),
1007                          proc->criteria * 100.0);
1008
1009   tab_vline (self->t, TAL_GAP, 6, 1, 1);
1010   tab_hline (self->t, TAL_1, 5, 6, 2);
1011   tab_text (self->t, 1, 2, TAB_CENTER | TAT_TITLE, _("t"));
1012   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("df"));
1013   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
1014   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
1015   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
1016   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
1017 }
1018
1019 /* Populate the one sample trbox */
1020 static void
1021 trbox_one_sample_populate (struct trbox *trb, struct t_test_proc *proc)
1022 {
1023   int i;
1024
1025   assert (trb->t);
1026
1027   for (i = 0; i < proc->n_vars; i++)
1028     {
1029       double t;
1030       double p, q;
1031       double df;
1032       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1033
1034       tab_text (trb->t, 0, i + 3, TAB_LEFT, var_get_name (proc->vars[i]));
1035
1036       t = (gs->mean - proc->testval) * sqrt (gs->n) / gs->std_dev;
1037
1038       tab_double (trb->t, 1, i + 3, TAB_RIGHT, t, NULL);
1039
1040       /* degrees of freedom */
1041       df = gs->n - 1;
1042
1043       tab_double (trb->t, 2, i + 3, TAB_RIGHT, df, &proc->weight_format);
1044
1045       p = gsl_cdf_tdist_P (t, df);
1046       q = gsl_cdf_tdist_Q (t, df);
1047
1048       /* Multiply by 2 to get 2-tailed significance, makeing sure we've got
1049          the correct tail*/
1050       tab_double (trb->t, 3, i + 3, TAB_RIGHT, 2.0 * (t > 0 ? q : p), NULL);
1051       tab_double (trb->t, 4, i + 3, TAB_RIGHT, gs->mean_diff, NULL);
1052
1053
1054       q = (1 - proc->criteria) / 2.0;  /* 2-tailed test */
1055       t = gsl_cdf_tdist_Qinv (q, df);
1056
1057       tab_double (trb->t, 5, i + 3, TAB_RIGHT,
1058                  gs->mean_diff - t * gs->se_mean, NULL);
1059       tab_double (trb->t, 6, i + 3, TAB_RIGHT,
1060                  gs->mean_diff + t * gs->se_mean, NULL);
1061     }
1062 }
1063
1064 /* Base initializer for the generalized trbox */
1065 static void
1066 trbox_base_init (struct trbox *self, size_t data_rows, int cols)
1067 {
1068   const size_t rows = 3 + data_rows;
1069
1070   self->finalize = trbox_base_finalize;
1071   self->t = tab_create (cols, rows, 0);
1072   tab_headers (self->t, 0, 0, 3, 0);
1073   tab_box (self->t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols - 1, rows - 1);
1074   tab_hline (self->t, TAL_2, 0, cols- 1, 3);
1075   tab_dim (self->t, tab_natural_dimensions, NULL);
1076 }
1077
1078 /* Base finalizer for the trbox */
1079 static void
1080 trbox_base_finalize (struct trbox *trb)
1081 {
1082   tab_submit (trb->t);
1083 }
1084
1085 /* Create, populate and submit the Paired Samples Correlation box */
1086 static void
1087 pscbox (struct t_test_proc *proc)
1088 {
1089   const int rows=1+proc->n_pairs;
1090   const int cols=5;
1091   int i;
1092
1093   struct tab_table *table;
1094
1095   table = tab_create (cols, rows, 0);
1096
1097   tab_columns (table, SOM_COL_DOWN, 1);
1098   tab_headers (table, 0, 0, 1, 0);
1099   tab_box (table, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols - 1, rows - 1);
1100   tab_hline (table, TAL_2, 0, cols - 1, 1);
1101   tab_vline (table, TAL_2, 2, 0, rows - 1);
1102   tab_dim (table, tab_natural_dimensions, NULL);
1103   tab_title (table, _("Paired Samples Correlations"));
1104
1105   /* column headings */
1106   tab_text (table, 2, 0, TAB_CENTER | TAT_TITLE, _("N"));
1107   tab_text (table, 3, 0, TAB_CENTER | TAT_TITLE, _("Correlation"));
1108   tab_text (table, 4, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
1109
1110   for (i = 0; i < proc->n_pairs; i++)
1111     {
1112       struct pair *pair = &proc->pairs[i];
1113       double df = pair->n - 2;
1114       double p, q;
1115
1116       /* corr2 will mathematically always be in the range [0, 1.0].  Inaccurate
1117          calculations sometimes cause it to be slightly greater than 1.0, so
1118          force it into the correct range to avoid NaN from sqrt(). */
1119       double corr2 = MIN (1.0, pow2 (pair->correlation));
1120       double correlation_t = pair->correlation * sqrt (df) / sqrt (1 - corr2);
1121
1122       /* row headings */
1123       tab_text_format (table, 0, i + 1, TAB_LEFT | TAT_TITLE,
1124                        _("Pair %d"), i);
1125       tab_text_format (table, 1, i + 1, TAB_LEFT | TAT_TITLE,
1126                        _("%s & %s"),
1127                        var_get_name (pair->v[0]),
1128                        var_get_name (pair->v[1]));
1129
1130       /* row data */
1131       tab_double (table, 2, i + 1, TAB_RIGHT, pair->n, &proc->weight_format);
1132       tab_double (table, 3, i + 1, TAB_RIGHT, pair->correlation, NULL);
1133
1134       p = gsl_cdf_tdist_P (correlation_t, df);
1135       q = gsl_cdf_tdist_Q (correlation_t, df);
1136       tab_double (table, 4, i + 1, TAB_RIGHT,
1137                  2.0 * (correlation_t > 0 ? q : p), NULL);
1138     }
1139
1140   tab_submit (table);
1141 }
1142 \f
1143 /* Calculation Implementation */
1144
1145 /* Calculations common to all variants of the T test. */
1146 static void
1147 common_calc (const struct dictionary *dict,
1148              struct t_test_proc *proc,
1149              struct casereader *reader)
1150 {
1151   struct ccase *c;
1152   int i;
1153
1154   for (i = 0; i < proc->n_vars; i++)
1155     {
1156       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1157       gs->sum = 0;
1158       gs->n = 0;
1159       gs->ssq = 0;
1160       gs->sum_diff = 0;
1161     }
1162
1163   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
1164     {
1165       double weight = dict_get_case_weight (dict, c, NULL);
1166
1167       /* Listwise has to be implicit if the independent variable
1168          is missing ?? */
1169       if (proc->mode == T_IND_SAMPLES)
1170         {
1171           if (var_is_value_missing (proc->indep_var,
1172                                     case_data (c, proc->indep_var),
1173                                     proc->exclude))
1174             continue;
1175         }
1176
1177       for (i = 0; i < proc->n_vars; i++)
1178         {
1179           const struct variable *v = proc->vars[i];
1180           const union value *val = case_data (c, v);
1181
1182           if (!var_is_value_missing (v, val, proc->exclude))
1183             {
1184               struct group_statistics *gs;
1185               gs = &group_proc_get (v)->ugs;
1186
1187               gs->n += weight;
1188               gs->sum += weight * val->f;
1189               gs->ssq += weight * pow2 (val->f);
1190             }
1191         }
1192     }
1193   casereader_destroy (reader);
1194
1195   for (i = 0; i < proc->n_vars; i++)
1196     {
1197       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1198
1199       gs->mean = gs->sum / gs->n;
1200       gs->s_std_dev = sqrt (((gs->ssq / gs->n) - pow2 (gs->mean)));
1201       gs->std_dev = sqrt (gs->n / (gs->n- 1)
1202                           * ((gs->ssq / gs->n) - pow2 (gs->mean)));
1203       gs->se_mean = gs->std_dev / sqrt (gs->n);
1204       gs->mean_diff = gs->sum_diff / gs->n;
1205     }
1206 }
1207
1208 /* Calculations for one sample T test. */
1209 static int
1210 one_sample_calc (const struct dictionary *dict, struct t_test_proc *proc,
1211                  struct casereader *reader)
1212 {
1213   struct ccase *c;
1214   int i;
1215
1216   for (i = 0; i < proc->n_vars; i++)
1217     {
1218       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1219       gs->sum_diff = 0;
1220     }
1221
1222   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
1223     {
1224       double weight = dict_get_case_weight (dict, c, NULL);
1225       for (i = 0; i < proc->n_vars; i++)
1226         {
1227           const struct variable *v = proc->vars[i];
1228           struct group_statistics *gs = &group_proc_get (v)->ugs;
1229           const union value *val = case_data (c, v);
1230           if (!var_is_value_missing (v, val, proc->exclude))
1231             gs->sum_diff += weight * (val->f - proc->testval);
1232         }
1233     }
1234
1235   for (i = 0; i < proc->n_vars; i++)
1236     {
1237       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1238       gs->mean_diff = gs->sum_diff / gs->n;
1239     }
1240
1241   casereader_destroy (reader);
1242
1243   return 0;
1244 }
1245
1246 static int
1247 paired_calc (const struct dictionary *dict, struct t_test_proc *proc,
1248              struct casereader *reader)
1249 {
1250   struct ccase *c;
1251   int i;
1252
1253   for (i = 0; i < proc->n_pairs; i++)
1254     {
1255       struct pair *pair = &proc->pairs[i];
1256       pair->n = 0;
1257       pair->sum[0] = pair->sum[1] = 0;
1258       pair->ssq[0] = pair->ssq[1] = 0;
1259       pair->sum_of_prod = 0;
1260       pair->correlation = 0;
1261       pair->sum_of_diffs = 0;
1262       pair->ssq_diffs = 0;
1263     }
1264
1265   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
1266     {
1267       double weight = dict_get_case_weight (dict, c, NULL);
1268       for (i = 0; i < proc->n_pairs; i++)
1269         {
1270           struct pair *pair = &proc->pairs[i];
1271           const struct variable *v0 = pair->v[0];
1272           const struct variable *v1 = pair->v[1];
1273
1274           const union value *val0 = case_data (c, v0);
1275           const union value *val1 = case_data (c, v1);
1276
1277           if (!var_is_value_missing (v0, val0, proc->exclude)
1278               && !var_is_value_missing (v1, val1, proc->exclude))
1279             {
1280               pair->n += weight;
1281               pair->sum[0] += weight * val0->f;
1282               pair->sum[1] += weight * val1->f;
1283               pair->ssq[0] += weight * pow2 (val0->f);
1284               pair->ssq[1] += weight * pow2 (val1->f);
1285               pair->sum_of_prod += weight * val0->f * val1->f;
1286               pair->sum_of_diffs += weight * (val0->f - val1->f);
1287               pair->ssq_diffs += weight * pow2 (val0->f - val1->f);
1288             }
1289         }
1290     }
1291
1292   for (i = 0; i < proc->n_pairs; i++)
1293     {
1294       struct pair *pair = &proc->pairs[i];
1295       const double n = pair->n;
1296       int j;
1297
1298       for (j=0; j < 2; j++)
1299         {
1300           pair->mean[j] = pair->sum[j] / n;
1301           pair->s_std_dev[j] = sqrt ((pair->ssq[j] / n
1302                                       - pow2 (pair->mean[j])));
1303           pair->std_dev[j] = sqrt (n / (n- 1) * (pair->ssq[j] / n
1304                                                 - pow2 (pair->mean[j])));
1305         }
1306
1307       pair->correlation = (pair->sum_of_prod / pair->n
1308                            - pair->mean[0] * pair->mean[1]);
1309       /* correlation now actually contains the covariance */
1310       pair->correlation /= pair->std_dev[0] * pair->std_dev[1];
1311       pair->correlation *= pair->n / (pair->n - 1);
1312
1313       pair->mean_diff = pair->sum_of_diffs / n;
1314       pair->std_dev_diff = sqrt (n / (n - 1) * ((pair->ssq_diffs / n)
1315                                                 - pow2 (pair->mean_diff)));
1316     }
1317
1318   casereader_destroy (reader);
1319   return 0;
1320 }
1321
1322 static int
1323 group_calc (const struct dictionary *dict, struct t_test_proc *proc,
1324             struct casereader *reader)
1325 {
1326   struct ccase *c;
1327   int i;
1328
1329   for (i = 0; i < proc->n_vars; i++)
1330     {
1331       struct group_proc *ttpr = group_proc_get (proc->vars[i]);
1332       int j;
1333
1334       /* There's always 2 groups for a T - TEST */
1335       ttpr->n_groups = 2;
1336       ttpr->group_hash = hsh_create (2,
1337                                      (hsh_compare_func *) compare_group_binary,
1338                                      (hsh_hash_func *) hash_group_binary,
1339                                      (hsh_free_func *) free_group,
1340                                      proc);
1341
1342       for (j = 0; j < 2; j++)
1343         {
1344           struct group_statistics *gs = xmalloc (sizeof *gs);
1345           gs->sum = 0;
1346           gs->n = 0;
1347           gs->ssq = 0;
1348           if (proc->criterion == CMP_EQ)
1349             gs->id = proc->g_value[j];
1350           else
1351             {
1352               if (j == 0)
1353                 gs->id.f = proc->critical_value - 1.0;
1354               else
1355                 gs->id.f = proc->critical_value + 1.0;
1356             }
1357
1358           hsh_insert (ttpr->group_hash, gs);
1359         }
1360     }
1361
1362   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
1363     {
1364       const double weight = dict_get_case_weight (dict, c, NULL);
1365       const union value *gv;
1366
1367       if (var_is_value_missing (proc->indep_var,
1368                                 case_data (c, proc->indep_var), proc->exclude))
1369         continue;
1370
1371       gv = case_data (c, proc->indep_var);
1372       for (i = 0; i < proc->n_vars; i++)
1373         {
1374           const struct variable *var = proc->vars[i];
1375           const union value *val = case_data (c, var);
1376           struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1377           struct group_statistics *gs = hsh_find (grp_hash, gv);
1378
1379           /* If the independent variable doesn't match either of the values
1380              for this case then move on to the next case. */
1381           if (gs == NULL)
1382             break;
1383
1384           if (!var_is_value_missing (var, val, proc->exclude))
1385             {
1386               gs->n += weight;
1387               gs->sum += weight * val->f;
1388               gs->ssq += weight * pow2 (val->f);
1389             }
1390         }
1391     }
1392
1393   for (i = 0; i < proc->n_vars; i++)
1394     {
1395       const struct variable *var = proc->vars[i];
1396       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1397       struct hsh_iterator g;
1398       struct group_statistics *gs;
1399       int count = 0;
1400
1401       for (gs = hsh_first (grp_hash, &g); gs != NULL;
1402            gs = hsh_next (grp_hash, &g))
1403         {
1404           gs->mean = gs->sum / gs->n;
1405           gs->s_std_dev = sqrt (((gs->ssq / gs->n) - pow2 (gs->mean)));
1406           gs->std_dev = sqrt (gs->n / (gs->n- 1)
1407                               * ((gs->ssq / gs->n) - pow2 (gs->mean)));
1408           gs->se_mean = gs->std_dev / sqrt (gs->n);
1409           count++;
1410         }
1411       assert (count == 2);
1412     }
1413
1414   casereader_destroy (reader);
1415
1416   return 0;
1417 }
1418
1419 static void
1420 calculate (struct t_test_proc *proc,
1421            struct casereader *input, const struct dataset *ds)
1422 {
1423   const struct dictionary *dict = dataset_dict (ds);
1424   struct ssbox stat_summary_box;
1425   struct trbox test_results_box;
1426   struct taint *taint;
1427   struct ccase *c;
1428
1429   c = casereader_peek (input, 0);
1430   if (c == NULL)
1431     {
1432       casereader_destroy (input);
1433       return;
1434     }
1435   output_split_file_values (ds, c);
1436   case_unref (c);
1437
1438   if (proc->listwise_missing)
1439     input = casereader_create_filter_missing (input,
1440                                               proc->vars,
1441                                               proc->n_vars,
1442                                               proc->exclude, NULL, NULL);
1443   input = casereader_create_filter_weight (input, dict, NULL, NULL);
1444   taint = taint_clone (casereader_get_taint (input));
1445
1446   common_calc (dict, proc, casereader_clone (input));
1447   switch (proc->mode)
1448     {
1449     case T_1_SAMPLE:
1450       one_sample_calc (dict, proc, input);
1451       break;
1452     case T_PAIRED:
1453       paired_calc (dict, proc, input);
1454       break;
1455     case T_IND_SAMPLES:
1456       group_calc (dict, proc, casereader_clone (input));
1457       levene (dict, input, proc->indep_var, proc->n_vars, proc->vars,
1458               proc->exclude);
1459       break;
1460     default:
1461       NOT_REACHED ();
1462     }
1463
1464   if (!taint_has_tainted_successor (taint))
1465     {
1466       ssbox_create (&stat_summary_box, proc);
1467       ssbox_populate (&stat_summary_box, proc);
1468       ssbox_finalize (&stat_summary_box);
1469
1470       if (proc->mode == T_PAIRED)
1471         pscbox (proc);
1472
1473       trbox_create (&test_results_box, proc);
1474       trbox_populate (&test_results_box, proc);
1475       trbox_finalize (&test_results_box);
1476     }
1477
1478   taint_destroy (taint);
1479 }
1480
1481 /* return 0 if G belongs to group 0,
1482           1 if it belongs to group 1,
1483           2 if it belongs to neither group */
1484 static int
1485 which_group (const struct group_statistics *g,
1486              const struct t_test_proc *proc)
1487 {
1488   int width = var_get_width (proc->indep_var);
1489
1490   if (0 == value_compare_3way (&g->id, &proc->g_value[0], width))
1491     return 0;
1492
1493   if (0 == value_compare_3way (&g->id, &proc->g_value[1], width))
1494     return 1;
1495
1496   return 2;
1497 }
1498
1499 /* Return -1 if the id of a is less than b; +1 if greater than and
1500    0 if equal */
1501 static int
1502 compare_group_binary (const struct group_statistics *a,
1503                       const struct group_statistics *b,
1504                       const struct t_test_proc *proc)
1505 {
1506   int flag_a;
1507   int flag_b;
1508
1509   if (proc->criterion == CMP_LE)
1510     {
1511       flag_a = (a->id.f < proc->critical_value);
1512       flag_b = (b->id.f < proc->critical_value);
1513     }
1514   else
1515     {
1516       flag_a = which_group (a, proc);
1517       flag_b = which_group (b, proc);
1518     }
1519
1520   if (flag_a < flag_b)
1521     return - 1;
1522
1523   return (flag_a > flag_b);
1524 }
1525
1526 /* This is a degenerate case of a hash, since it can only return three possible
1527    values.  It's really a comparison, being used as a hash function */
1528
1529 static unsigned
1530 hash_group_binary (const struct group_statistics *g,
1531                    const struct t_test_proc *proc)
1532 {
1533   return (proc->criterion == CMP_LE
1534           ? g->id.f < proc->critical_value
1535           : which_group (g, proc));
1536 }
1537
1538 /*
1539   Local Variables:
1540   mode: c
1541   End:
1542 */