src/language/stats/wilcoxon.c

   1 /* Pspp - a program for statistical analysis.
   2    Copyright (C) 2008, 2009 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
  19 #include <config.h>
  20 #include "wilcoxon.h"
  21 #include <data/variable.h>
  22 #include <data/casereader.h>
  23 #include <data/casewriter.h>
  24 #include <data/subcase.h>
  25 #include <math/sort.h>
  26 #include <libpspp/message.h>
  27 #include <xalloc.h>
  28 #include <output/table.h>
  29 #include <data/procedure.h>
  30 #include <data/dictionary.h>
  31 #include <math/wilcoxon-sig.h>
  32 #include <gsl/gsl_cdf.h>
  33 #include <unistd.h>
  34 #include <signal.h>
  35 #include <libpspp/assertion.h>
  36
  37 static double
  38 append_difference (const struct ccase *c, casenumber n UNUSED, void *aux)
  39 {
  40   const variable_pair *vp = aux;
  41
  42   return case_data (c, (*vp)[0])->f - case_data (c, (*vp)[1])->f;
  43 }
  44
  45 static void show_ranks_box (const struct wilcoxon_state *,
  46                             const struct two_sample_test *);
  47
  48 static void show_tests_box (const struct wilcoxon_state *,
  49                             const struct two_sample_test *,
  50                             bool exact, double timer);
  51
  52
  53
  54 static void
  55 distinct_callback (double v UNUSED, casenumber n, double w UNUSED, void *aux)
  56 {
  57   struct wilcoxon_state *ws = aux;
  58
  59   ws->tiebreaker += pow3 (n) - n;
  60 }
  61
  62 #define WEIGHT_IDX 2
  63
  64 void
  65 wilcoxon_execute (const struct dataset *ds,
  66                   struct casereader *input,
  67                   enum mv_class exclude,
  68                   const struct npar_test *test,
  69                   bool exact,
  70                   double timer)
  71 {
  72   int i;
  73   bool warn = true;
  74   const struct dictionary *dict = dataset_dict (ds);
  75   const struct two_sample_test *t2s = (struct two_sample_test *) test;
  76
  77   struct wilcoxon_state *ws = xcalloc (sizeof (*ws), t2s->n_pairs);
  78   const struct variable *weight = dict_get_weight (dict);
  79   struct variable *weightx = var_create_internal (WEIGHT_IDX);
  80
  81   input =
  82     casereader_create_filter_weight (input, dict, &warn, NULL);
  83
  84   for (i = 0 ; i < t2s->n_pairs; ++i )
  85     {
  86       struct casereader *r = casereader_clone (input);
  87       struct casewriter *writer;
  88       struct ccase *c;
  89       struct subcase ordering;
  90       variable_pair *vp = &t2s->pairs[i];
  91
  92       const int reader_width = weight ? 3 : 2;
  93
  94       ws[i].sign = var_create_internal (0);
  95       ws[i].absdiff = var_create_internal (1);
  96
  97       r = casereader_create_filter_missing (r, *vp, 2,
  98                                             exclude,
  99                                             NULL, NULL);
 100
 101       subcase_init_var (&ordering, ws[i].absdiff, SC_ASCEND);
 102       writer = sort_create_writer (&ordering, reader_width);
 103       subcase_destroy (&ordering);
 104
 105       for (; (c = casereader_read (r)) != NULL; case_unref (c))
 106         {
 107           struct ccase *output = case_create (reader_width);
 108           double d = append_difference (c, 0, vp);
 109
 110           if (d > 0)
 111             {
 112               case_data_rw (output, ws[i].sign)->f = 1.0;
 113
 114             }
 115           else if (d < 0)
 116             {
 117               case_data_rw (output, ws[i].sign)->f = -1.0;
 118             }
 119           else
 120             {
 121               double w = 1.0;
 122               if (weight)
 123                 w = case_data (c, weight)->f;
 124
 125               /* Central point values should be dropped */
 126               ws[i].n_zeros += w;
 127               case_unref (output);
 128               continue;
 129             }
 130
 131           case_data_rw (output, ws[i].absdiff)->f = fabs (d);
 132
 133           if (weight)
 134            case_data_rw (output, weightx)->f = case_data (c, weight)->f;
 135
 136           casewriter_write (writer, output);
 137         }
 138       casereader_destroy (r);
 139       ws[i].reader = casewriter_make_reader (writer);
 140     }
 141
 142   for (i = 0 ; i < t2s->n_pairs; ++i )
 143     {
 144       struct casereader *rr ;
 145       struct ccase *c;
 146       enum rank_error err = 0;
 147
 148       rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff,
 149                                           weight ? weightx : NULL, &err,
 150                                           distinct_callback, &ws[i]
 151                                           );
 152
 153       for (; (c = casereader_read (rr)) != NULL; case_unref (c))
 154         {
 155           double sign = case_data (c, ws[i].sign)->f;
 156           double rank = case_data_idx (c, weight ? 3 : 2)->f;
 157           double w = 1.0;
 158           if (weight)
 159             w = case_data (c, weightx)->f;
 160
 161           if ( sign > 0 )
 162             {
 163               ws[i].positives.sum += rank * w;
 164               ws[i].positives.n += w;
 165             }
 166           else if (sign < 0)
 167             {
 168               ws[i].negatives.sum += rank * w;
 169               ws[i].negatives.n += w;
 170             }
 171           else
 172             NOT_REACHED ();
 173         }
 174
 175       casereader_destroy (rr);
 176     }
 177
 178   casereader_destroy (input);
 179
 180   var_destroy (weightx);
 181
 182   show_ranks_box (ws, t2s);
 183   show_tests_box (ws, t2s, exact, timer);
 184
 185   for (i = 0 ; i < t2s->n_pairs; ++i )
 186     {
 187       var_destroy (ws[i].sign);
 188       var_destroy (ws[i].absdiff);
 189     }
 190
 191   free (ws);
 192 }
 193
 194
 195 \f
 196
 197 #include "gettext.h"
 198 #define _(msgid) gettext (msgid)
 199
 200 static void
 201 show_ranks_box (const struct wilcoxon_state *ws, const struct two_sample_test *t2s)
 202 {
 203   size_t i;
 204   struct tab_table *table = tab_create (5, 1 + 4 * t2s->n_pairs, 0);
 205
 206   tab_dim (table, tab_natural_dimensions);
 207
 208   tab_title (table, _("Ranks"));
 209
 210   tab_headers (table, 2, 0, 1, 0);
 211
 212   /* Vertical lines inside the box */
 213   tab_box (table, 0, 0, -1, TAL_1,
 214            1, 0, table->nc - 1, tab_nr (table) - 1 );
 215
 216   /* Box around entire table */
 217   tab_box (table, TAL_2, TAL_2, -1, -1,
 218            0, 0, table->nc - 1, tab_nr (table) - 1 );
 219
 220
 221   tab_text (table,  2, 0,  TAB_CENTER, _("N"));
 222   tab_text (table,  3, 0,  TAB_CENTER, _("Mean Rank"));
 223   tab_text (table,  4, 0,  TAB_CENTER, _("Sum of Ranks"));
 224
 225
 226   for (i = 0 ; i < t2s->n_pairs; ++i)
 227     {
 228       variable_pair *vp = &t2s->pairs[i];
 229
 230       struct string pair_name;
 231       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
 232       ds_put_cstr (&pair_name, " - ");
 233       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
 234
 235       tab_text (table, 1, 1 + i * 4, TAB_LEFT, _("Negative Ranks"));
 236       tab_text (table, 1, 2 + i * 4, TAB_LEFT, _("Positive Ranks"));
 237       tab_text (table, 1, 3 + i * 4, TAB_LEFT, _("Ties"));
 238       tab_text (table, 1, 4 + i * 4, TAB_LEFT, _("Total"));
 239
 240       tab_hline (table, TAL_1, 0, table->nc - 1, 1 + i * 4);
 241
 242
 243       tab_text (table, 0, 1 + i * 4, TAB_LEFT, ds_cstr (&pair_name));
 244       ds_destroy (&pair_name);
 245
 246
 247       /* N */
 248       tab_float (table, 2, 1 + i * 4, TAB_RIGHT, ws[i].negatives.n, 8, 0);
 249       tab_float (table, 2, 2 + i * 4, TAB_RIGHT, ws[i].positives.n, 8, 0);
 250       tab_float (table, 2, 3 + i * 4, TAB_RIGHT, ws[i].n_zeros, 8, 0);
 251
 252       tab_float (table, 2, 4 + i * 4, TAB_RIGHT,
 253                  ws[i].n_zeros + ws[i].positives.n + ws[i].negatives.n, 8, 0);
 254
 255       /* Sums */
 256       tab_float (table, 4, 1 + i * 4, TAB_RIGHT, ws[i].negatives.sum, 8, 2);
 257       tab_float (table, 4, 2 + i * 4, TAB_RIGHT, ws[i].positives.sum, 8, 2);
 258
 259
 260       /* Means */
 261       tab_float (table, 3, 1 + i * 4, TAB_RIGHT,
 262                  ws[i].negatives.sum / (double) ws[i].negatives.n, 8, 2);
 263
 264       tab_float (table, 3, 2 + i * 4, TAB_RIGHT,
 265                  ws[i].positives.sum / (double) ws[i].positives.n, 8, 2);
 266
 267     }
 268
 269   tab_hline (table, TAL_2, 0, table->nc - 1, 1);
 270   tab_vline (table, TAL_2, 2, 0, table->nr - 1);
 271
 272
 273   tab_submit (table);
 274 }
 275
 276
 277 static void
 278 show_tests_box (const struct wilcoxon_state *ws,
 279                 const struct two_sample_test *t2s,
 280                 bool exact,
 281                 double timer UNUSED
 282                 )
 283 {
 284   size_t i;
 285   struct tab_table *table = tab_create (1 + t2s->n_pairs, exact ? 5 : 3, 0);
 286
 287   tab_dim (table, tab_natural_dimensions);
 288
 289   tab_title (table, _("Test Statistics"));
 290
 291   tab_headers (table, 1, 0, 1, 0);
 292
 293   /* Vertical lines inside the box */
 294   tab_box (table, 0, 0, -1, TAL_1,
 295            0, 0, table->nc - 1, tab_nr (table) - 1 );
 296
 297   /* Box around entire table */
 298   tab_box (table, TAL_2, TAL_2, -1, -1,
 299            0, 0, table->nc - 1, tab_nr (table) - 1 );
 300
 301
 302   tab_text (table,  0, 1,  TAB_LEFT, _("Z"));
 303   tab_text (table,  0, 2,  TAB_LEFT, _("Asymp. Sig (2-tailed)"));
 304
 305   if ( exact )
 306     {
 307       tab_text (table,  0, 3,  TAB_LEFT, _("Exact Sig (2-tailed)"));
 308       tab_text (table,  0, 4,  TAB_LEFT, _("Exact Sig (1-tailed)"));
 309
 310 #if 0
 311       tab_text (table,  0, 5,  TAB_LEFT, _("Point Probability"));
 312 #endif
 313     }
 314
 315   for (i = 0 ; i < t2s->n_pairs; ++i)
 316     {
 317       double z;
 318       double n = ws[i].positives.n + ws[i].negatives.n;
 319       variable_pair *vp = &t2s->pairs[i];
 320
 321       struct string pair_name;
 322       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
 323       ds_put_cstr (&pair_name, " - ");
 324       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
 325
 326
 327       tab_text (table, 1 + i, 0, TAB_CENTER, ds_cstr (&pair_name));
 328       ds_destroy (&pair_name);
 329
 330       z = MIN (ws[i].positives.sum, ws[i].negatives.sum);
 331       z -= n * (n + 1)/ 4.0;
 332
 333       z /= sqrt (n * (n + 1) * (2*n + 1)/24.0 - ws[i].tiebreaker / 48.0);
 334
 335       tab_float (table, 1 + i, 1, TAB_RIGHT, z, 8, 3);
 336
 337       tab_float (table, 1 + i, 2, TAB_RIGHT,
 338                  2.0 * gsl_cdf_ugaussian_P (z),
 339                  8, 3);
 340
 341       if (exact)
 342         {
 343           double p = LevelOfSignificanceWXMPSR (ws[i].positives.sum, n);
 344           if (p < 0)
 345             {
 346               msg (MW, ("Too many pairs to calculate exact significance."));
 347             }
 348           else
 349             {
 350               tab_float (table, 1 + i, 3, TAB_RIGHT, p, 8, 3);
 351               tab_float (table, 1 + i, 4, TAB_RIGHT, p / 2.0, 8, 3);
 352             }
 353         }
 354     }
 355
 356   tab_hline (table, TAL_2, 0, table->nc - 1, 1);
 357   tab_vline (table, TAL_2, 1, 0, table->nr - 1);
 358
 359
 360   tab_submit (table);
 361 }