pintos-os.org Git - pspp/blob - src/language/stats/oneway.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-9, 2000, 2007, 2009, 2010 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <data/case.h>
  20 #include <data/casegrouper.h>
  21 #include <data/casereader.h>
  22
  23 #include <libpspp/ll.h>
  24
  25 #include <language/lexer/lexer.h>
  26 #include <language/lexer/variable-parser.h>
  27 #include <language/lexer/value-parser.h>
  28 #include <language/command.h>
  29
  30 #include <data/procedure.h>
  31 #include <data/dictionary.h>
  32
  33
  34 #include <language/dictionary/split-file.h>
  35 #include <libpspp/hash.h>
  36 #include <libpspp/taint.h>
  37 #include <math/group-proc.h>
  38 #include <math/levene.h>
  39 #include <libpspp/misc.h>
  40
  41 #include <output/tab.h>
  42
  43 #include <gsl/gsl_cdf.h>
  44 #include <math.h>
  45 #include <data/format.h>
  46
  47 #include <libpspp/message.h>
  48
  49 #include "gettext.h"
  50 #define _(msgid) gettext (msgid)
  51
  52 enum missing_type
  53   {
  54     MISS_LISTWISE,
  55     MISS_ANALYSIS,
  56   };
  57
  58 enum statistics
  59   {
  60     STATS_DESCRIPTIVES = 0x0001,
  61     STATS_HOMOGENEITY = 0x0002
  62   };
  63
  64 struct coeff_node
  65 {
  66   struct ll ll;
  67   double coeff;
  68 };
  69
  70
  71 struct contrasts_node
  72 {
  73   struct ll ll;
  74   struct ll_list coefficient_list;
  75
  76   bool bad_count; /* True if the number of coefficients does not equal the number of groups */
  77 };
  78
  79 struct oneway
  80 {
  81   size_t n_vars;
  82   const struct variable **vars;
  83
  84   const struct variable *indep_var;
  85
  86   enum statistics stats;
  87
  88   enum missing_type missing_type;
  89   enum mv_class exclude;
  90
  91   /* The number of distinct values of the independent variable, when all
  92      missing values are disregarded */
  93   int actual_number_of_groups;
  94
  95   /* A  hash table containing all the distinct values of the independent
  96      variable */
  97   struct hsh_table *group_hash;
  98
  99   /* List of contrasts */
 100   struct ll_list contrast_list;
 101 };
 102
 103 /* Routines to show the output tables */
 104 static void show_anova_table (const struct oneway *);
 105 static void show_descriptives (const struct oneway *, const struct dictionary *dict);
 106 static void show_homogeneity (const struct oneway *);
 107
 108 static void output_oneway (const struct oneway *, const struct dictionary *dict);
 109 static void run_oneway (struct oneway *cmd, struct casereader *input, const struct dataset *ds);
 110
 111 int
 112 cmd_oneway (struct lexer *lexer, struct dataset *ds)
 113 {
 114   const struct dictionary *dict = dataset_dict (ds);
 115
 116   struct oneway oneway ;
 117   oneway.n_vars = 0;
 118   oneway.vars = NULL;
 119   oneway.indep_var = NULL;
 120   oneway.stats = 0;
 121   oneway.missing_type = MISS_ANALYSIS;
 122   oneway.exclude = MV_ANY;
 123   oneway.actual_number_of_groups = 0;
 124   oneway.group_hash = NULL;
 125
 126   ll_init (&oneway.contrast_list);
 127
 128
 129   if ( lex_match (lexer, '/'))
 130     {
 131       if (!lex_force_match_id (lexer, "VARIABLES"))
 132         {
 133           goto error;
 134         }
 135       lex_match (lexer, '=');
 136     }
 137
 138   if (!parse_variables_const (lexer, dict,
 139                               &oneway.vars, &oneway.n_vars,
 140                               PV_NO_DUPLICATE | PV_NUMERIC))
 141     goto error;
 142
 143   lex_force_match (lexer, T_BY);
 144
 145   oneway.indep_var = parse_variable_const (lexer, dict);
 146
 147   while (lex_token (lexer) != '.')
 148     {
 149       lex_match (lexer, '/');
 150
 151       if (lex_match_id (lexer, "STATISTICS"))
 152         {
 153           lex_match (lexer, '=');
 154           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 155             {
 156               if (lex_match_id (lexer, "DESCRIPTIVES"))
 157                 {
 158                   oneway.stats |= STATS_DESCRIPTIVES;
 159                 }
 160               else if (lex_match_id (lexer, "HOMOGENEITY"))
 161                 {
 162                   oneway.stats |= STATS_HOMOGENEITY;
 163                 }
 164               else
 165                 {
 166                   lex_error (lexer, NULL);
 167                   goto error;
 168                 }
 169             }
 170         }
 171       else if (lex_match_id (lexer, "CONTRAST"))
 172         {
 173           struct contrasts_node *cl = xzalloc (sizeof *cl);
 174
 175           struct ll_list *coefficient_list = &cl->coefficient_list;
 176           lex_match (lexer, '=');
 177
 178           ll_init (coefficient_list);
 179
 180           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 181             {
 182               union value val;
 183               if ( parse_value (lexer, &val, 0))
 184                 {
 185                   struct coeff_node *cc = xmalloc (sizeof *cc);
 186                   cc->coeff = val.f;
 187
 188                   ll_push_tail (coefficient_list, &cc->ll);
 189                 }
 190               else
 191                 {
 192                   lex_error (lexer, NULL);
 193                   goto error;
 194                 }
 195             }
 196
 197           ll_push_tail (&oneway.contrast_list, &cl->ll);
 198         }
 199       else if (lex_match_id (lexer, "MISSING"))
 200         {
 201           lex_match (lexer, '=');
 202           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 203             {
 204               if (lex_match_id (lexer, "INCLUDE"))
 205                 {
 206                   oneway.exclude = MV_SYSTEM;
 207                 }
 208               else if (lex_match_id (lexer, "EXCLUDE"))
 209                 {
 210                   oneway.exclude = MV_ANY;
 211                 }
 212               else if (lex_match_id (lexer, "LISTWISE"))
 213                 {
 214                   oneway.missing_type = MISS_LISTWISE;
 215                 }
 216               else if (lex_match_id (lexer, "ANALYSIS"))
 217                 {
 218                   oneway.missing_type = MISS_ANALYSIS;
 219                 }
 220               else
 221                 {
 222                   lex_error (lexer, NULL);
 223                   goto error;
 224                 }
 225             }
 226         }
 227       else
 228         {
 229           lex_error (lexer, NULL);
 230           goto error;
 231         }
 232     }
 233
 234
 235   {
 236     struct casegrouper *grouper;
 237     struct casereader *group;
 238     bool ok;
 239
 240     grouper = casegrouper_create_splits (proc_open (ds), dict);
 241     while (casegrouper_get_next_group (grouper, &group))
 242       run_oneway (&oneway, group, ds);
 243     ok = casegrouper_destroy (grouper);
 244     ok = proc_commit (ds) && ok;
 245   }
 246
 247   free (oneway.vars);
 248   return CMD_SUCCESS;
 249
 250  error:
 251   free (oneway.vars);
 252   return CMD_FAILURE;
 253 }
 254
 255
 256 \f
 257
 258 static int
 259 compare_double_3way (const void *a_, const void *b_, const void *aux UNUSED)
 260 {
 261   const double *a = a_;
 262   const double *b = b_;
 263   return *a < *b ? -1 : *a > *b;
 264 }
 265
 266 static unsigned
 267 do_hash_double (const void *value_, const void *aux UNUSED)
 268 {
 269   const double *value = value_;
 270   return hash_double (*value, 0);
 271 }
 272
 273 static void
 274 free_double (void *value_, const void *aux UNUSED)
 275 {
 276   double *value = value_;
 277   free (value);
 278 }
 279
 280 static void postcalc (const struct oneway *cmd);
 281 static void  precalc (const struct oneway *cmd);
 282
 283
 284 static void
 285 run_oneway (struct oneway *cmd,
 286             struct casereader *input,
 287             const struct dataset *ds)
 288 {
 289   struct taint *taint;
 290   struct dictionary *dict = dataset_dict (ds);
 291   struct casereader *reader;
 292   struct ccase *c;
 293
 294   c = casereader_peek (input, 0);
 295   if (c == NULL)
 296     {
 297       casereader_destroy (input);
 298       return;
 299     }
 300   output_split_file_values (ds, c);
 301   case_unref (c);
 302
 303   taint = taint_clone (casereader_get_taint (input));
 304
 305   cmd->group_hash = hsh_create (4,
 306                                   compare_double_3way,
 307                                   do_hash_double,
 308                                   free_double,
 309                                   cmd->indep_var);
 310
 311   precalc (cmd);
 312
 313   input = casereader_create_filter_missing (input, &cmd->indep_var, 1,
 314                                             cmd->exclude, NULL, NULL);
 315   if (cmd->missing_type == MISS_LISTWISE)
 316     input = casereader_create_filter_missing (input, cmd->vars, cmd->n_vars,
 317                                               cmd->exclude, NULL, NULL);
 318   input = casereader_create_filter_weight (input, dict, NULL, NULL);
 319
 320   reader = casereader_clone (input);
 321   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
 322     {
 323       size_t i;
 324
 325       const double weight = dict_get_case_weight (dict, c, NULL);
 326
 327       const union value *indep_val = case_data (c, cmd->indep_var);
 328       void **p = hsh_probe (cmd->group_hash, &indep_val->f);
 329       if (*p == NULL)
 330         {
 331           double *value = *p = xmalloc (sizeof *value);
 332           *value = indep_val->f;
 333         }
 334
 335       for (i = 0; i < cmd->n_vars; ++i)
 336         {
 337           const struct variable *v = cmd->vars[i];
 338
 339           const union value *val = case_data (c, v);
 340
 341           struct group_proc *gp = group_proc_get (cmd->vars[i]);
 342           struct hsh_table *group_hash = gp->group_hash;
 343
 344           struct group_statistics *gs;
 345
 346           gs = hsh_find (group_hash, indep_val );
 347
 348           if ( ! gs )
 349             {
 350               gs = xmalloc (sizeof *gs);
 351               gs->id = *indep_val;
 352               gs->sum = 0;
 353               gs->n = 0;
 354               gs->ssq = 0;
 355               gs->sum_diff = 0;
 356               gs->minimum = DBL_MAX;
 357               gs->maximum = -DBL_MAX;
 358
 359               hsh_insert ( group_hash, gs );
 360             }
 361
 362           if (!var_is_value_missing (v, val, cmd->exclude))
 363             {
 364               struct group_statistics *totals = &gp->ugs;
 365
 366               totals->n += weight;
 367               totals->sum += weight * val->f;
 368               totals->ssq += weight * pow2 (val->f);
 369
 370               if ( val->f * weight  < totals->minimum )
 371                 totals->minimum = val->f * weight;
 372
 373               if ( val->f * weight  > totals->maximum )
 374                 totals->maximum = val->f * weight;
 375
 376               gs->n += weight;
 377               gs->sum += weight * val->f;
 378               gs->ssq += weight * pow2 (val->f);
 379
 380               if ( val->f * weight  < gs->minimum )
 381                 gs->minimum = val->f * weight;
 382
 383               if ( val->f * weight  > gs->maximum )
 384                 gs->maximum = val->f * weight;
 385             }
 386
 387           gp->n_groups = hsh_count (group_hash );
 388         }
 389
 390     }
 391   casereader_destroy (reader);
 392
 393   postcalc (cmd);
 394
 395   if ( cmd->stats & STATS_HOMOGENEITY )
 396     levene (dict, casereader_clone (input), cmd->indep_var,
 397             cmd->n_vars, cmd->vars, cmd->exclude);
 398
 399   casereader_destroy (input);
 400
 401   cmd->actual_number_of_groups = hsh_count (cmd->group_hash);
 402
 403   if (!taint_has_tainted_successor (taint))
 404     output_oneway (cmd, dict);
 405
 406   taint_destroy (taint);
 407 }
 408
 409 /* Pre calculations */
 410 static void
 411 precalc (const struct oneway *cmd)
 412 {
 413   size_t i = 0;
 414
 415   for (i = 0; i < cmd->n_vars; ++i)
 416     {
 417       struct group_proc *gp = group_proc_get (cmd->vars[i]);
 418       struct group_statistics *totals = &gp->ugs;
 419
 420       /* Create a hash for each of the dependent variables.
 421          The hash contains a group_statistics structure,
 422          and is keyed by value of the independent variable */
 423
 424       gp->group_hash = hsh_create (4, compare_group, hash_group,
 425                                    (hsh_free_func *) free_group,
 426                                    cmd->indep_var);
 427
 428       totals->sum = 0;
 429       totals->n = 0;
 430       totals->ssq = 0;
 431       totals->sum_diff = 0;
 432       totals->maximum = -DBL_MAX;
 433       totals->minimum = DBL_MAX;
 434     }
 435 }
 436
 437 /* Post calculations for the ONEWAY command */
 438 static void
 439 postcalc (const struct oneway *cmd)
 440 {
 441   size_t i = 0;
 442
 443   for (i = 0; i < cmd->n_vars; ++i)
 444     {
 445       struct group_proc *gp = group_proc_get (cmd->vars[i]);
 446       struct hsh_table *group_hash = gp->group_hash;
 447       struct group_statistics *totals = &gp->ugs;
 448
 449       struct hsh_iterator g;
 450       struct group_statistics *gs;
 451
 452       for (gs =  hsh_first (group_hash, &g);
 453            gs != 0;
 454            gs = hsh_next (group_hash, &g))
 455         {
 456           gs->mean = gs->sum / gs->n;
 457           gs->s_std_dev = sqrt (gs->ssq / gs->n - pow2 (gs->mean));
 458
 459           gs->std_dev = sqrt (
 460                               gs->n / (gs->n - 1) *
 461                               ( gs->ssq / gs->n - pow2 (gs->mean))
 462                               );
 463
 464           gs->se_mean = gs->std_dev / sqrt (gs->n);
 465           gs->mean_diff = gs->sum_diff / gs->n;
 466         }
 467
 468       totals->mean = totals->sum / totals->n;
 469       totals->std_dev = sqrt (
 470                               totals->n / (totals->n - 1) *
 471                               (totals->ssq / totals->n - pow2 (totals->mean))
 472                               );
 473
 474       totals->se_mean = totals->std_dev / sqrt (totals->n);
 475     }
 476 }
 477
 478 static void show_contrast_coeffs (const struct oneway *cmd);
 479 static void show_contrast_tests (const struct oneway *cmd);
 480
 481 static void
 482 output_oneway (const struct oneway *cmd, const struct dictionary *dict)
 483 {
 484   size_t i = 0;
 485
 486   /* Check the sanity of the given contrast values */
 487   struct contrasts_node *coeff_list  = NULL;
 488   ll_for_each (coeff_list, struct contrasts_node, ll, &cmd->contrast_list)
 489     {
 490       struct coeff_node *cn = NULL;
 491       double sum = 0;
 492       struct ll_list *cl = &coeff_list->coefficient_list;
 493       ++i;
 494
 495       if (ll_count (cl) != cmd->actual_number_of_groups)
 496         {
 497           msg (SW,
 498                _("Number of contrast coefficients must equal the number of groups"));
 499           coeff_list->bad_count = true;
 500           continue;
 501         }
 502
 503       ll_for_each (cn, struct coeff_node, ll, cl)
 504         sum += cn->coeff;
 505
 506       if ( sum != 0.0 )
 507         msg (SW, _("Coefficients for contrast %zu do not total zero"), i);
 508     }
 509
 510   if (cmd->stats & STATS_DESCRIPTIVES)
 511     show_descriptives (cmd, dict);
 512
 513   if (cmd->stats & STATS_HOMOGENEITY)
 514     show_homogeneity (cmd);
 515
 516   show_anova_table (cmd);
 517
 518
 519   if (ll_count (&cmd->contrast_list) > 0)
 520     {
 521       show_contrast_coeffs (cmd);
 522       show_contrast_tests (cmd);
 523     }
 524
 525
 526   /* Clean up */
 527   for (i = 0; i < cmd->n_vars; ++i )
 528     {
 529       struct hsh_table *group_hash = group_proc_get (cmd->vars[i])->group_hash;
 530
 531       hsh_destroy (group_hash);
 532     }
 533
 534   hsh_destroy (cmd->group_hash);
 535 }
 536
 537
 538 /* Show the ANOVA table */
 539 static void
 540 show_anova_table (const struct oneway *cmd)
 541 {
 542   size_t i;
 543   int n_cols =7;
 544   size_t n_rows = cmd->n_vars * 3 + 1;
 545
 546   struct tab_table *t = tab_create (n_cols, n_rows);
 547
 548   tab_headers (t, 2, 0, 1, 0);
 549
 550   tab_box (t,
 551            TAL_2, TAL_2,
 552            -1, TAL_1,
 553            0, 0,
 554            n_cols - 1, n_rows - 1);
 555
 556   tab_hline (t, TAL_2, 0, n_cols - 1, 1 );
 557   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 558   tab_vline (t, TAL_0, 1, 0, 0);
 559
 560   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 561   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 562   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 563   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 564   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 565
 566
 567   for (i = 0; i < cmd->n_vars; ++i)
 568     {
 569       struct group_statistics *totals = &group_proc_get (cmd->vars[i])->ugs;
 570       struct hsh_table *group_hash = group_proc_get (cmd->vars[i])->group_hash;
 571       struct hsh_iterator g;
 572       struct group_statistics *gs;
 573       double ssa = 0;
 574       const char *s = var_to_string (cmd->vars[i]);
 575
 576       for (gs =  hsh_first (group_hash, &g);
 577            gs != 0;
 578            gs = hsh_next (group_hash, &g))
 579         {
 580           ssa += pow2 (gs->sum) / gs->n;
 581         }
 582
 583       ssa -= pow2 (totals->sum) / totals->n;
 584
 585       tab_text (t, 0, i * 3 + 1, TAB_LEFT | TAT_TITLE, s);
 586       tab_text (t, 1, i * 3 + 1, TAB_LEFT | TAT_TITLE, _("Between Groups"));
 587       tab_text (t, 1, i * 3 + 2, TAB_LEFT | TAT_TITLE, _("Within Groups"));
 588       tab_text (t, 1, i * 3 + 3, TAB_LEFT | TAT_TITLE, _("Total"));
 589
 590       if (i > 0)
 591         tab_hline (t, TAL_1, 0, n_cols - 1, i * 3 + 1);
 592
 593       {
 594         struct group_proc *gp = group_proc_get (cmd->vars[i]);
 595         const double sst = totals->ssq - pow2 (totals->sum) / totals->n;
 596         const double df1 = gp->n_groups - 1;
 597         const double df2 = totals->n - gp->n_groups;
 598         const double msa = ssa / df1;
 599
 600         gp->mse  = (sst - ssa) / df2;
 601
 602
 603         /* Sums of Squares */
 604         tab_double (t, 2, i * 3 + 1, 0, ssa, NULL);
 605         tab_double (t, 2, i * 3 + 3, 0, sst, NULL);
 606         tab_double (t, 2, i * 3 + 2, 0, sst - ssa, NULL);
 607
 608
 609         /* Degrees of freedom */
 610         tab_fixed (t, 3, i * 3 + 1, 0, df1, 4, 0);
 611         tab_fixed (t, 3, i * 3 + 2, 0, df2, 4, 0);
 612         tab_fixed (t, 3, i * 3 + 3, 0, totals->n - 1, 4, 0);
 613
 614         /* Mean Squares */
 615         tab_double (t, 4, i * 3 + 1, TAB_RIGHT, msa, NULL);
 616         tab_double (t, 4, i * 3 + 2, TAB_RIGHT, gp->mse, NULL);
 617
 618         {
 619           const double F = msa / gp->mse ;
 620
 621           /* The F value */
 622           tab_double (t, 5, i * 3 + 1, 0,  F, NULL);
 623
 624           /* The significance */
 625           tab_double (t, 6, i * 3 + 1, 0, gsl_cdf_fdist_Q (F, df1, df2), NULL);
 626         }
 627       }
 628     }
 629
 630
 631   tab_title (t, _("ANOVA"));
 632   tab_submit (t);
 633 }
 634
 635
 636 /* Show the descriptives table */
 637 static void
 638 show_descriptives (const struct oneway *cmd, const struct dictionary *dict)
 639 {
 640   size_t v;
 641   int n_cols = 10;
 642   struct tab_table *t;
 643   int row;
 644
 645   const double confidence = 0.95;
 646   const double q = (1.0 - confidence) / 2.0;
 647
 648   const struct variable *wv = dict_get_weight (dict);
 649   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : & F_8_0;
 650
 651   int n_rows = 2;
 652
 653   for ( v = 0; v < cmd->n_vars; ++v )
 654     n_rows += group_proc_get (cmd->vars[v])->n_groups + 1;
 655
 656   t = tab_create (n_cols, n_rows);
 657   tab_headers (t, 2, 0, 2, 0);
 658
 659
 660   /* Put a frame around the entire box, and vertical lines inside */
 661   tab_box (t,
 662            TAL_2, TAL_2,
 663            -1, TAL_1,
 664            0, 0,
 665            n_cols - 1, n_rows - 1);
 666
 667   /* Underline headers */
 668   tab_hline (t, TAL_2, 0, n_cols - 1, 2);
 669   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 670
 671   tab_text (t, 2, 1, TAB_CENTER | TAT_TITLE, _("N"));
 672   tab_text (t, 3, 1, TAB_CENTER | TAT_TITLE, _("Mean"));
 673   tab_text (t, 4, 1, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 674   tab_text (t, 5, 1, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 675
 676
 677   tab_vline (t, TAL_0, 7, 0, 0);
 678   tab_hline (t, TAL_1, 6, 7, 1);
 679   tab_joint_text_format (t, 6, 0, 7, 0, TAB_CENTER | TAT_TITLE,
 680                          _("%g%% Confidence Interval for Mean"),
 681                          confidence*100.0);
 682
 683   tab_text (t, 6, 1, TAB_CENTER | TAT_TITLE, _("Lower Bound"));
 684   tab_text (t, 7, 1, TAB_CENTER | TAT_TITLE, _("Upper Bound"));
 685
 686   tab_text (t, 8, 1, TAB_CENTER | TAT_TITLE, _("Minimum"));
 687   tab_text (t, 9, 1, TAB_CENTER | TAT_TITLE, _("Maximum"));
 688
 689
 690   tab_title (t, _("Descriptives"));
 691
 692
 693   row = 2;
 694   for (v = 0; v < cmd->n_vars; ++v)
 695     {
 696       double T;
 697       double std_error;
 698
 699       struct group_proc *gp = group_proc_get (cmd->vars[v]);
 700
 701       struct group_statistics *gs;
 702       struct group_statistics *totals = &gp->ugs;
 703
 704       const char *s = var_to_string (cmd->vars[v]);
 705       const struct fmt_spec *fmt = var_get_print_format (cmd->vars[v]);
 706
 707       struct group_statistics *const *gs_array =
 708         (struct group_statistics *const *) hsh_sort (gp->group_hash);
 709       int count = 0;
 710
 711       tab_text (t, 0, row, TAB_LEFT | TAT_TITLE, s);
 712       if ( v > 0)
 713         tab_hline (t, TAL_1, 0, n_cols - 1, row);
 714
 715       for (count = 0; count < hsh_count (gp->group_hash); ++count)
 716         {
 717           struct string vstr;
 718           ds_init_empty (&vstr);
 719           gs = gs_array[count];
 720
 721           var_append_value_name (cmd->indep_var, &gs->id, &vstr);
 722
 723           tab_text (t, 1, row + count,
 724                     TAB_LEFT | TAT_TITLE,
 725                     ds_cstr (&vstr));
 726
 727           ds_destroy (&vstr);
 728
 729           /* Now fill in the numbers ... */
 730
 731           tab_fixed (t, 2, row + count, 0, gs->n, 8, 0);
 732
 733           tab_double (t, 3, row + count, 0, gs->mean, NULL);
 734
 735           tab_double (t, 4, row + count, 0, gs->std_dev, NULL);
 736
 737           std_error = gs->std_dev / sqrt (gs->n) ;
 738           tab_double (t, 5, row + count, 0,
 739                       std_error, NULL);
 740
 741           /* Now the confidence interval */
 742
 743           T = gsl_cdf_tdist_Qinv (q, gs->n - 1);
 744
 745           tab_double (t, 6, row + count, 0,
 746                       gs->mean - T * std_error, NULL);
 747
 748           tab_double (t, 7, row + count, 0,
 749                       gs->mean + T * std_error, NULL);
 750
 751           /* Min and Max */
 752
 753           tab_double (t, 8, row + count, 0,  gs->minimum, fmt);
 754           tab_double (t, 9, row + count, 0,  gs->maximum, fmt);
 755         }
 756
 757       tab_text (t, 1, row + count,
 758                 TAB_LEFT | TAT_TITLE, _("Total"));
 759
 760       tab_double (t, 2, row + count, 0, totals->n, wfmt);
 761
 762       tab_double (t, 3, row + count, 0, totals->mean, NULL);
 763
 764       tab_double (t, 4, row + count, 0, totals->std_dev, NULL);
 765
 766       std_error = totals->std_dev / sqrt (totals->n) ;
 767
 768       tab_double (t, 5, row + count, 0, std_error, NULL);
 769
 770       /* Now the confidence interval */
 771
 772       T = gsl_cdf_tdist_Qinv (q, totals->n - 1);
 773
 774       tab_double (t, 6, row + count, 0,
 775                   totals->mean - T * std_error, NULL);
 776
 777       tab_double (t, 7, row + count, 0,
 778                   totals->mean + T * std_error, NULL);
 779
 780       /* Min and Max */
 781
 782       tab_double (t, 8, row + count, 0,  totals->minimum, fmt);
 783       tab_double (t, 9, row + count, 0,  totals->maximum, fmt);
 784
 785       row += gp->n_groups + 1;
 786     }
 787
 788   tab_submit (t);
 789 }
 790
 791 /* Show the homogeneity table */
 792 static void
 793 show_homogeneity (const struct oneway *cmd)
 794 {
 795   size_t v;
 796   int n_cols = 5;
 797   size_t n_rows = cmd->n_vars + 1;
 798
 799   struct tab_table *t;
 800
 801
 802   t = tab_create (n_cols, n_rows);
 803   tab_headers (t, 1, 0, 1, 0);
 804
 805
 806   /* Put a frame around the entire box, and vertical lines inside */
 807   tab_box (t,
 808            TAL_2, TAL_2,
 809            -1, TAL_1,
 810            0, 0,
 811            n_cols - 1, n_rows - 1);
 812
 813
 814   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 815   tab_vline (t, TAL_2, 1, 0, n_rows - 1);
 816
 817
 818   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Levene Statistic"));
 819   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df1"));
 820   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df2"));
 821   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 822
 823   tab_title (t, _("Test of Homogeneity of Variances"));
 824
 825   for (v = 0; v < cmd->n_vars; ++v)
 826     {
 827       double F;
 828       const struct variable *var = cmd->vars[v];
 829       const struct group_proc *gp = group_proc_get (cmd->vars[v]);
 830       const char *s = var_to_string (var);
 831       const struct group_statistics *totals = &gp->ugs;
 832
 833       const double df1 = gp->n_groups - 1;
 834       const double df2 = totals->n - gp->n_groups;
 835
 836       tab_text (t, 0, v + 1, TAB_LEFT | TAT_TITLE, s);
 837
 838       F = gp->levene;
 839       tab_double (t, 1, v + 1, TAB_RIGHT, F, NULL);
 840       tab_fixed (t, 2, v + 1, TAB_RIGHT, df1, 8, 0);
 841       tab_fixed (t, 3, v + 1, TAB_RIGHT, df2, 8, 0);
 842
 843       /* Now the significance */
 844       tab_double (t, 4, v + 1, TAB_RIGHT,gsl_cdf_fdist_Q (F, df1, df2), NULL);
 845     }
 846
 847   tab_submit (t);
 848 }
 849
 850
 851 /* Show the contrast coefficients table */
 852 static void
 853 show_contrast_coeffs (const struct oneway *cmd)
 854 {
 855   int c_num = 0;
 856   struct ll *cli;
 857
 858   int n_contrasts = ll_count (&cmd->contrast_list);
 859   int n_cols = 2 + cmd->actual_number_of_groups;
 860   int n_rows = 2 + n_contrasts;
 861
 862   void *const *group_values;
 863
 864   struct tab_table *t;
 865
 866   t = tab_create (n_cols, n_rows);
 867   tab_headers (t, 2, 0, 2, 0);
 868
 869   /* Put a frame around the entire box, and vertical lines inside */
 870   tab_box (t,
 871            TAL_2, TAL_2,
 872            -1, TAL_1,
 873            0, 0,
 874            n_cols - 1, n_rows - 1);
 875
 876   tab_box (t,
 877            -1, -1,
 878            TAL_0, TAL_0,
 879            2, 0,
 880            n_cols - 1, 0);
 881
 882   tab_box (t,
 883            -1, -1,
 884            TAL_0, TAL_0,
 885            0, 0,
 886            1, 1);
 887
 888   tab_hline (t, TAL_1, 2, n_cols - 1, 1);
 889   tab_hline (t, TAL_2, 0, n_cols - 1, 2);
 890
 891   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 892
 893   tab_title (t, _("Contrast Coefficients"));
 894
 895   tab_text (t,  0, 2, TAB_LEFT | TAT_TITLE, _("Contrast"));
 896
 897
 898   tab_joint_text (t, 2, 0, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
 899                   var_to_string (cmd->indep_var));
 900
 901   group_values = hsh_sort (cmd->group_hash);
 902
 903   for ( cli = ll_head (&cmd->contrast_list);
 904         cli != ll_null (&cmd->contrast_list);
 905         cli = ll_next (cli))
 906     {
 907       int count = 0;
 908       struct contrasts_node *cn = ll_data (cli, struct contrasts_node, ll);
 909       struct ll *coeffi = ll_head (&cn->coefficient_list);
 910
 911       tab_text_format (t, 1, c_num + 2, TAB_CENTER, "%d", c_num + 1);
 912
 913       for (count = 0;
 914            count < hsh_count (cmd->group_hash) && coeffi != ll_null (&cn->coefficient_list);
 915            ++count)
 916         {
 917           double *group_value_p;
 918           union value group_value;
 919           struct string vstr;
 920
 921           ds_init_empty (&vstr);
 922
 923           group_value_p = group_values[count];
 924           group_value.f = *group_value_p;
 925           var_append_value_name (cmd->indep_var, &group_value, &vstr);
 926
 927           tab_text (t, count + 2, 1, TAB_CENTER | TAT_TITLE,
 928                     ds_cstr (&vstr));
 929
 930           ds_destroy (&vstr);
 931
 932           if (cn->bad_count)
 933             tab_text (t, count + 2, c_num + 2, TAB_RIGHT, "?" );
 934           else
 935             {
 936               struct coeff_node *coeffn = ll_data (coeffi, struct coeff_node, ll);
 937
 938               tab_text_format (t, count + 2, c_num + 2, TAB_RIGHT, "%g", coeffn->coeff);
 939             }
 940
 941           coeffi = ll_next (coeffi);
 942         }
 943       ++c_num;
 944     }
 945
 946   tab_submit (t);
 947 }
 948
 949
 950 /* Show the results of the contrast tests */
 951 static void
 952 show_contrast_tests (const struct oneway *cmd)
 953 {
 954   int n_contrasts = ll_count (&cmd->contrast_list);
 955   size_t v;
 956   int n_cols = 8;
 957   size_t n_rows = 1 + cmd->n_vars * 2 * n_contrasts;
 958
 959   struct tab_table *t;
 960
 961   t = tab_create (n_cols, n_rows);
 962   tab_headers (t, 3, 0, 1, 0);
 963
 964   /* Put a frame around the entire box, and vertical lines inside */
 965   tab_box (t,
 966            TAL_2, TAL_2,
 967            -1, TAL_1,
 968            0, 0,
 969            n_cols - 1, n_rows - 1);
 970
 971   tab_box (t,
 972            -1, -1,
 973            TAL_0, TAL_0,
 974            0, 0,
 975            2, 0);
 976
 977   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 978   tab_vline (t, TAL_2, 3, 0, n_rows - 1);
 979
 980
 981   tab_title (t, _("Contrast Tests"));
 982
 983   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Contrast"));
 984   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Value of Contrast"));
 985   tab_text (t,  4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 986   tab_text (t,  5, 0, TAB_CENTER | TAT_TITLE, _("t"));
 987   tab_text (t,  6, 0, TAB_CENTER | TAT_TITLE, _("df"));
 988   tab_text (t,  7, 0, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
 989
 990   for (v = 0; v < cmd->n_vars; ++v)
 991     {
 992       struct ll *cli;
 993       int i = 0;
 994       int lines_per_variable = 2 * n_contrasts;
 995
 996       tab_text (t,  0, (v * lines_per_variable) + 1, TAB_LEFT | TAT_TITLE,
 997                 var_to_string (cmd->vars[v]));
 998
 999       for ( cli = ll_head (&cmd->contrast_list);
1000             cli != ll_null (&cmd->contrast_list);
1001             ++i, cli = ll_next (cli))
1002         {
1003           struct contrasts_node *cn = ll_data (cli, struct contrasts_node, ll);
1004           struct ll *coeffi = ll_head (&cn->coefficient_list);
1005           int ci;
1006           double contrast_value = 0.0;
1007           double coef_msq = 0.0;
1008           struct group_proc *grp_data = group_proc_get (cmd->vars[v]);
1009           struct hsh_table *group_hash = grp_data->group_hash;
1010
1011           void *const *group_stat_array;
1012
1013           double T;
1014           double std_error_contrast;
1015           double df;
1016           double sec_vneq = 0.0;
1017
1018           /* Note: The calculation of the degrees of freedom in the
1019              "variances not equal" case is painfull!!
1020              The following formula may help to understand it:
1021              \frac{\left (\sum_{i=1}^k{c_i^2\frac{s_i^2}{n_i}}\right)^2}
1022              {
1023              \sum_{i=1}^k\left (
1024              \frac{\left (c_i^2\frac{s_i^2}{n_i}\right)^2}  {n_i-1}
1025              \right)
1026              }
1027           */
1028
1029           double df_denominator = 0.0;
1030           double df_numerator = 0.0;
1031           if ( i == 0 )
1032             {
1033               tab_text (t,  1, (v * lines_per_variable) + i + 1,
1034                         TAB_LEFT | TAT_TITLE,
1035                         _("Assume equal variances"));
1036
1037               tab_text (t,  1, (v * lines_per_variable) + i + 1 + n_contrasts,
1038                         TAB_LEFT | TAT_TITLE,
1039                         _("Does not assume equal"));
1040             }
1041
1042           tab_text_format (t,  2, (v * lines_per_variable) + i + 1,
1043                            TAB_CENTER | TAT_TITLE, "%d", i + 1);
1044
1045
1046           tab_text_format (t,  2,
1047                            (v * lines_per_variable) + i + 1 + n_contrasts,
1048                            TAB_CENTER | TAT_TITLE, "%d", i + 1);
1049
1050           if (cn->bad_count)
1051             continue;
1052
1053           group_stat_array = hsh_sort (group_hash);
1054
1055           for (ci = 0;
1056                coeffi != ll_null (&cn->coefficient_list) &&
1057                  ci < hsh_count (group_hash);
1058                ++ci, coeffi = ll_next (coeffi))
1059             {
1060               struct coeff_node *cn = ll_data (coeffi, struct coeff_node, ll);
1061               const double coef = cn->coeff;
1062               struct group_statistics *gs = group_stat_array[ci];
1063
1064               const double winv = pow2 (gs->std_dev) / gs->n;
1065
1066               contrast_value += coef * gs->mean;
1067
1068               coef_msq += (coef * coef) / gs->n;
1069
1070               sec_vneq += (coef * coef) * pow2 (gs->std_dev) /gs->n;
1071
1072               df_numerator += (coef * coef) * winv;
1073               df_denominator += pow2((coef * coef) * winv) / (gs->n - 1);
1074             }
1075
1076           sec_vneq = sqrt (sec_vneq);
1077
1078           df_numerator = pow2 (df_numerator);
1079
1080           tab_double (t,  3, (v * lines_per_variable) + i + 1,
1081                       TAB_RIGHT, contrast_value, NULL);
1082
1083           tab_double (t,  3, (v * lines_per_variable) + i + 1 +
1084                       n_contrasts,
1085                       TAB_RIGHT, contrast_value, NULL);
1086
1087           std_error_contrast = sqrt (grp_data->mse * coef_msq);
1088
1089           /* Std. Error */
1090           tab_double (t,  4, (v * lines_per_variable) + i + 1,
1091                       TAB_RIGHT, std_error_contrast,
1092                       NULL);
1093
1094           T = fabs (contrast_value / std_error_contrast);
1095
1096           /* T Statistic */
1097
1098           tab_double (t,  5, (v * lines_per_variable) + i + 1,
1099                       TAB_RIGHT, T,
1100                       NULL);
1101
1102           df = grp_data->ugs.n - grp_data->n_groups;
1103
1104           /* Degrees of Freedom */
1105           tab_fixed (t,  6, (v * lines_per_variable) + i + 1,
1106                      TAB_RIGHT,  df,
1107                      8, 0);
1108
1109
1110           /* Significance TWO TAILED !!*/
1111           tab_double (t,  7, (v * lines_per_variable) + i + 1,
1112                       TAB_RIGHT,  2 * gsl_cdf_tdist_Q (T, df),
1113                       NULL);
1114
1115           /* Now for the Variances NOT Equal case */
1116
1117           /* Std. Error */
1118           tab_double (t,  4,
1119                       (v * lines_per_variable) + i + 1 + n_contrasts,
1120                       TAB_RIGHT, sec_vneq,
1121                       NULL);
1122
1123           T = contrast_value / sec_vneq;
1124           tab_double (t,  5,
1125                       (v * lines_per_variable) + i + 1 + n_contrasts,
1126                       TAB_RIGHT, T,
1127                       NULL);
1128
1129           df = df_numerator / df_denominator;
1130
1131           tab_double (t,  6,
1132                       (v * lines_per_variable) + i + 1 + n_contrasts,
1133                       TAB_RIGHT, df,
1134                       NULL);
1135
1136           /* The Significance */
1137           tab_double (t, 7, (v * lines_per_variable) + i + 1 + n_contrasts,
1138                       TAB_RIGHT,  2 * gsl_cdf_tdist_Q (T,df),
1139                       NULL);
1140         }
1141
1142       if ( v > 0 )
1143         tab_hline (t, TAL_1, 0, n_cols - 1, (v * lines_per_variable) + 1);
1144     }
1145
1146   tab_submit (t);
1147 }
1148
1149