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