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