1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2010, 2011 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 "language/stats/mann-whitney.h"
21 #include <gsl/gsl_cdf.h>
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"
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);
59 belongs_to_test (const struct ccase *c, void *aux)
61 const struct n_sample_test *nst = aux;
63 const union value *group = case_data (c, nst->indep_var);
64 const size_t group_var_width = var_get_width (nst->indep_var);
66 if ( value_equal (group, &nst->val1, group_var_width))
69 if ( value_equal (group, &nst->val2, group_var_width))
78 mann_whitney_execute (const struct dataset *ds,
79 struct casereader *input,
80 enum mv_class exclude,
81 const struct npar_test *test,
86 const struct dictionary *dict = dataset_dict (ds);
87 const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
89 const struct caseproto *proto = casereader_get_proto (input);
90 size_t rank_idx = caseproto_get_n_widths (proto);
92 struct mw *mw = xcalloc (nst->n_vars, sizeof *mw);
94 for (i = 0; i < nst->n_vars; ++i)
96 double tiebreaker = 0.0;
98 enum rank_error rerr = 0;
99 struct casereader *rr;
101 const struct variable *var = nst->vars[i];
103 struct casereader *reader =
104 casereader_create_filter_func (casereader_clone (input),
107 CONST_CAST (struct n_sample_test *, nst),
110 reader = casereader_create_filter_missing (reader, &var, 1,
114 reader = sort_execute_1var (reader, var);
116 rr = casereader_create_append_rank (reader, var,
117 dict_get_weight (dict),
119 distinct_callback, &tiebreaker);
121 for (; (c = casereader_read (rr)); case_unref (c))
123 const union value *group = case_data (c, nst->indep_var);
124 const size_t group_var_width = var_get_width (nst->indep_var);
125 const double rank = case_data_idx (c, rank_idx)->f;
127 if ( value_equal (group, &nst->val1, group_var_width))
129 mw[i].rank_sum[0] += rank;
130 mw[i].n[0] += dict_get_case_weight (dict, c, &warn);
132 else if ( value_equal (group, &nst->val2, group_var_width))
134 mw[i].rank_sum[1] += rank;
135 mw[i].n[1] += dict_get_case_weight (dict, c, &warn);
138 casereader_destroy (rr);
143 struct mw *mwv = &mw[i];
145 mwv->u = mwv->n[0] * mwv->n[1] ;
146 mwv->u += mwv->n[0] * (mwv->n[0] + 1) / 2.0;
147 mwv->u -= mwv->rank_sum[0];
149 mwv->w = mwv->rank_sum[1];
150 if ( mwv->u > mwv->n[0] * mwv->n[1] / 2.0)
152 mwv->u = mwv->n[0] * mwv->n[1] - mwv->u;
153 mwv->w = mwv->rank_sum[0];
155 mwv->z = mwv->u - mwv->n[0] * mwv->n[1] / 2.0;
156 n = mwv->n[0] + mwv->n[1];
157 denominator = pow3(n) - n;
159 denominator -= tiebreaker;
160 denominator *= mwv->n[0] * mwv->n[1];
161 denominator /= n * (n - 1);
163 mwv->z /= sqrt (denominator);
166 casereader_destroy (input);
168 show_ranks_box (nst, mw);
169 show_statistics_box (nst, mw, exact);
177 #define _(msgid) gettext (msgid)
180 show_ranks_box (const struct n_sample_test *nst, const struct mw *mwv)
183 const int row_headers = 1;
184 const int column_headers = 2;
185 struct tab_table *table =
186 tab_create (row_headers + 7, column_headers + nst->n_vars);
188 struct string g1str, g2str;;
189 ds_init_empty (&g1str);
190 var_append_value_name (nst->indep_var, &nst->val1, &g1str);
192 ds_init_empty (&g2str);
193 var_append_value_name (nst->indep_var, &nst->val2, &g2str);
195 tab_headers (table, row_headers, 0, column_headers, 0);
197 tab_title (table, _("Ranks"));
199 /* Vertical lines inside the box */
200 tab_box (table, 1, 0, -1, TAL_1,
201 row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
203 /* Box around the table */
204 tab_box (table, TAL_2, TAL_2, -1, -1,
205 0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
207 tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
208 tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
210 tab_hline (table, TAL_1, row_headers, tab_nc (table) -1, 1);
212 tab_text (table, 1, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
213 tab_text (table, 2, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
214 tab_text (table, 3, 1, TAT_TITLE | TAB_CENTER, _("Total"));
215 tab_joint_text (table, 1, 0, 3, 0,
216 TAT_TITLE | TAB_CENTER, _("N"));
217 tab_vline (table, TAL_2, 4, 0, tab_nr (table) - 1);
219 tab_text (table, 4, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
220 tab_text (table, 5, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
221 tab_joint_text (table, 4, 0, 5, 0,
222 TAT_TITLE | TAB_CENTER, _("Mean Rank"));
223 tab_vline (table, TAL_2, 6, 0, tab_nr (table) - 1);
225 tab_text (table, 6, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
226 tab_text (table, 7, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
227 tab_joint_text (table, 6, 0, 7, 0,
228 TAT_TITLE | TAB_CENTER, _("Sum of Ranks"));
233 for (i = 0 ; i < nst->n_vars ; ++i)
235 const struct mw *mw = &mwv[i];
236 tab_text (table, 0, column_headers + i, TAT_TITLE,
237 var_to_string (nst->vars[i]));
239 tab_double (table, 1, column_headers + i, 0,
240 mw->n[0], NULL, RC_OTHER);
242 tab_double (table, 2, column_headers + i, 0,
243 mw->n[1], NULL, RC_OTHER);
245 tab_double (table, 3, column_headers + i, 0,
246 mw->n[1] + mw->n[0], NULL, RC_OTHER);
249 tab_double (table, 4, column_headers + i, 0,
250 mw->rank_sum[0] / mw->n[0], NULL, RC_OTHER);
252 tab_double (table, 5, column_headers + i, 0,
253 mw->rank_sum[1] / mw->n[1], NULL, RC_OTHER);
256 tab_double (table, 6, column_headers + i, 0,
257 mw->rank_sum[0], NULL, RC_OTHER);
259 tab_double (table, 7, column_headers + i, 0,
260 mw->rank_sum[1], NULL, RC_OTHER);
267 show_statistics_box (const struct n_sample_test *nst, const struct mw *mwv, bool exact)
270 const int row_headers = 1;
271 const int column_headers = 1;
272 struct tab_table *table =
273 tab_create (row_headers + (exact ? 6 : 4), column_headers + nst->n_vars);
275 tab_headers (table, row_headers, 0, column_headers, 0);
277 tab_title (table, _("Test Statistics"));
279 /* Vertical lines inside the box */
280 tab_box (table, 1, 0, -1, TAL_1,
281 row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
283 /* Box around the table */
284 tab_box (table, TAL_2, TAL_2, -1, -1,
285 0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
287 tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
288 tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
290 tab_text (table, 1, 0, TAT_TITLE | TAB_CENTER, _("Mann-Whitney U"));
291 tab_text (table, 2, 0, TAT_TITLE | TAB_CENTER, _("Wilcoxon W"));
292 tab_text (table, 3, 0, TAT_TITLE | TAB_CENTER, _("Z"));
293 tab_text (table, 4, 0, TAT_TITLE | TAB_CENTER, _("Asymp. Sig. (2-tailed)"));
297 tab_text (table, 5, 0, TAT_TITLE | TAB_CENTER, _("Exact Sig. (2-tailed)"));
298 tab_text (table, 6, 0, TAT_TITLE | TAB_CENTER, _("Point Probability"));
301 for (i = 0 ; i < nst->n_vars ; ++i)
303 const struct mw *mw = &mwv[i];
305 tab_text (table, 0, column_headers + i, TAT_TITLE,
306 var_to_string (nst->vars[i]));
308 tab_double (table, 1, column_headers + i, 0,
309 mw->u, NULL, RC_OTHER);
311 tab_double (table, 2, column_headers + i, 0,
312 mw->w, NULL, RC_OTHER);
314 tab_double (table, 3, column_headers + i, 0,
315 mw->z, NULL, RC_OTHER);
317 tab_double (table, 4, column_headers + i, 0,
318 2.0 * gsl_cdf_ugaussian_P (mw->z), NULL, RC_PVALUE);