src/language/stats/kruskal-wallis.c

   1 /* Pspp - a program for statistical analysis.
   2    Copyright (C) 2010, 2011 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
  18 #include <config.h>
  19
  20 #include "kruskal-wallis.h"
  21
  22 #include <gsl/gsl_cdf.h>
  23 #include <math.h>
  24
  25 #include "data/casereader.h"
  26 #include "data/casewriter.h"
  27 #include "data/dataset.h"
  28 #include "data/dictionary.h"
  29 #include "data/format.h"
  30 #include "data/subcase.h"
  31 #include "data/variable.h"
  32 #include "libpspp/assertion.h"
  33 #include "libpspp/hmap.h"
  34 #include "libpspp/message.h"
  35 #include "libpspp/misc.h"
  36 #include "math/sort.h"
  37 #include "output/tab.h"
  38
  39 #include "gl/minmax.h"
  40 #include "gl/xalloc.h"
  41
  42
  43 /* Returns true iff the independent variable lies in the range [nst->val1, nst->val2] */
  44 static bool
  45 include_func (const struct ccase *c, void *aux)
  46 {
  47   const struct n_sample_test *nst = aux;
  48
  49   if (0 < value_compare_3way (&nst->val1, case_data (c, nst->indep_var), var_get_width (nst->indep_var)))
  50     return false;
  51
  52   if (0 > value_compare_3way (&nst->val2, case_data (c, nst->indep_var), var_get_width (nst->indep_var)))
  53     return false;
  54
  55   return true;
  56 }
  57
  58
  59 struct rank_entry
  60 {
  61   struct hmap_node node;
  62   union value group;
  63
  64   double sum_of_ranks;
  65   double n;
  66 };
  67
  68 /* Return the entry with the key GROUP or null if there is no such entry */
  69 static struct rank_entry *
  70 find_rank_entry (const struct hmap *map, const union value *group, size_t width)
  71 {
  72   struct rank_entry *re = NULL;
  73   size_t hash  = value_hash (group, width, 0);
  74
  75   HMAP_FOR_EACH_WITH_HASH (re, struct rank_entry, node, hash, map)
  76     {
  77       if (0 == value_compare_3way (group, &re->group, width))
  78         return re;
  79     }
  80
  81   return re;
  82 }
  83
  84 /* Calculates the adjustment necessary for tie compensation */
  85 static void
  86 distinct_callback (double v UNUSED, casenumber t, double w UNUSED, void *aux)
  87 {
  88   double *tiebreaker = aux;
  89
  90   *tiebreaker += pow3 (t) - t;
  91 }
  92
  93
  94 struct kw
  95 {
  96   struct hmap map;
  97   double h;
  98 };
  99
 100 static void show_ranks_box (const struct n_sample_test *nst, const struct kw *kw, int n_groups);
 101 static void show_sig_box (const struct n_sample_test *nst, const struct kw *kw);
 102
 103 void
 104 kruskal_wallis_execute (const struct dataset *ds,
 105                         struct casereader *input,
 106                         enum mv_class exclude,
 107                         const struct npar_test *test,
 108                         bool exact UNUSED,
 109                         double timer UNUSED)
 110 {
 111   int i;
 112   struct ccase *c;
 113   bool warn = true;
 114   const struct dictionary *dict = dataset_dict (ds);
 115   const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
 116   const struct caseproto *proto ;
 117   size_t rank_idx ;
 118
 119   int total_n_groups = 0.0;
 120
 121   struct kw *kw = xcalloc (nst->n_vars, sizeof *kw);
 122
 123   /* If the independent variable is missing, then we ignore the case */
 124   input = casereader_create_filter_missing (input,
 125                                             &nst->indep_var, 1,
 126                                             exclude,
 127                                             NULL, NULL);
 128
 129   input = casereader_create_filter_weight (input, dict, &warn, NULL);
 130
 131   /* Remove all those cases which are outside the range (val1, val2) */
 132   input = casereader_create_filter_func (input, include_func, NULL,
 133         CONST_CAST (struct n_sample_test *, nst), NULL);
 134
 135   proto = casereader_get_proto (input);
 136   rank_idx = caseproto_get_n_widths (proto);
 137
 138   /* Rank cases by the v value */
 139   for (i = 0; i < nst->n_vars; ++i)
 140     {
 141       double tiebreaker = 0.0;
 142       bool warn = true;
 143       enum rank_error rerr = 0;
 144       struct casereader *rr;
 145       struct casereader *r = casereader_clone (input);
 146
 147       r = sort_execute_1var (r, nst->vars[i]);
 148
 149       /* Ignore missings in the test variable */
 150       r = casereader_create_filter_missing (r, &nst->vars[i], 1,
 151                                             exclude,
 152                                             NULL, NULL);
 153
 154       rr = casereader_create_append_rank (r,
 155                                           nst->vars[i],
 156                                           dict_get_weight (dict),
 157                                           &rerr,
 158                                           distinct_callback, &tiebreaker);
 159
 160       hmap_init (&kw[i].map);
 161       for (; (c = casereader_read (rr)); case_unref (c))
 162         {
 163           const union value *group = case_data (c, nst->indep_var);
 164           const size_t group_var_width = var_get_width (nst->indep_var);
 165           struct rank_entry *rank = find_rank_entry (&kw[i].map, group, group_var_width);
 166
 167           if ( NULL == rank)
 168             {
 169               rank = xzalloc (sizeof *rank);
 170               value_clone (&rank->group, group, group_var_width);
 171
 172               hmap_insert (&kw[i].map, &rank->node,
 173                            value_hash (&rank->group, group_var_width, 0));
 174             }
 175
 176           rank->sum_of_ranks += case_data_idx (c, rank_idx)->f;
 177           rank->n += dict_get_case_weight (dict, c, &warn);
 178
 179           /* If this assertion fires, then either the data wasn't sorted or some other
 180              problem occured */
 181           assert (rerr == 0);
 182         }
 183
 184       casereader_destroy (rr);
 185
 186       /* Calculate the value of h */
 187       {
 188         struct rank_entry *mre;
 189         double n = 0.0;
 190
 191         HMAP_FOR_EACH (mre, struct rank_entry, node, &kw[i].map)
 192           {
 193             kw[i].h += pow2 (mre->sum_of_ranks) / mre->n;
 194             n += mre->n;
 195
 196             total_n_groups ++;
 197           }
 198         kw[i].h *= 12 / (n * ( n + 1));
 199         kw[i].h -= 3 * (n + 1) ;
 200
 201         kw[i].h /= 1 - tiebreaker/ (pow3 (n) - n);
 202       }
 203     }
 204
 205   casereader_destroy (input);
 206
 207   show_ranks_box (nst, kw, total_n_groups);
 208   show_sig_box (nst, kw);
 209
 210   /* Cleanup allocated memory */
 211   for (i = 0 ; i < nst->n_vars; ++i)
 212     {
 213       struct rank_entry *mre, *next;
 214       HMAP_FOR_EACH_SAFE (mre, next, struct rank_entry, node, &kw[i].map)
 215         {
 216           hmap_delete (&kw[i].map, &mre->node);
 217           free (mre);
 218         }
 219       hmap_destroy (&kw[i].map);
 220     }
 221
 222   free (kw);
 223 }
 224
 225 \f
 226 #include "gettext.h"
 227 #define _(msgid) gettext (msgid)
 228
 229
 230 static void
 231 show_ranks_box (const struct n_sample_test *nst, const struct kw *kw, int n_groups)
 232 {
 233   int row;
 234   int i;
 235   const int row_headers = 2;
 236   const int column_headers = 1;
 237   struct tab_table *table =
 238     tab_create (row_headers + 2, column_headers + n_groups + nst->n_vars);
 239
 240   tab_headers (table, row_headers, 0, column_headers, 0);
 241
 242   tab_title (table, _("Ranks"));
 243
 244   /* Vertical lines inside the box */
 245   tab_box (table, 1, 0, -1, TAL_1,
 246            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
 247
 248   /* Box around the table */
 249   tab_box (table, TAL_2, TAL_2, -1, -1,
 250            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 251
 252   tab_text (table, 1, 0, TAT_TITLE,
 253             var_to_string (nst->indep_var)
 254             );
 255
 256   tab_text (table, 3, 0, 0, _("Mean Rank"));
 257   tab_text (table, 2, 0, 0, _("N"));
 258
 259   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 260   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 261
 262
 263   row = column_headers;
 264   for (i = 0 ; i < nst->n_vars ; ++i)
 265     {
 266       int tot = 0;
 267       const struct rank_entry *re;
 268
 269       if (i > 0)
 270         tab_hline (table, TAL_1, 0, tab_nc (table) -1, row);
 271
 272       tab_text (table,  0, row,
 273                 TAT_TITLE, var_to_string (nst->vars[i]));
 274
 275       HMAP_FOR_EACH (re, const struct rank_entry, node, &kw[i].map)
 276         {
 277           struct string str;
 278           ds_init_empty (&str);
 279
 280           var_append_value_name (nst->indep_var, &re->group, &str);
 281
 282           tab_text   (table, 1, row, TAB_LEFT, ds_cstr (&str));
 283           tab_double (table, 2, row, TAB_LEFT, re->n, &F_8_0);
 284           tab_double (table, 3, row, TAB_LEFT, re->sum_of_ranks / re->n, 0);
 285
 286           tot += re->n;
 287           row++;
 288           ds_destroy (&str);
 289         }
 290       tab_double (table, 2, row, TAB_LEFT,
 291                   tot, &F_8_0);
 292       tab_text (table, 1, row++, TAB_LEFT, _("Total"));
 293     }
 294
 295   tab_submit (table);
 296 }
 297
 298
 299 static void
 300 show_sig_box (const struct n_sample_test *nst, const struct kw *kw)
 301 {
 302   int i;
 303   const int row_headers = 1;
 304   const int column_headers = 1;
 305   struct tab_table *table =
 306     tab_create (row_headers + nst->n_vars * 2, column_headers + 3);
 307
 308   tab_headers (table, row_headers, 0, column_headers, 0);
 309
 310   tab_title (table, _("Test Statistics"));
 311
 312   tab_text (table,  0, column_headers,
 313             TAT_TITLE | TAB_LEFT , _("Chi-Square"));
 314
 315   tab_text (table,  0, 1 + column_headers,
 316             TAT_TITLE | TAB_LEFT, _("df"));
 317
 318   tab_text (table,  0, 2 + column_headers,
 319             TAT_TITLE | TAB_LEFT, _("Asymp. Sig."));
 320
 321   /* Box around the table */
 322   tab_box (table, TAL_2, TAL_2, -1, -1,
 323            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 324
 325
 326   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 327   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 328
 329   for (i = 0 ; i < nst->n_vars; ++i)
 330     {
 331       const double df = hmap_count (&kw[i].map) - 1;
 332       tab_text (table, column_headers + 1 + i, 0, TAT_TITLE,
 333                 var_to_string (nst->vars[i])
 334                 );
 335
 336       tab_double (table, column_headers + 1 + i, 1, 0,
 337                   kw[i].h, 0);
 338
 339       tab_double (table, column_headers + 1 + i, 2, 0,
 340                   df, &F_8_0);
 341
 342       tab_double (table, column_headers + 1 + i, 3, 0,
 343                   gsl_cdf_chisq_Q (kw[i].h, df),
 344                   0);
 345     }
 346
 347   tab_submit (table);
 348 }