src/language/stats/t-test.q

   1 /* PSPP - computes sample statistics. -*-c-*-
   2
   3    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   4    Written by John Williams <johnr.williams@stonebow.otago.ac.nz>.
   5    Almost completly re-written by John Darrington 2004
   6
   7    This program is free software; you can redistribute it and/or
   8    modify it under the terms of the GNU General Public License as
   9    published by the Free Software Foundation; either version 2 of the
  10    License, or (at your option) any later version.
  11
  12    This program is distributed in the hope that it will be useful, but
  13    WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15    General Public License for more details.
  16
  17    You should have received a copy of the GNU General Public License
  18    along with this program; if not, write to the Free Software
  19    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  20    02110-1301, USA. */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_cdf.h>
  24 #include "message.h"
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27 #include <math.h>
  28 #include "alloc.h"
  29 #include "case.h"
  30 #include "casefile.h"
  31 #include "command.h"
  32 #include "dictionary.h"
  33 #include "message.h"
  34 #include "group-proc.h"
  35 #include "hash.h"
  36 #include "levene.h"
  37 #include "lexer.h"
  38 #include "magic.h"
  39 #include "misc.h"
  40 #include "size_max.h"
  41 #include "manager.h"
  42 #include "str.h"
  43 #include "table.h"
  44 #include "value-labels.h"
  45 #include "variable.h"
  46 #include "procedure.h"
  47
  48 #include "gettext.h"
  49 #define _(msgid) gettext (msgid)
  50
  51 /* (headers) */
  52
  53 /* (specification)
  54    "T-TEST" (tts_):
  55      +groups=custom;
  56      testval=double;
  57      variables=varlist("PV_NO_SCRATCH | PV_NUMERIC");
  58      pairs=custom;
  59      +missing=miss:!analysis/listwise,
  60              incl:include/!exclude;
  61      format=fmt:!labels/nolabels;
  62      criteria=:cin(d:criteria,"%s > 0. && %s < 1.").
  63 */
  64 /* (declarations) */
  65 /* (functions) */
  66
  67
  68
  69
  70 /* Function to use for testing for missing values */
  71 static is_missing_func *value_is_missing;
  72
  73 /* Variable for the GROUPS subcommand, if given. */
  74 static struct variable *indep_var;
  75
  76 enum comparison
  77   {
  78     CMP_LE = -2,
  79     CMP_EQ = 0,
  80   };
  81
  82 struct group_properties
  83 {
  84   /* The comparison criterion */
  85   enum comparison criterion;
  86
  87   /* The width of the independent variable */
  88   int indep_width ;
  89
  90   union {
  91     /* The value of the independent variable at which groups are determined to
  92        belong to one group or the other */
  93     double critical_value;
  94
  95
  96     /* The values of the independent variable for each group */
  97     union value g_value[2];
  98   } v ;
  99
 100 };
 101
 102
 103 static struct group_properties gp ;
 104
 105
 106
 107 /* PAIRS: Number of pairs to be compared ; each pair. */
 108 static int n_pairs = 0 ;
 109 struct pair
 110 {
 111   /* The variables comprising the pair */
 112   struct variable *v[2];
 113
 114   /* The number of valid variable pairs */
 115   double n;
 116
 117   /* The sum of the members */
 118   double sum[2];
 119
 120   /* sum of squares of the members */
 121   double ssq[2];
 122
 123   /* Std deviation of the members */
 124   double std_dev[2];
 125
 126
 127   /* Sample Std deviation of the members */
 128   double s_std_dev[2];
 129
 130   /* The means of the members */
 131   double mean[2];
 132
 133   /* The correlation coefficient between the variables */
 134   double correlation;
 135
 136   /* The sum of the differences */
 137   double sum_of_diffs;
 138
 139   /* The sum of the products */
 140   double sum_of_prod;
 141
 142   /* The mean of the differences */
 143   double mean_diff;
 144
 145   /* The sum of the squares of the differences */
 146   double ssq_diffs;
 147
 148   /* The std deviation of the differences */
 149   double std_dev_diff;
 150 };
 151
 152 static struct pair *pairs=0;
 153
 154 static int parse_value (union value * v, int type) ;
 155
 156 /* Structures and Functions for the Statistics Summary Box */
 157 struct ssbox;
 158 typedef void populate_ssbox_func(struct ssbox *ssb,
 159                                             struct cmd_t_test *cmd);
 160 typedef void finalize_ssbox_func(struct ssbox *ssb);
 161
 162 struct ssbox
 163 {
 164   struct tab_table *t;
 165
 166   populate_ssbox_func *populate;
 167   finalize_ssbox_func *finalize;
 168
 169 };
 170
 171 /* Create a ssbox */
 172 void ssbox_create(struct ssbox *ssb,   struct cmd_t_test *cmd, int mode);
 173
 174 /* Populate a ssbox according to cmd */
 175 void ssbox_populate(struct ssbox *ssb, struct cmd_t_test *cmd);
 176
 177 /* Submit and destroy a ssbox */
 178 void ssbox_finalize(struct ssbox *ssb);
 179
 180 /* A function to create, populate and submit the Paired Samples Correlation
 181    box */
 182 void pscbox(void);
 183
 184
 185 /* Structures and Functions for the Test Results Box */
 186 struct trbox;
 187
 188 typedef void populate_trbox_func(struct trbox *trb,
 189                                  struct cmd_t_test *cmd);
 190 typedef void finalize_trbox_func(struct trbox *trb);
 191
 192 struct trbox {
 193   struct tab_table *t;
 194   populate_trbox_func *populate;
 195   finalize_trbox_func *finalize;
 196 };
 197
 198 /* Create a trbox */
 199 void trbox_create(struct trbox *trb,   struct cmd_t_test *cmd, int mode);
 200
 201 /* Populate a ssbox according to cmd */
 202 void trbox_populate(struct trbox *trb, struct cmd_t_test *cmd);
 203
 204 /* Submit and destroy a ssbox */
 205 void trbox_finalize(struct trbox *trb);
 206
 207 /* Which mode was T-TEST invoked */
 208 enum {
 209   T_1_SAMPLE = 0 ,
 210   T_IND_SAMPLES,
 211   T_PAIRED
 212 };
 213
 214
 215 static int common_calc (const struct ccase *, void *);
 216 static void common_precalc (struct cmd_t_test *);
 217 static void common_postcalc (struct cmd_t_test *);
 218
 219 static int one_sample_calc (const struct ccase *, void *);
 220 static void one_sample_precalc (struct cmd_t_test *);
 221 static void one_sample_postcalc (struct cmd_t_test *);
 222
 223 static int  paired_calc (const struct ccase *, void *);
 224 static void paired_precalc (struct cmd_t_test *);
 225 static void paired_postcalc (struct cmd_t_test *);
 226
 227 static void group_precalc (struct cmd_t_test *);
 228 static int  group_calc (const struct ccase *, struct cmd_t_test *);
 229 static void group_postcalc (struct cmd_t_test *);
 230
 231
 232 static bool calculate(const struct casefile *cf, void *_mode);
 233
 234 static  int mode;
 235
 236 static struct cmd_t_test cmd;
 237
 238 static int bad_weight_warn;
 239
 240
 241 static int compare_group_binary(const struct group_statistics *a,
 242                                 const struct group_statistics *b,
 243                                 const struct group_properties *p);
 244
 245
 246 static unsigned  hash_group_binary(const struct group_statistics *g,
 247                                    const struct group_properties *p);
 248
 249
 250
 251 int
 252 cmd_t_test(void)
 253 {
 254   bool ok;
 255
 256   if ( !parse_t_test(&cmd) )
 257     return CMD_FAILURE;
 258
 259   if (! cmd.sbc_criteria)
 260     cmd.criteria=0.95;
 261
 262   {
 263     int m=0;
 264     if (cmd.sbc_testval) ++m;
 265     if (cmd.sbc_groups) ++m;
 266     if (cmd.sbc_pairs) ++m;
 267
 268     if ( m != 1)
 269       {
 270         msg(SE,
 271             _("TESTVAL, GROUPS and PAIRS subcommands are mutually exclusive.")
 272             );
 273         free_t_test(&cmd);
 274         return CMD_FAILURE;
 275       }
 276   }
 277
 278   if (cmd.sbc_testval)
 279     mode=T_1_SAMPLE;
 280   else if (cmd.sbc_groups)
 281     mode=T_IND_SAMPLES;
 282   else
 283     mode=T_PAIRED;
 284
 285   if ( mode == T_PAIRED)
 286     {
 287       if (cmd.sbc_variables)
 288         {
 289           msg(SE, _("VARIABLES subcommand is not appropriate with PAIRS"));
 290           free_t_test(&cmd);
 291           return CMD_FAILURE;
 292         }
 293       else
 294         {
 295           /* Iterate through the pairs and put each variable that is a
 296              member of a pair into cmd.v_variables */
 297
 298           int i;
 299           struct hsh_iterator hi;
 300           struct hsh_table *hash;
 301           struct variable *v;
 302
 303           hash = hsh_create (n_pairs, compare_var_names, hash_var_name, 0, 0);
 304
 305           for (i=0; i < n_pairs; ++i)
 306             {
 307               hsh_insert(hash,pairs[i].v[0]);
 308               hsh_insert(hash,pairs[i].v[1]);
 309             }
 310
 311           assert(cmd.n_variables == 0);
 312           cmd.n_variables = hsh_count(hash);
 313
 314           cmd.v_variables = xnrealloc (cmd.v_variables, cmd.n_variables,
 315                                        sizeof *cmd.v_variables);
 316           /* Iterate through the hash */
 317           for (i=0,v = (struct variable *) hsh_first(hash,&hi);
 318                v != 0;
 319                v=hsh_next(hash,&hi) )
 320             cmd.v_variables[i++]=v;
 321
 322           hsh_destroy(hash);
 323         }
 324     }
 325   else if ( !cmd.sbc_variables)
 326     {
 327       msg(SE, _("One or more VARIABLES must be specified."));
 328       free_t_test(&cmd);
 329       return CMD_FAILURE;
 330     }
 331
 332
 333   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
 334   if (cmd.incl == TTS_INCLUDE )
 335     value_is_missing = mv_is_value_system_missing;
 336   else
 337     value_is_missing = mv_is_value_missing;
 338
 339   bad_weight_warn = 1;
 340
 341   ok = multipass_procedure_with_splits (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 cmd_t_test *cmd UNUSED)
 364 {
 365   int n_group_values=0;
 366
 367   lex_match('=');
 368
 369   indep_var = parse_variable ();
 370   if (!indep_var)
 371     {
 372       lex_error ("expecting variable name in GROUPS subcommand");
 373       return 0;
 374     }
 375
 376   if (indep_var->type == T_STRING && indep_var->width > MAX_SHORT_STRING)
 377     {
 378       msg (SE, _("Long string variable %s is not valid here."),
 379            indep_var->name);
 380       return 0;
 381     }
 382
 383   if (!lex_match ('('))
 384     {
 385       if (indep_var->type == NUMERIC)
 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 (&gp.v.g_value[0], indep_var->type))
 405       return 0;
 406
 407   lex_match (',');
 408   if (lex_match (')'))
 409     {
 410       if (indep_var->type != NUMERIC)
 411         {
 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 (&gp.v.g_value[1], indep_var->type))
 425     return 0;
 426
 427   n_group_values = 2;
 428   if (!lex_force_match (')'))
 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 cmd_t_test *cmd 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('=');
 453
 454   n_vars=0;
 455   if (!parse_variables (default_dict, &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 (T_WITH))
 465     {
 466       n_before_WITH = n_vars;
 467       if (!parse_variables (default_dict, &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 ('(') && lex_match_id ("PAIRED") && lex_match (')'));
 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 (union value * v, int type )
 568 {
 569   if (type == NUMERIC)
 570     {
 571       if (!lex_force_num ())
 572         return 0;
 573       v->f = tokval;
 574     }
 575   else
 576     {
 577       if (!lex_force_string ())
 578         return 0;
 579       strncpy (v->s, ds_c_str (&tokstr), ds_length (&tokstr));
 580     }
 581
 582   lex_get ();
 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         assert(0);
 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, 0, _("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_title (this->t, 0, _("Group Statistics"));
 702   tab_vline(this->t,0,1,0,vsize - 1);
 703   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE, indep_var->name);
 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   char *val_lab0=0;
 719   char *val_lab1=0;
 720   double indep_value[2];
 721
 722   char prefix[2][3]={"",""};
 723
 724   if ( indep_var->type == NUMERIC )
 725     {
 726       val_lab0 = val_labs_find( indep_var->val_labs,gp.v.g_value[0]);
 727       val_lab1 = val_labs_find( indep_var->val_labs,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, cmd->v_variables[i]->name);
 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, 0, _("Paired Sample Statistics"));
 826   tab_vline(this->t,TAL_0,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, pairs[i].v[j]->name);
 858
 859           /* Values */
 860           tab_float (ssb->t,2, i*2+j+1, TAB_RIGHT, pairs[i].mean[j], 8, 2);
 861           tab_float (ssb->t,3, i*2+j+1, TAB_RIGHT, pairs[i].n, 2, 0);
 862           tab_float (ssb->t,4, i*2+j+1, TAB_RIGHT, pairs[i].std_dev[j], 8, 3);
 863           tab_float (ssb->t,5, i*2+j+1, TAB_RIGHT, pairs[i].std_dev[j]/sqrt(pairs[i].n), 8, 3);
 864
 865         }
 866     }
 867 }
 868
 869 /* Populate the one sample ssbox */
 870 void
 871 ssbox_one_sample_populate(struct ssbox *ssb, struct cmd_t_test *cmd)
 872 {
 873   int i;
 874
 875   assert(ssb->t);
 876
 877   for (i=0; i < cmd->n_variables; ++i)
 878     {
 879       struct group_statistics *gs = &group_proc_get (cmd->v_variables[i])->ugs;
 880
 881       tab_text (ssb->t, 0, i+1, TAB_LEFT, cmd->v_variables[i]->name);
 882       tab_float (ssb->t,1, i+1, TAB_RIGHT, gs->n, 2, 0);
 883       tab_float (ssb->t,2, i+1, TAB_RIGHT, gs->mean, 8, 2);
 884       tab_float (ssb->t,3, i+1, TAB_RIGHT, gs->std_dev, 8, 2);
 885       tab_float (ssb->t,4, i+1, TAB_RIGHT, gs->se_mean, 8, 3);
 886     }
 887
 888 }
 889
 890
 891
 892 /* Implementation of the Test Results box struct */
 893
 894 void trbox_base_init(struct trbox *self,size_t n_vars, int cols);
 895 void trbox_base_finalize(struct trbox *trb);
 896
 897 void trbox_independent_samples_init(struct trbox *trb,
 898                                     struct cmd_t_test *cmd );
 899
 900 void trbox_independent_samples_populate(struct trbox *trb,
 901                                         struct cmd_t_test *cmd);
 902
 903 void trbox_one_sample_init(struct trbox *self,
 904                       struct cmd_t_test *cmd );
 905
 906 void trbox_one_sample_populate(struct trbox *trb,
 907                                struct cmd_t_test *cmd);
 908
 909 void trbox_paired_init(struct trbox *self,
 910                        struct cmd_t_test *cmd );
 911
 912 void trbox_paired_populate(struct trbox *trb,
 913                       struct cmd_t_test *cmd);
 914
 915
 916
 917 /* Create a trbox according to mode*/
 918 void
 919 trbox_create(struct trbox *trb,
 920              struct cmd_t_test *cmd, int mode)
 921 {
 922     switch (mode)
 923       {
 924       case T_1_SAMPLE:
 925         trbox_one_sample_init(trb,cmd);
 926         break;
 927       case T_IND_SAMPLES:
 928         trbox_independent_samples_init(trb,cmd);
 929         break;
 930       case T_PAIRED:
 931         trbox_paired_init(trb,cmd);
 932         break;
 933       default:
 934         assert(0);
 935       }
 936 }
 937
 938 /* Populate a trbox according to cmd */
 939 void
 940 trbox_populate(struct trbox *trb, struct cmd_t_test *cmd)
 941 {
 942   trb->populate(trb,cmd);
 943 }
 944
 945 /* Submit and destroy a trbox */
 946 void
 947 trbox_finalize(struct trbox *trb)
 948 {
 949   trb->finalize(trb);
 950 }
 951
 952 /* Initialize the independent samples trbox */
 953 void
 954 trbox_independent_samples_init(struct trbox *self,
 955                            struct cmd_t_test *cmd UNUSED)
 956 {
 957   const int hsize=11;
 958   const int vsize=cmd->n_variables*2+3;
 959
 960   assert(self);
 961   self->populate = trbox_independent_samples_populate;
 962
 963   trbox_base_init(self,cmd->n_variables*2,hsize);
 964   tab_title(self->t,0,_("Independent Samples Test"));
 965   tab_hline(self->t,TAL_1,2,hsize-1,1);
 966   tab_vline(self->t,TAL_2,2,0,vsize-1);
 967   tab_vline(self->t,TAL_1,4,0,vsize-1);
 968   tab_box(self->t,-1,-1,-1,TAL_1, 2,1,hsize-2,vsize-1);
 969   tab_hline(self->t,TAL_1, hsize-2,hsize-1,2);
 970   tab_box(self->t,-1,-1,-1,TAL_1, hsize-2,2,hsize-1,vsize-1);
 971   tab_joint_text(self->t, 2, 0, 3, 0,
 972                  TAB_CENTER,_("Levene's Test for Equality of Variances"));
 973   tab_joint_text(self->t, 4,0,hsize-1,0,
 974                  TAB_CENTER,_("t-test for Equality of Means"));
 975
 976   tab_text(self->t,2,2, TAB_CENTER | TAT_TITLE,_("F"));
 977   tab_text(self->t,3,2, TAB_CENTER | TAT_TITLE,_("Sig."));
 978   tab_text(self->t,4,2, TAB_CENTER | TAT_TITLE,_("t"));
 979   tab_text(self->t,5,2, TAB_CENTER | TAT_TITLE,_("df"));
 980   tab_text(self->t,6,2, TAB_CENTER | TAT_TITLE,_("Sig. (2-tailed)"));
 981   tab_text(self->t,7,2, TAB_CENTER | TAT_TITLE,_("Mean Difference"));
 982   tab_text(self->t,8,2, TAB_CENTER | TAT_TITLE,_("Std. Error Difference"));
 983   tab_text(self->t,9,2, TAB_CENTER | TAT_TITLE,_("Lower"));
 984   tab_text(self->t,10,2, TAB_CENTER | TAT_TITLE,_("Upper"));
 985
 986   tab_joint_text(self->t, 9, 1, 10, 1, TAB_CENTER | TAT_PRINTF,
 987                  _("%g%% Confidence Interval of the Difference"),
 988                  cmd->criteria*100.0);
 989
 990 }
 991
 992 /* Populate the independent samples trbox */
 993 void
 994 trbox_independent_samples_populate(struct trbox *self,
 995                                    struct cmd_t_test *cmd )
 996 {
 997   int i;
 998
 999   assert(self);
1000   for (i=0; i < cmd->n_variables; ++i)
1001     {
1002       double p,q;
1003
1004       double t;
1005       double df;
1006
1007       double df1, df2;
1008
1009       double pooled_variance;
1010       double std_err_diff;
1011       double mean_diff;
1012
1013       struct variable *var = cmd->v_variables[i];
1014       struct group_proc *grp_data = group_proc_get (var);
1015
1016       struct hsh_table *grp_hash = grp_data->group_hash;
1017
1018       struct group_statistics *gs0 ;
1019       struct group_statistics *gs1 ;
1020
1021       union value search_val;
1022
1023       if ( gp.criterion == CMP_LE )
1024         search_val.f = gp.v.critical_value - 1.0;
1025       else
1026         search_val = gp.v.g_value[0];
1027
1028       gs0 = hsh_find(grp_hash, (void *) &search_val);
1029       assert(gs0);
1030
1031       if ( gp.criterion == CMP_LE )
1032         search_val.f = gp.v.critical_value + 1.0;
1033       else
1034         search_val = gp.v.g_value[1];
1035
1036       gs1 = hsh_find(grp_hash, (void *) &search_val);
1037       assert(gs1);
1038
1039
1040       tab_text (self->t, 0, i*2+3, TAB_LEFT, cmd->v_variables[i]->name);
1041
1042       tab_text (self->t, 1, i*2+3, TAB_LEFT, _("Equal variances assumed"));
1043
1044
1045       tab_float(self->t, 2, i*2+3, TAB_CENTER, grp_data->levene, 8,3);
1046
1047       /* Now work out the significance of the Levene test */
1048       df1 = 1; df2 = grp_data->ugs.n - 2;
1049       q = gsl_cdf_fdist_Q(grp_data->levene, df1, df2);
1050
1051       tab_float(self->t, 3, i*2+3, TAB_CENTER, q, 8,3 );
1052
1053       df = gs0->n + gs1->n - 2.0 ;
1054       tab_float (self->t, 5, i*2+3, TAB_RIGHT, df, 2, 0);
1055
1056       pooled_variance = ( (gs0->n )*pow2(gs0->s_std_dev)
1057                           +
1058                           (gs1->n )*pow2(gs1->s_std_dev)
1059                         ) / df  ;
1060
1061       t = (gs0->mean - gs1->mean) / sqrt(pooled_variance) ;
1062       t /= sqrt((gs0->n + gs1->n)/(gs0->n*gs1->n));
1063
1064       tab_float (self->t, 4, i*2+3, TAB_RIGHT, t, 8, 3);
1065
1066       p = gsl_cdf_tdist_P(t, df);
1067       q = gsl_cdf_tdist_Q(t, df);
1068
1069       tab_float(self->t, 6, i*2+3, TAB_RIGHT, 2.0*(t>0?q:p) , 8, 3);
1070
1071       mean_diff = gs0->mean - gs1->mean;
1072       tab_float(self->t, 7, i*2+3, TAB_RIGHT, mean_diff, 8, 3);
1073
1074
1075       std_err_diff = sqrt( pow2(gs0->se_mean) + pow2(gs1->se_mean));
1076       tab_float(self->t, 8, i*2+3, TAB_RIGHT, std_err_diff, 8, 3);
1077
1078
1079       /* Now work out the confidence interval */
1080       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1081
1082       t = gsl_cdf_tdist_Qinv(q,df);
1083       tab_float(self->t, 9, i*2+3, TAB_RIGHT,
1084                 mean_diff - t * std_err_diff, 8, 3);
1085
1086       tab_float(self->t, 10, i*2+3, TAB_RIGHT,
1087                 mean_diff + t * std_err_diff, 8, 3);
1088
1089
1090       {
1091         double se2;
1092       /* Now for the \sigma_1 != \sigma_2 case */
1093       tab_text (self->t, 1, i*2+3+1,
1094                 TAB_LEFT, _("Equal variances not assumed"));
1095
1096
1097       se2 = (pow2(gs0->s_std_dev)/(gs0->n -1) ) +
1098         (pow2(gs1->s_std_dev)/(gs1->n -1) );
1099
1100       t = mean_diff / sqrt(se2) ;
1101       tab_float (self->t, 4, i*2+3+1, TAB_RIGHT, t, 8, 3);
1102
1103       df = pow2(se2) / (
1104                        (pow2(pow2(gs0->s_std_dev)/(gs0->n - 1 ))
1105                         /(gs0->n -1 )
1106                         )
1107                        +
1108                        (pow2(pow2(gs1->s_std_dev)/(gs1->n - 1 ))
1109                         /(gs1->n -1 )
1110                         )
1111                        ) ;
1112       tab_float (self->t, 5, i*2+3+1, TAB_RIGHT, df, 8, 3);
1113
1114       p = gsl_cdf_tdist_P(t, df);
1115       q = gsl_cdf_tdist_Q(t, df);
1116
1117       tab_float(self->t, 6, i*2+3+1, TAB_RIGHT, 2.0*(t>0?q:p) , 8, 3);
1118
1119       /* Now work out the confidence interval */
1120       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1121
1122       t = gsl_cdf_tdist_Qinv(q, df);
1123
1124       tab_float(self->t, 7, i*2+3+1, TAB_RIGHT, mean_diff, 8, 3);
1125
1126
1127       tab_float(self->t, 8, i*2+3+1, TAB_RIGHT, std_err_diff, 8, 3);
1128
1129
1130       tab_float(self->t, 9, i*2+3+1, TAB_RIGHT,
1131                 mean_diff - t * std_err_diff, 8, 3);
1132
1133       tab_float(self->t, 10, i*2+3+1, TAB_RIGHT,
1134                 mean_diff + t * std_err_diff, 8, 3);
1135
1136       }
1137     }
1138 }
1139
1140 /* Initialize the paired samples trbox */
1141 void
1142 trbox_paired_init(struct trbox *self,
1143                            struct cmd_t_test *cmd UNUSED)
1144 {
1145
1146   const int hsize=10;
1147   const int vsize=n_pairs+3;
1148
1149   self->populate = trbox_paired_populate;
1150
1151   trbox_base_init(self,n_pairs,hsize);
1152   tab_title (self->t, 0, _("Paired Samples Test"));
1153   tab_hline(self->t,TAL_1,2,6,1);
1154   tab_vline(self->t,TAL_2,2,0,vsize - 1);
1155   tab_joint_text(self->t,2,0,6,0,TAB_CENTER,_("Paired Differences"));
1156   tab_box(self->t,-1,-1,-1,TAL_1, 2,1,6,vsize-1);
1157   tab_box(self->t,-1,-1,-1,TAL_1, 6,0,hsize-1,vsize-1);
1158   tab_hline(self->t,TAL_1,5,6, 2);
1159   tab_vline(self->t,TAL_0,6,0,1);
1160
1161   tab_joint_text(self->t, 5, 1, 6, 1, TAB_CENTER | TAT_PRINTF,
1162                  _("%g%% Confidence Interval of the Difference"),
1163                  cmd->criteria*100.0);
1164
1165   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("Mean"));
1166   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
1167   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("Std. Error Mean"));
1168   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
1169   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
1170   tab_text (self->t, 7, 2, TAB_CENTER | TAT_TITLE, _("t"));
1171   tab_text (self->t, 8, 2, TAB_CENTER | TAT_TITLE, _("df"));
1172   tab_text (self->t, 9, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
1173 }
1174
1175 /* Populate the paired samples trbox */
1176 void
1177 trbox_paired_populate(struct trbox *trb,
1178                               struct cmd_t_test *cmd UNUSED)
1179 {
1180   int i;
1181
1182   for (i=0; i < n_pairs; ++i)
1183     {
1184       double p,q;
1185       double se_mean;
1186
1187       double n = pairs[i].n;
1188       double t;
1189       double df = n - 1;
1190
1191       tab_text (trb->t, 0, i+3, TAB_LEFT | TAT_PRINTF, _("Pair %d"),i);
1192
1193       tab_text (trb->t, 1, i+3, TAB_LEFT | TAT_PRINTF, "%s - %s",
1194                 pairs[i].v[0]->name, pairs[i].v[1]->name);
1195
1196       tab_float(trb->t, 2, i+3, TAB_RIGHT, pairs[i].mean_diff, 8, 4);
1197
1198       tab_float(trb->t, 3, i+3, TAB_RIGHT, pairs[i].std_dev_diff, 8, 5);
1199
1200       /* SE Mean */
1201       se_mean = pairs[i].std_dev_diff / sqrt(n) ;
1202       tab_float(trb->t, 4, i+3, TAB_RIGHT, se_mean, 8,5 );
1203
1204       /* Now work out the confidence interval */
1205       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1206
1207       t = gsl_cdf_tdist_Qinv(q, df);
1208
1209       tab_float(trb->t, 5, i+3, TAB_RIGHT,
1210                 pairs[i].mean_diff - t * se_mean , 8, 4);
1211
1212       tab_float(trb->t, 6, i+3, TAB_RIGHT,
1213                 pairs[i].mean_diff + t * se_mean , 8, 4);
1214
1215       t = (pairs[i].mean[0] - pairs[i].mean[1])
1216         / sqrt (
1217                 ( pow2 (pairs[i].s_std_dev[0]) + pow2 (pairs[i].s_std_dev[1]) -
1218                   2 * pairs[i].correlation *
1219                   pairs[i].s_std_dev[0] * pairs[i].s_std_dev[1] )
1220                 / (n - 1)
1221                 );
1222
1223       tab_float(trb->t, 7, i+3, TAB_RIGHT, t , 8,3 );
1224
1225       /* Degrees of freedom */
1226       tab_float(trb->t, 8, i+3, TAB_RIGHT, df , 2, 0 );
1227
1228       p = gsl_cdf_tdist_P(t,df);
1229       q = gsl_cdf_tdist_P(t,df);
1230
1231       tab_float(trb->t, 9, i+3, TAB_RIGHT, 2.0*(t>0?q:p) , 8, 3);
1232
1233     }
1234 }
1235
1236 /* Initialize the one sample trbox */
1237 void
1238 trbox_one_sample_init(struct trbox *self, struct cmd_t_test *cmd )
1239 {
1240   const int hsize=7;
1241   const int vsize=cmd->n_variables+3;
1242
1243   self->populate = trbox_one_sample_populate;
1244
1245   trbox_base_init(self, cmd->n_variables,hsize);
1246   tab_title (self->t, 0, _("One-Sample Test"));
1247   tab_hline(self->t, TAL_1, 1, hsize - 1, 1);
1248   tab_vline(self->t, TAL_2, 1, 0, vsize - 1);
1249
1250   tab_joint_text(self->t, 1, 0, hsize-1,0, TAB_CENTER | TAT_PRINTF,
1251                  _("Test Value = %f"), cmd->n_testval[0]);
1252
1253   tab_box(self->t, -1, -1, -1, TAL_1, 1,1,hsize-1,vsize-1);
1254
1255
1256   tab_joint_text(self->t,5,1,6,1,TAB_CENTER  | TAT_PRINTF,
1257                  _("%g%% Confidence Interval of the Difference"),
1258                  cmd->criteria*100.0);
1259
1260   tab_vline(self->t,TAL_0,6,1,1);
1261   tab_hline(self->t,TAL_1,5,6,2);
1262   tab_text (self->t, 1, 2, TAB_CENTER | TAT_TITLE, _("t"));
1263   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("df"));
1264   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
1265   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
1266   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
1267   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
1268
1269 }
1270
1271
1272 /* Populate the one sample trbox */
1273 void
1274 trbox_one_sample_populate(struct trbox *trb, struct cmd_t_test *cmd)
1275 {
1276   int i;
1277
1278   assert(trb->t);
1279
1280   for (i=0; i < cmd->n_variables; ++i)
1281     {
1282       double t;
1283       double p,q;
1284       double df;
1285       struct group_statistics *gs = &group_proc_get (cmd->v_variables[i])->ugs;
1286
1287
1288       tab_text (trb->t, 0, i+3, TAB_LEFT, cmd->v_variables[i]->name);
1289
1290       t = (gs->mean - cmd->n_testval[0] ) * sqrt(gs->n) / gs->std_dev ;
1291
1292       tab_float (trb->t, 1, i+3, TAB_RIGHT, t, 8,3);
1293
1294       /* degrees of freedom */
1295       df = gs->n - 1;
1296
1297       tab_float (trb->t, 2, i+3, TAB_RIGHT, df, 8,0);
1298
1299       p = gsl_cdf_tdist_P(t, df);
1300       q = gsl_cdf_tdist_Q(t, df);
1301
1302       /* Multiply by 2 to get 2-tailed significance, makeing sure we've got
1303          the correct tail*/
1304       tab_float (trb->t, 3, i+3, TAB_RIGHT, 2.0*(t>0?q:p), 8,3);
1305
1306       tab_float (trb->t, 4, i+3, TAB_RIGHT, gs->mean_diff, 8,3);
1307
1308
1309       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1310       t = gsl_cdf_tdist_Qinv(q, df);
1311
1312       tab_float (trb->t, 5, i+3, TAB_RIGHT,
1313                  gs->mean_diff - t * gs->se_mean, 8,4);
1314
1315       tab_float (trb->t, 6, i+3, TAB_RIGHT,
1316                  gs->mean_diff + t * gs->se_mean, 8,4);
1317     }
1318 }
1319
1320 /* Base initializer for the generalized trbox */
1321 void
1322 trbox_base_init(struct trbox *self, size_t data_rows, int cols)
1323 {
1324   const size_t rows = 3 + data_rows;
1325
1326   self->finalize = trbox_base_finalize;
1327   self->t = tab_create (cols, rows, 0);
1328   tab_headers (self->t,0,0,3,0);
1329   tab_box (self->t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols -1, rows -1);
1330   tab_hline(self->t, TAL_2,0,cols-1,3);
1331   tab_dim (self->t, tab_natural_dimensions);
1332 }
1333
1334
1335 /* Base finalizer for the trbox */
1336 void
1337 trbox_base_finalize(struct trbox *trb)
1338 {
1339   tab_submit(trb->t);
1340 }
1341
1342
1343 /* Create , populate and submit the Paired Samples Correlation box */
1344 void
1345 pscbox(void)
1346 {
1347   const int rows=1+n_pairs;
1348   const int cols=5;
1349   int i;
1350
1351   struct tab_table *table;
1352
1353   table = tab_create (cols,rows,0);
1354
1355   tab_columns (table, SOM_COL_DOWN, 1);
1356   tab_headers (table,0,0,1,0);
1357   tab_box (table, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols -1, rows -1 );
1358   tab_hline(table, TAL_2, 0, cols - 1, 1);
1359   tab_vline(table, TAL_2, 2, 0, rows - 1);
1360   tab_dim(table, tab_natural_dimensions);
1361   tab_title(table, 0, _("Paired Samples Correlations"));
1362
1363   /* column headings */
1364   tab_text(table, 2,0, TAB_CENTER | TAT_TITLE, _("N"));
1365   tab_text(table, 3,0, TAB_CENTER | TAT_TITLE, _("Correlation"));
1366   tab_text(table, 4,0, TAB_CENTER | TAT_TITLE, _("Sig."));
1367
1368   for (i=0; i < n_pairs; ++i)
1369     {
1370       double p,q;
1371
1372       double df = pairs[i].n -2;
1373
1374       double correlation_t =
1375         pairs[i].correlation * sqrt(df) /
1376         sqrt(1 - pow2(pairs[i].correlation));
1377
1378
1379       /* row headings */
1380       tab_text(table, 0,i+1, TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1381                _("Pair %d"), i);
1382
1383       tab_text(table, 1,i+1, TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1384                _("%s & %s"), pairs[i].v[0]->name, pairs[i].v[1]->name);
1385
1386
1387       /* row data */
1388       tab_float(table, 2, i+1, TAB_RIGHT, pairs[i].n, 4, 0);
1389       tab_float(table, 3, i+1, TAB_RIGHT, pairs[i].correlation, 8, 3);
1390
1391       p = gsl_cdf_tdist_P(correlation_t, df);
1392       q = gsl_cdf_tdist_Q(correlation_t, df);
1393
1394       tab_float(table, 4, i+1, TAB_RIGHT, 2.0*(correlation_t>0?q:p), 8, 3);
1395     }
1396
1397   tab_submit(table);
1398 }
1399
1400
1401
1402
1403 /* Calculation Implementation */
1404
1405 /* Per case calculations common to all variants of the T test */
1406 static int
1407 common_calc (const struct ccase *c, void *_cmd)
1408 {
1409   int i;
1410   struct cmd_t_test *cmd = (struct cmd_t_test *)_cmd;
1411
1412   double weight = dict_get_case_weight(default_dict,c,&bad_weight_warn);
1413
1414
1415   /* Skip the entire case if /MISSING=LISTWISE is set */
1416   if ( cmd->miss == TTS_LISTWISE )
1417     {
1418       for(i=0; i< cmd->n_variables ; ++i)
1419         {
1420           struct variable *v = cmd->v_variables[i];
1421           const union value *val = case_data (c, v->fv);
1422
1423           if (value_is_missing(&v->miss, val) )
1424             {
1425               return 0;
1426             }
1427         }
1428     }
1429
1430   /* Listwise has to be implicit if the independent variable is missing ?? */
1431   if ( cmd->sbc_groups )
1432     {
1433       const union value *gv = case_data (c, indep_var->fv);
1434       if ( value_is_missing(&indep_var->miss, gv) )
1435         {
1436           return 0;
1437         }
1438     }
1439
1440
1441   for(i=0; i< cmd->n_variables ; ++i)
1442     {
1443       struct group_statistics *gs;
1444       struct variable *v = cmd->v_variables[i];
1445       const union value *val = case_data (c, v->fv);
1446
1447       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1448
1449       if (! value_is_missing(&v->miss, val) )
1450         {
1451           gs->n+=weight;
1452           gs->sum+=weight * val->f;
1453           gs->ssq+=weight * val->f * val->f;
1454         }
1455     }
1456   return 0;
1457 }
1458
1459 /* Pre calculations common to all variants of the T test */
1460 static void
1461 common_precalc ( struct cmd_t_test *cmd )
1462 {
1463   int i=0;
1464
1465   for(i=0; i< cmd->n_variables ; ++i)
1466     {
1467       struct group_statistics *gs;
1468       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1469
1470       gs->sum=0;
1471       gs->n=0;
1472       gs->ssq=0;
1473       gs->sum_diff=0;
1474     }
1475 }
1476
1477 /* Post calculations common to all variants of the T test */
1478 void
1479 common_postcalc (  struct cmd_t_test *cmd )
1480 {
1481   int i=0;
1482
1483
1484   for(i=0; i< cmd->n_variables ; ++i)
1485     {
1486       struct group_statistics *gs;
1487       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1488
1489       gs->mean=gs->sum / gs->n;
1490       gs->s_std_dev= sqrt(
1491                          ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1492                          ) ;
1493
1494       gs->std_dev= sqrt(
1495                          gs->n/(gs->n-1) *
1496                          ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1497                          ) ;
1498
1499       gs->se_mean = gs->std_dev / sqrt(gs->n);
1500       gs->mean_diff= gs->sum_diff / gs->n;
1501     }
1502 }
1503
1504 /* Per case calculations for one sample t test  */
1505 static int
1506 one_sample_calc (const struct ccase *c, void *cmd_)
1507 {
1508   int i;
1509   struct cmd_t_test *cmd = (struct cmd_t_test *)cmd_;
1510
1511
1512   double weight = dict_get_case_weight(default_dict,c,&bad_weight_warn);
1513
1514   /* Skip the entire case if /MISSING=LISTWISE is set */
1515   if ( cmd->miss == TTS_LISTWISE )
1516     {
1517       for(i=0; i< cmd->n_variables ; ++i)
1518         {
1519           struct variable *v = cmd->v_variables[i];
1520           const union value *val = case_data (c, v->fv);
1521
1522           if (value_is_missing(&v->miss, val) )
1523             {
1524               return 0;
1525             }
1526         }
1527     }
1528
1529   for(i=0; i< cmd->n_variables ; ++i)
1530     {
1531       struct group_statistics *gs;
1532       struct variable *v = cmd->v_variables[i];
1533       const union value *val = case_data (c, v->fv);
1534
1535       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1536
1537       if ( ! value_is_missing(&v->miss, val))
1538         gs->sum_diff += weight * (val->f - cmd->n_testval[0]);
1539     }
1540
1541   return 0;
1542 }
1543
1544 /* Pre calculations for one sample t test */
1545 static void
1546 one_sample_precalc ( struct cmd_t_test *cmd )
1547 {
1548   int i=0;
1549
1550   for(i=0; i< cmd->n_variables ; ++i)
1551     {
1552       struct group_statistics *gs;
1553       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1554
1555       gs->sum_diff=0;
1556     }
1557 }
1558
1559 /* Post calculations for one sample t test */
1560 static void
1561 one_sample_postcalc (struct cmd_t_test *cmd)
1562 {
1563   int i=0;
1564
1565   for(i=0; i< cmd->n_variables ; ++i)
1566     {
1567       struct group_statistics *gs;
1568       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1569
1570       gs->mean_diff = gs->sum_diff / gs->n ;
1571     }
1572 }
1573
1574
1575
1576 static void
1577 paired_precalc (struct cmd_t_test *cmd UNUSED)
1578 {
1579   int i;
1580
1581   for(i=0; i < n_pairs ; ++i )
1582     {
1583       pairs[i].n = 0;
1584       pairs[i].sum[0] = 0;      pairs[i].sum[1] = 0;
1585       pairs[i].ssq[0] = 0;      pairs[i].ssq[1] = 0;
1586       pairs[i].sum_of_prod = 0;
1587       pairs[i].correlation = 0;
1588       pairs[i].sum_of_diffs = 0;
1589       pairs[i].ssq_diffs = 0;
1590     }
1591
1592 }
1593
1594
1595 static int
1596 paired_calc (const struct ccase *c, void *cmd_)
1597 {
1598   int i;
1599
1600   struct cmd_t_test *cmd  = (struct cmd_t_test *) cmd_;
1601
1602   double weight = dict_get_case_weight(default_dict,c,&bad_weight_warn);
1603
1604   /* Skip the entire case if /MISSING=LISTWISE is set ,
1605    AND one member of a pair is missing */
1606   if ( cmd->miss == TTS_LISTWISE )
1607     {
1608       for(i=0; i < n_pairs ; ++i )
1609         {
1610           struct variable *v0 = pairs[i].v[0];
1611           struct variable *v1 = pairs[i].v[1];
1612
1613           const union value *val0 = case_data (c, v0->fv);
1614           const union value *val1 = case_data (c, v1->fv);
1615
1616           if ( value_is_missing(&v0->miss, val0) ||
1617                value_is_missing(&v1->miss, val1) )
1618             {
1619               return 0;
1620             }
1621         }
1622     }
1623
1624   for(i=0; i < n_pairs ; ++i )
1625     {
1626       struct variable *v0 = pairs[i].v[0];
1627       struct variable *v1 = pairs[i].v[1];
1628
1629       const union value *val0 = case_data (c, v0->fv);
1630       const union value *val1 = case_data (c, v1->fv);
1631
1632       if ( ( !value_is_missing(&v0->miss, val0)
1633              && !value_is_missing(&v1->miss, val1) ) )
1634       {
1635         pairs[i].n += weight;
1636         pairs[i].sum[0] += weight * val0->f;
1637         pairs[i].sum[1] += weight * val1->f;
1638
1639         pairs[i].ssq[0] += weight * pow2(val0->f);
1640         pairs[i].ssq[1] += weight * pow2(val1->f);
1641
1642         pairs[i].sum_of_prod += weight * val0->f * val1->f ;
1643
1644         pairs[i].sum_of_diffs += weight * ( val0->f - val1->f ) ;
1645         pairs[i].ssq_diffs += weight * pow2(val0->f - val1->f);
1646       }
1647     }
1648
1649   return 0;
1650 }
1651
1652 static void
1653 paired_postcalc (struct cmd_t_test *cmd UNUSED)
1654 {
1655   int i;
1656
1657   for(i=0; i < n_pairs ; ++i )
1658     {
1659       int j;
1660       const double n = pairs[i].n;
1661
1662       for (j=0; j < 2 ; ++j)
1663         {
1664           pairs[i].mean[j] = pairs[i].sum[j] / n ;
1665           pairs[i].s_std_dev[j] = sqrt((pairs[i].ssq[j] / n -
1666                                               pow2(pairs[i].mean[j]))
1667                                      );
1668
1669           pairs[i].std_dev[j] = sqrt(n/(n-1)*(pairs[i].ssq[j] / n -
1670                                               pow2(pairs[i].mean[j]))
1671                                      );
1672         }
1673
1674       pairs[i].correlation = pairs[i].sum_of_prod / pairs[i].n -
1675         pairs[i].mean[0] * pairs[i].mean[1] ;
1676       /* correlation now actually contains the covariance */
1677
1678       pairs[i].correlation /= pairs[i].std_dev[0] * pairs[i].std_dev[1];
1679       pairs[i].correlation *= pairs[i].n / ( pairs[i].n - 1 );
1680
1681       pairs[i].mean_diff = pairs[i].sum_of_diffs / n ;
1682
1683       pairs[i].std_dev_diff = sqrt (  n / (n - 1) * (
1684                                     ( pairs[i].ssq_diffs / n )
1685                                     -
1686                                     pow2(pairs[i].mean_diff )
1687                                     ) );
1688     }
1689 }
1690
1691 static void
1692 group_precalc (struct cmd_t_test *cmd )
1693 {
1694   int i;
1695   int j;
1696
1697   for(i=0; i< cmd->n_variables ; ++i)
1698     {
1699       struct group_proc *ttpr = group_proc_get (cmd->v_variables[i]);
1700
1701       /* There's always 2 groups for a T - TEST */
1702       ttpr->n_groups = 2;
1703
1704       gp.indep_width = indep_var->width;
1705
1706       ttpr->group_hash = hsh_create(2,
1707                                     (hsh_compare_func *) compare_group_binary,
1708                                     (hsh_hash_func *) hash_group_binary,
1709                                     (hsh_free_func *) free_group,
1710                                     (void *) &gp );
1711
1712       for (j=0 ; j < 2 ; ++j)
1713         {
1714
1715           struct group_statistics *gs = xmalloc (sizeof *gs);
1716
1717           gs->sum = 0;
1718           gs->n = 0;
1719           gs->ssq = 0;
1720
1721           if ( gp.criterion == CMP_EQ )
1722             {
1723               gs->id = gp.v.g_value[j];
1724             }
1725           else
1726             {
1727               if ( j == 0 )
1728                 gs->id.f = gp.v.critical_value - 1.0 ;
1729               else
1730                 gs->id.f = gp.v.critical_value + 1.0 ;
1731             }
1732
1733           hsh_insert ( ttpr->group_hash, (void *) gs );
1734
1735         }
1736     }
1737
1738 }
1739
1740 static int
1741 group_calc (const struct ccase *c, struct cmd_t_test *cmd)
1742 {
1743   int i;
1744
1745   const union value *gv = case_data (c, indep_var->fv);
1746
1747   const double weight = dict_get_case_weight(default_dict,c,&bad_weight_warn);
1748
1749   if ( value_is_missing(&indep_var->miss, gv) )
1750     {
1751       return 0;
1752     }
1753
1754   if ( cmd->miss == TTS_LISTWISE )
1755     {
1756       for(i=0; i< cmd->n_variables ; ++i)
1757         {
1758           struct variable *v = cmd->v_variables[i];
1759           const union value *val = case_data (c, v->fv);
1760
1761           if (value_is_missing(&v->miss, val) )
1762             {
1763               return 0;
1764             }
1765         }
1766     }
1767
1768   gv = case_data (c, indep_var->fv);
1769
1770   for(i=0; i< cmd->n_variables ; ++i)
1771     {
1772       struct variable *var = cmd->v_variables[i];
1773       const union value *val = case_data (c, var->fv);
1774       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1775       struct group_statistics *gs;
1776
1777       gs = hsh_find(grp_hash, (void *) gv);
1778
1779       /* If the independent variable doesn't match either of the values
1780          for this case then move on to the next case */
1781       if ( ! gs )
1782         return 0;
1783
1784       if ( !value_is_missing(&var->miss, val) )
1785         {
1786           gs->n+=weight;
1787           gs->sum+=weight * val->f;
1788           gs->ssq+=weight * pow2(val->f);
1789         }
1790     }
1791
1792   return 0;
1793 }
1794
1795
1796 static void
1797 group_postcalc ( struct cmd_t_test *cmd )
1798 {
1799   int i;
1800
1801   for(i=0; i< cmd->n_variables ; ++i)
1802     {
1803       struct variable *var = cmd->v_variables[i];
1804       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1805       struct hsh_iterator g;
1806       struct group_statistics *gs;
1807       int count=0;
1808
1809       for (gs =  hsh_first (grp_hash,&g);
1810            gs != 0;
1811            gs = hsh_next(grp_hash,&g))
1812         {
1813           gs->mean = gs->sum / gs->n;
1814
1815           gs->s_std_dev= sqrt(
1816                               ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1817                               ) ;
1818
1819           gs->std_dev= sqrt(
1820                             gs->n/(gs->n-1) *
1821                             ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1822                             ) ;
1823
1824           gs->se_mean = gs->std_dev / sqrt(gs->n);
1825           count ++;
1826         }
1827       assert(count == 2);
1828     }
1829 }
1830
1831
1832
1833 static bool
1834 calculate(const struct casefile *cf, void *cmd_)
1835 {
1836   struct ssbox stat_summary_box;
1837   struct trbox test_results_box;
1838
1839   struct casereader *r;
1840   struct ccase c;
1841
1842   struct cmd_t_test *cmd = (struct cmd_t_test *) cmd_;
1843
1844   common_precalc(cmd);
1845   for(r = casefile_get_reader (cf);
1846       casereader_read (r, &c) ;
1847       case_destroy (&c))
1848     {
1849       common_calc(&c,cmd);
1850     }
1851   casereader_destroy (r);
1852   common_postcalc(cmd);
1853
1854   switch(mode)
1855     {
1856     case T_1_SAMPLE:
1857       one_sample_precalc(cmd);
1858       for(r = casefile_get_reader (cf);
1859           casereader_read (r, &c) ;
1860           case_destroy (&c))
1861         {
1862           one_sample_calc(&c,cmd);
1863         }
1864       casereader_destroy (r);
1865       one_sample_postcalc(cmd);
1866
1867       break;
1868     case T_PAIRED:
1869       paired_precalc(cmd);
1870       for(r = casefile_get_reader (cf);
1871           casereader_read (r, &c) ;
1872           case_destroy (&c))
1873         {
1874           paired_calc(&c,cmd);
1875         }
1876       casereader_destroy (r);
1877       paired_postcalc(cmd);
1878
1879       break;
1880     case T_IND_SAMPLES:
1881
1882       group_precalc(cmd);
1883       for(r = casefile_get_reader (cf);
1884           casereader_read (r, &c) ;
1885           case_destroy (&c))
1886         {
1887           group_calc(&c,cmd);
1888         }
1889       casereader_destroy (r);
1890       group_postcalc(cmd);
1891
1892       levene(cf, indep_var, cmd->n_variables, cmd->v_variables,
1893              (cmd->miss == TTS_LISTWISE)?LEV_LISTWISE:LEV_ANALYSIS ,
1894              value_is_missing);
1895       break;
1896     }
1897
1898   ssbox_create(&stat_summary_box,cmd,mode);
1899   ssbox_populate(&stat_summary_box,cmd);
1900   ssbox_finalize(&stat_summary_box);
1901
1902   if ( mode == T_PAIRED)
1903       pscbox();
1904
1905   trbox_create(&test_results_box,cmd,mode);
1906   trbox_populate(&test_results_box,cmd);
1907   trbox_finalize(&test_results_box);
1908
1909   return true;
1910 }
1911
1912 short which_group(const struct group_statistics *g,
1913                   const struct group_properties *p);
1914
1915 /* Return -1 if the id of a is less than b; +1 if greater than and
1916    0 if equal */
1917 static int
1918 compare_group_binary(const struct group_statistics *a,
1919                      const struct group_statistics *b,
1920                      const struct group_properties *p)
1921 {
1922   short flag_a;
1923   short flag_b;
1924
1925   if ( p->criterion == CMP_LE )
1926     {
1927       /* less-than-or-equal comparision is not meaningfull for
1928          alpha variables, so we shouldn't ever arrive here */
1929       assert(p->indep_width == 0 ) ;
1930
1931       flag_a = ( a->id.f < p->v.critical_value ) ;
1932       flag_b = ( b->id.f < p->v.critical_value ) ;
1933     }
1934   else
1935     {
1936       flag_a = which_group(a, p);
1937       flag_b = which_group(b, p);
1938     }
1939
1940   if (flag_a < flag_b )
1941     return -1;
1942
1943   return (flag_a > flag_b);
1944 }
1945
1946 /* This is a degenerate case of a hash, since it can only return three possible
1947    values.  It's really a comparison, being used as a hash function */
1948
1949 static unsigned
1950 hash_group_binary(const struct group_statistics *g,
1951                   const struct group_properties *p)
1952 {
1953   short flag = -1;
1954
1955   if ( p->criterion == CMP_LE )
1956     {
1957       /* Not meaningfull to do a less than compare for alpha values ? */
1958       assert(p->indep_width == 0 ) ;
1959       flag = ( g->id.f < p->v.critical_value ) ;
1960     }
1961   else if ( p->criterion == CMP_EQ)
1962     {
1963       flag = which_group(g,p);
1964     }
1965   else
1966     assert(0);
1967
1968   return flag;
1969 }
1970
1971 /* return 0 if G belongs to group 0,
1972           1 if it belongs to group 1,
1973           2 if it belongs to neither group */
1974 short
1975 which_group(const struct group_statistics *g,
1976             const struct group_properties *p)
1977 {
1978
1979   if ( 0 == compare_values (&g->id, &p->v.g_value[0], p->indep_width))
1980     return 0;
1981
1982   if ( 0 == compare_values (&g->id, &p->v.g_value[1], p->indep_width))
1983     return 1;
1984
1985   return 2;
1986 }
1987