GLM: Disable the interactions option until it actually gives the right results
[pspp] / 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/dataset.h"
29 #include "data/dictionary.h"
30 #include "data/format.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 }