1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2010 Free Software Foundation, Inc.
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.
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.
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/>. */
19 #include "mann-whitney.h"
21 #include <gsl/gsl_cdf.h>
23 #include <libpspp/cast.h>
24 #include <libpspp/misc.h>
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>
34 /* Calculates the adjustment necessary for tie compensation */
36 distinct_callback (double v UNUSED, casenumber t, double w UNUSED, void *aux)
38 double *tiebreaker = aux;
40 *tiebreaker += (pow3 (t) - t) / 12.0;
48 double u; /* The Mann-Whitney U statistic */
49 double w; /* The Wilcoxon Rank Sum W statistic */
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);
58 mann_whitney_execute (const struct dataset *ds,
59 struct casereader *input,
60 enum mv_class exclude,
61 const struct npar_test *test,
66 const struct dictionary *dict = dataset_dict (ds);
67 const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
69 const struct caseproto *proto = casereader_get_proto (input);
70 size_t rank_idx = caseproto_get_n_widths (proto);
72 struct mw *mw = xcalloc (nst->n_vars, sizeof *mw);
74 for (i = 0; i < nst->n_vars; ++i)
76 double tiebreaker = 0.0;
78 enum rank_error rerr = 0;
79 struct casereader *rr;
81 const struct variable *var = nst->vars[i];
83 struct casereader *reader =
84 sort_execute_1var (casereader_clone (input), var);
86 rr = casereader_create_append_rank (reader, var,
87 dict_get_weight (dict),
89 distinct_callback, &tiebreaker);
91 for (; (c = casereader_read (rr)); case_unref (c))
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;
98 if ( var_is_value_missing (var, val, exclude))
101 if ( value_equal (group, &nst->val1, group_var_width))
103 mw[i].rank_sum[0] += rank;
104 mw[i].n[0] += dict_get_case_weight (dict, c, &warn);
106 else if ( value_equal (group, &nst->val2, group_var_width))
108 mw[i].rank_sum[1] += rank;
109 mw[i].n[1] += dict_get_case_weight (dict, c, &warn);
112 casereader_destroy (rr);
117 struct mw *mwv = &mw[i];
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];
123 mwv->w = mwv->rank_sum[1];
124 if ( mwv->u > mwv->n[0] * mwv->n[1] / 2.0)
126 mwv->u = mwv->n[0] * mwv->n[1] - mwv->u;
127 mwv->w = mwv->rank_sum[0];
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;
133 denominator -= tiebreaker;
134 denominator *= mwv->n[0] * mwv->n[1];
135 denominator /= n * (n - 1);
137 mwv->z /= sqrt (denominator);
140 casereader_destroy (input);
142 show_ranks_box (nst, mw);
143 show_statistics_box (nst, mw, exact);
150 #include <output/tab.h>
152 #define _(msgid) gettext (msgid)
155 show_ranks_box (const struct n_sample_test *nst, const struct mw *mwv)
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);
163 struct string g1str, g2str;;
164 ds_init_empty (&g1str);
165 var_append_value_name (nst->indep_var, &nst->val1, &g1str);
167 ds_init_empty (&g2str);
168 var_append_value_name (nst->indep_var, &nst->val2, &g2str);
170 tab_headers (table, row_headers, 0, column_headers, 0);
172 tab_title (table, _("Ranks"));
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 );
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 );
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);
185 tab_hline (table, TAL_1, row_headers, tab_nc (table) -1, 1);
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);
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);
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"));
208 for (i = 0 ; i < nst->n_vars ; ++i)
210 const struct mw *mw = &mwv[i];
211 tab_text (table, 0, column_headers + i, TAT_TITLE,
212 var_to_string (nst->vars[i]));
214 tab_double (table, 1, column_headers + i, 0,
217 tab_double (table, 2, column_headers + i, 0,
220 tab_double (table, 3, column_headers + i, 0,
221 mw->n[1] + mw->n[0], 0);
224 tab_double (table, 4, column_headers + i, 0,
225 mw->rank_sum[0] / mw->n[0], 0);
227 tab_double (table, 5, column_headers + i, 0,
228 mw->rank_sum[1] / mw->n[1], 0);
231 tab_double (table, 6, column_headers + i, 0,
234 tab_double (table, 7, column_headers + i, 0,
242 show_statistics_box (const struct n_sample_test *nst, const struct mw *mwv, bool exact)
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);
250 tab_headers (table, row_headers, 0, column_headers, 0);
252 tab_title (table, _("Test Statistics"));
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 );
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 );
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);
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)"));
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"));
276 for (i = 0 ; i < nst->n_vars ; ++i)
278 const struct mw *mw = &mwv[i];
280 tab_text (table, 0, column_headers + i, TAT_TITLE,
281 var_to_string (nst->vars[i]));
283 tab_double (table, 1, column_headers + i, 0,
286 tab_double (table, 2, column_headers + i, 0,
289 tab_double (table, 3, column_headers + i, 0,
292 tab_double (table, 4, column_headers + i, 0,
293 2.0 * gsl_cdf_ugaussian_P (mw->z), 0);