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