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               continue;
 128             }
 129
 130           case_data_rw (output, ws[i].absdiff)->f = fabs (d);
 131
 132           if (weight)
 133            case_data_rw (output, weightx)->f = case_data (c, weight)->f;
 134
 135           casewriter_write (writer, output);
 136         }
 137       casereader_destroy (r);
 138       ws[i].reader = casewriter_make_reader (writer);
 139     }
 140
 141   for (i = 0 ; i < t2s->n_pairs; ++i )
 142     {
 143       struct casereader *rr ;
 144       struct ccase *c;
 145       enum rank_error err = 0;
 146
 147       rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff,
 148                                           weight ? weightx : NULL, &err,
 149                                           distinct_callback, &ws[i]
 150                                           );
 151
 152       for (; (c = casereader_read (rr)) != NULL; case_unref (c))
 153         {
 154           double sign = case_data (c, ws[i].sign)->f;
 155           double rank = case_data_idx (c, weight ? 3 : 2)->f;
 156           double w = 1.0;
 157           if (weight)
 158             w = case_data (c, weightx)->f;
 159
 160           if ( sign > 0 )
 161             {
 162               ws[i].positives.sum += rank * w;
 163               ws[i].positives.n += w;
 164             }
 165           else if (sign < 0)
 166             {
 167               ws[i].negatives.sum += rank * w;
 168               ws[i].negatives.n += w;
 169             }
 170           else
 171             NOT_REACHED ();
 172         }
 173
 174       casereader_destroy (rr);
 175     }
 176
 177   casereader_destroy (input);
 178
 179   var_destroy (weightx);
 180
 181   show_ranks_box (ws, t2s);
 182   show_tests_box (ws, t2s, exact, timer);
 183
 184   for (i = 0 ; i < t2s->n_pairs; ++i )
 185     {
 186       var_destroy (ws[i].sign);
 187       var_destroy (ws[i].absdiff);
 188     }
 189
 190   free (ws);
 191 }
 192
 193
 194 \f
 195
 196 #include "gettext.h"
 197 #define _(msgid) gettext (msgid)
 198
 199 static void
 200 show_ranks_box (const struct wilcoxon_state *ws, const struct two_sample_test *t2s)
 201 {
 202   size_t i;
 203   struct tab_table *table = tab_create (5, 1 + 4 * t2s->n_pairs, 0);
 204
 205   tab_dim (table, tab_natural_dimensions);
 206
 207   tab_title (table, _("Ranks"));
 208
 209   tab_headers (table, 2, 0, 1, 0);
 210
 211   /* Vertical lines inside the box */
 212   tab_box (table, 0, 0, -1, TAL_1,
 213            1, 0, table->nc - 1, tab_nr (table) - 1 );
 214
 215   /* Box around entire table */
 216   tab_box (table, TAL_2, TAL_2, -1, -1,
 217            0, 0, table->nc - 1, tab_nr (table) - 1 );
 218
 219
 220   tab_text (table,  2, 0,  TAB_CENTER, _("N"));
 221   tab_text (table,  3, 0,  TAB_CENTER, _("Mean Rank"));
 222   tab_text (table,  4, 0,  TAB_CENTER, _("Sum of Ranks"));
 223
 224
 225   for (i = 0 ; i < t2s->n_pairs; ++i)
 226     {
 227       variable_pair *vp = &t2s->pairs[i];
 228
 229       struct string pair_name;
 230       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
 231       ds_put_cstr (&pair_name, " - ");
 232       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
 233
 234       tab_text (table, 1, 1 + i * 4, TAB_LEFT, _("Negative Ranks"));
 235       tab_text (table, 1, 2 + i * 4, TAB_LEFT, _("Positive Ranks"));
 236       tab_text (table, 1, 3 + i * 4, TAB_LEFT, _("Ties"));
 237       tab_text (table, 1, 4 + i * 4, TAB_LEFT, _("Total"));
 238
 239       tab_hline (table, TAL_1, 0, table->nc - 1, 1 + i * 4);
 240
 241
 242       tab_text (table, 0, 1 + i * 4, TAB_LEFT, ds_cstr (&pair_name));
 243       ds_destroy (&pair_name);
 244
 245
 246       /* N */
 247       tab_float (table, 2, 1 + i * 4, TAB_RIGHT, ws[i].negatives.n, 8, 0);
 248       tab_float (table, 2, 2 + i * 4, TAB_RIGHT, ws[i].positives.n, 8, 0);
 249       tab_float (table, 2, 3 + i * 4, TAB_RIGHT, ws[i].n_zeros, 8, 0);
 250
 251       tab_float (table, 2, 4 + i * 4, TAB_RIGHT,
 252                  ws[i].n_zeros + ws[i].positives.n + ws[i].negatives.n, 8, 0);
 253
 254       /* Sums */
 255       tab_float (table, 4, 1 + i * 4, TAB_RIGHT, ws[i].negatives.sum, 8, 2);
 256       tab_float (table, 4, 2 + i * 4, TAB_RIGHT, ws[i].positives.sum, 8, 2);
 257
 258
 259       /* Means */
 260       tab_float (table, 3, 1 + i * 4, TAB_RIGHT,
 261                  ws[i].negatives.sum / (double) ws[i].negatives.n, 8, 2);
 262
 263       tab_float (table, 3, 2 + i * 4, TAB_RIGHT,
 264                  ws[i].positives.sum / (double) ws[i].positives.n, 8, 2);
 265
 266     }
 267
 268   tab_hline (table, TAL_2, 0, table->nc - 1, 1);
 269   tab_vline (table, TAL_2, 2, 0, table->nr - 1);
 270
 271
 272   tab_submit (table);
 273 }
 274
 275
 276 static void
 277 show_tests_box (const struct wilcoxon_state *ws,
 278                 const struct two_sample_test *t2s,
 279                 bool exact,
 280                 double timer UNUSED
 281                 )
 282 {
 283   size_t i;
 284   struct tab_table *table = tab_create (1 + t2s->n_pairs, exact ? 5 : 3, 0);
 285
 286   tab_dim (table, tab_natural_dimensions);
 287
 288   tab_title (table, _("Test Statistics"));
 289
 290   tab_headers (table, 1, 0, 1, 0);
 291
 292   /* Vertical lines inside the box */
 293   tab_box (table, 0, 0, -1, TAL_1,
 294            0, 0, table->nc - 1, tab_nr (table) - 1 );
 295
 296   /* Box around entire table */
 297   tab_box (table, TAL_2, TAL_2, -1, -1,
 298            0, 0, table->nc - 1, tab_nr (table) - 1 );
 299
 300
 301   tab_text (table,  0, 1,  TAB_LEFT, _("Z"));
 302   tab_text (table,  0, 2,  TAB_LEFT, _("Asymp. Sig (2-tailed)"));
 303
 304   if ( exact )
 305     {
 306       tab_text (table,  0, 3,  TAB_LEFT, _("Exact Sig (2-tailed)"));
 307       tab_text (table,  0, 4,  TAB_LEFT, _("Exact Sig (1-tailed)"));
 308
 309 #if 0
 310       tab_text (table,  0, 5,  TAB_LEFT, _("Point Probability"));
 311 #endif
 312     }
 313
 314   for (i = 0 ; i < t2s->n_pairs; ++i)
 315     {
 316       double z;
 317       double n = ws[i].positives.n + ws[i].negatives.n;
 318       variable_pair *vp = &t2s->pairs[i];
 319
 320       struct string pair_name;
 321       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
 322       ds_put_cstr (&pair_name, " - ");
 323       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
 324
 325
 326       tab_text (table, 1 + i, 0, TAB_CENTER, ds_cstr (&pair_name));
 327       ds_destroy (&pair_name);
 328
 329       z = MIN (ws[i].positives.sum, ws[i].negatives.sum);
 330       z -= n * (n + 1)/ 4.0;
 331
 332       z /= sqrt (n * (n + 1) * (2*n + 1)/24.0 - ws[i].tiebreaker / 48.0);
 333
 334       tab_float (table, 1 + i, 1, TAB_RIGHT, z, 8, 3);
 335
 336       tab_float (table, 1 + i, 2, TAB_RIGHT,
 337                  2.0 * gsl_cdf_ugaussian_P (z),
 338                  8, 3);
 339
 340       if (exact)
 341         {
 342           double p = LevelOfSignificanceWXMPSR (ws[i].positives.sum, n);
 343           if (p < 0)
 344             {
 345               msg (MW, ("Too many pairs to calculate exact significance."));
 346             }
 347           else
 348             {
 349               tab_float (table, 1 + i, 3, TAB_RIGHT, p, 8, 3);
 350               tab_float (table, 1 + i, 4, TAB_RIGHT, p / 2.0, 8, 3);
 351             }
 352         }
 353     }
 354
 355   tab_hline (table, TAL_2, 0, table->nc - 1, 1);
 356   tab_vline (table, TAL_2, 1, 0, table->nr - 1);
 357
 358
 359   tab_submit (table);
 360 }