src/language/stats/mann-whitney.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2010 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 #include <config.h>
  18
  19 #include "mann-whitney.h"
  20
  21 #include <gsl/gsl_cdf.h>
  22
  23 #include <libpspp/cast.h>
  24 #include <libpspp/misc.h>
  25
  26 #include <math/sort.h>
  27 #include <data/case.h>
  28 #include <data/dictionary.h>
  29 #include <data/procedure.h>
  30 #include <data/casereader.h>
  31 #include <data/format.h>
  32 #include <data/variable.h>
  33
  34 /* Calculates the adjustment necessary for tie compensation */
  35 static void
  36 distinct_callback (double v UNUSED, casenumber t, double w UNUSED, void *aux)
  37 {
  38   double *tiebreaker = aux;
  39
  40   *tiebreaker += (pow3 (t) - t) / 12.0;
  41 }
  42
  43 struct mw
  44 {
  45   double rank_sum[2];
  46   double n[2];
  47
  48   double u;  /* The Mann-Whitney U statistic */
  49   double w;  /* The Wilcoxon Rank Sum W statistic */
  50   double z;
  51 };
  52
  53 static void show_ranks_box (const struct n_sample_test *nst, const struct mw *mw);
  54 static void show_statistics_box (const struct n_sample_test *nst, const struct mw *mw, bool exact);
  55
  56
  57 void
  58 mann_whitney_execute (const struct dataset *ds,
  59                       struct casereader *input,
  60                       enum mv_class exclude,
  61                       const struct npar_test *test,
  62                       bool exact,
  63                       double timer UNUSED)
  64 {
  65   int i;
  66   const struct dictionary *dict = dataset_dict (ds);
  67   const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
  68
  69   const struct caseproto *proto = casereader_get_proto (input);
  70   size_t rank_idx = caseproto_get_n_widths (proto);
  71
  72   struct mw *mw = xcalloc (nst->n_vars, sizeof *mw);
  73
  74   for (i = 0; i < nst->n_vars; ++i)
  75     {
  76       double tiebreaker = 0.0;
  77       bool warn = true;
  78       enum rank_error rerr = 0;
  79       struct casereader *rr;
  80       struct ccase *c;
  81       const struct variable *var = nst->vars[i];
  82
  83       struct casereader *reader =
  84         sort_execute_1var (casereader_clone (input), var);
  85
  86       rr = casereader_create_append_rank (reader, var,
  87                                           dict_get_weight (dict),
  88                                           &rerr,
  89                                           distinct_callback, &tiebreaker);
  90
  91       for (; (c = casereader_read (rr)); case_unref (c))
  92         {
  93           const union value *val = case_data (c, var);
  94           const union value *group = case_data (c, nst->indep_var);
  95           const size_t group_var_width = var_get_width (nst->indep_var);
  96           const double rank = case_data_idx (c, rank_idx)->f;
  97
  98           if ( var_is_value_missing (var, val, exclude))
  99             continue;
 100
 101           if ( value_equal (group, &nst->val1, group_var_width))
 102             {
 103               mw[i].rank_sum[0] += rank;
 104               mw[i].n[0] += dict_get_case_weight (dict, c, &warn);
 105             }
 106           else if ( value_equal (group, &nst->val2, group_var_width))
 107             {
 108               mw[i].rank_sum[1] += rank;
 109               mw[i].n[1] += dict_get_case_weight (dict, c, &warn);
 110             }
 111         }
 112       casereader_destroy (rr);
 113
 114       {
 115         double n;
 116         double denominator;
 117         struct mw *mwv = &mw[i];
 118
 119         mwv->u = mwv->n[0] * mwv->n[1] ;
 120         mwv->u += mwv->n[0] * (mwv->n[0] + 1) / 2.0;
 121         mwv->u -= mwv->rank_sum[0];
 122
 123         mwv->w = mwv->rank_sum[1];
 124         if ( mwv->u > mwv->n[0] * mwv->n[1] / 2.0)
 125           {
 126             mwv->u =  mwv->n[0] * mwv->n[1] - mwv->u;
 127             mwv->w = mwv->rank_sum[0];
 128           }
 129         mwv->z = mwv->u - mwv->n[0] * mwv->n[1] / 2.0;
 130         n = mwv->n[0] + mwv->n[1];
 131         denominator = pow3(n) - n;
 132         denominator /= 12;
 133         denominator -= tiebreaker;
 134         denominator *= mwv->n[0] * mwv->n[1];
 135         denominator /= n * (n - 1);
 136
 137         mwv->z /= sqrt (denominator);
 138       }
 139     }
 140   casereader_destroy (input);
 141
 142   show_ranks_box (nst, mw);
 143   show_statistics_box (nst, mw, exact);
 144
 145   free (mw);
 146 }
 147
 148 \f
 149
 150 #include <output/tab.h>
 151 #include "gettext.h"
 152 #define _(msgid) gettext (msgid)
 153
 154 static void
 155 show_ranks_box (const struct n_sample_test *nst, const struct mw *mwv)
 156 {
 157   int i;
 158   const int row_headers = 1;
 159   const int column_headers = 2;
 160   struct tab_table *table =
 161     tab_create (row_headers + 7, column_headers + nst->n_vars);
 162
 163   struct string g1str, g2str;;
 164   ds_init_empty (&g1str);
 165   var_append_value_name (nst->indep_var, &nst->val1, &g1str);
 166
 167   ds_init_empty (&g2str);
 168   var_append_value_name (nst->indep_var, &nst->val2, &g2str);
 169
 170   tab_headers (table, row_headers, 0, column_headers, 0);
 171
 172   tab_title (table, _("Ranks"));
 173
 174   /* Vertical lines inside the box */
 175   tab_box (table, 1, 0, -1, TAL_1,
 176            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
 177
 178   /* Box around the table */
 179   tab_box (table, TAL_2, TAL_2, -1, -1,
 180            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 181
 182   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 183   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 184
 185   tab_hline (table, TAL_1, row_headers, tab_nc (table) -1, 1);
 186
 187   tab_text (table, 1, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
 188   tab_text (table, 2, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
 189   tab_text (table, 3, 1, TAT_TITLE | TAB_CENTER, _("Total"));
 190   tab_joint_text (table, 1, 0, 3, 0,
 191                   TAT_TITLE | TAB_CENTER, _("N"));
 192   tab_vline (table, TAL_2, 4, 0, tab_nr (table) - 1);
 193
 194   tab_text (table, 4, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
 195   tab_text (table, 5, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
 196   tab_joint_text (table, 4, 0, 5, 0,
 197                   TAT_TITLE | TAB_CENTER, _("Mean Rank"));
 198   tab_vline (table, TAL_2, 6, 0, tab_nr (table) - 1);
 199
 200   tab_text (table, 6, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
 201   tab_text (table, 7, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
 202   tab_joint_text (table, 6, 0, 7, 0,
 203                   TAT_TITLE | TAB_CENTER, _("Sum of Ranks"));
 204
 205   ds_destroy (&g1str);
 206   ds_destroy (&g2str);
 207
 208   for (i = 0 ; i < nst->n_vars ; ++i)
 209     {
 210       const struct mw *mw = &mwv[i];
 211       tab_text (table, 0, column_headers + i, TAT_TITLE,
 212                 var_to_string (nst->vars[i]));
 213
 214       tab_double (table, 1, column_headers + i, 0,
 215                   mw->n[0], 0);
 216
 217       tab_double (table, 2, column_headers + i, 0,
 218                   mw->n[1], 0);
 219
 220       tab_double (table, 3, column_headers + i, 0,
 221                   mw->n[1] + mw->n[0], 0);
 222
 223       /* Mean Ranks */
 224       tab_double (table, 4, column_headers + i, 0,
 225                   mw->rank_sum[0] / mw->n[0], 0);
 226
 227       tab_double (table, 5, column_headers + i, 0,
 228                   mw->rank_sum[1] / mw->n[1], 0);
 229
 230       /* Sum of Ranks */
 231       tab_double (table, 6, column_headers + i, 0,
 232                   mw->rank_sum[0], 0);
 233
 234       tab_double (table, 7, column_headers + i, 0,
 235                   mw->rank_sum[1], 0);
 236     }
 237
 238   tab_submit (table);
 239 }
 240
 241 static void
 242 show_statistics_box (const struct n_sample_test *nst, const struct mw *mwv, bool exact)
 243 {
 244   int i;
 245   const int row_headers = 1;
 246   const int column_headers = 1;
 247   struct tab_table *table =
 248     tab_create (row_headers + (exact ? 6 : 4), column_headers + nst->n_vars);
 249
 250   tab_headers (table, row_headers, 0, column_headers, 0);
 251
 252   tab_title (table, _("Test Statistics"));
 253
 254   /* Vertical lines inside the box */
 255   tab_box (table, 1, 0, -1, TAL_1,
 256            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
 257
 258   /* Box around the table */
 259   tab_box (table, TAL_2, TAL_2, -1, -1,
 260            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 261
 262   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 263   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 264
 265   tab_text (table, 1, 0, TAT_TITLE | TAB_CENTER, _("Mann-Whitney U"));
 266   tab_text (table, 2, 0, TAT_TITLE | TAB_CENTER, _("Wilcoxon W"));
 267   tab_text (table, 3, 0, TAT_TITLE | TAB_CENTER, _("Z"));
 268   tab_text (table, 4, 0, TAT_TITLE | TAB_CENTER, _("Asymp. Sig. (2-tailed)"));
 269
 270   if (exact)
 271     {
 272       tab_text (table, 5, 0, TAT_TITLE | TAB_CENTER, _("Exact Sig. (2-tailed)"));
 273       tab_text (table, 6, 0, TAT_TITLE | TAB_CENTER, _("Point Probability"));
 274     }
 275
 276   for (i = 0 ; i < nst->n_vars ; ++i)
 277     {
 278       const struct mw *mw = &mwv[i];
 279
 280       tab_text (table, 0, column_headers + i, TAT_TITLE,
 281                 var_to_string (nst->vars[i]));
 282
 283       tab_double (table, 1, column_headers + i, 0,
 284                   mw->u, 0);
 285
 286       tab_double (table, 2, column_headers + i, 0,
 287                   mw->w, 0);
 288
 289       tab_double (table, 3, column_headers + i, 0,
 290                   mw->z, 0);
 291
 292       tab_double (table, 4, column_headers + i, 0,
 293                   2.0 * gsl_cdf_ugaussian_P (mw->z), 0);
 294     }
 295
 296   tab_submit (table);
 297 }