src/language/stats/wilcoxon.c

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