src/language/stats/t-test.q

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