pintos-os.org Git - pspp/blob - src/language/stats/mann-whitney.c

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