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