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