Added result_class parameter to tab_double and updated all callers. Removed tab_fixed
[pspp] / src / language / stats / mcnemar.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 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 "mcnemar.h"
20
21 #include <gsl/gsl_cdf.h>
22 #include <gsl/gsl_randist.h>
23
24 #include "data/casereader.h"
25 #include "data/dataset.h"
26 #include "data/dictionary.h"
27 #include "data/format.h"
28 #include "data/missing-values.h"
29 #include "data/variable.h"
30 #include "language/stats/npar.h"
31 #include "libpspp/str.h"
32 #include "output/tab.h"
33 #include "libpspp/message.h"
34  
35 #include "gl/minmax.h"
36 #include "gl/xalloc.h"
37
38 #include "data/value.h"
39
40 #include "gettext.h"
41 #define _(msgid) gettext (msgid)
42
43
44 struct mcnemar
45 {
46   union value val0;
47   union value val1;
48
49   double n00;
50   double n01;
51   double n10;
52   double n11;
53 };
54
55 static void
56 output_freq_table (variable_pair *vp,
57                    const struct mcnemar *param,
58                    const struct dictionary *dict);
59
60
61 static void
62 output_statistics_table (const struct two_sample_test *t2s,
63                          const struct mcnemar *param,
64                          const struct dictionary *dict);
65
66
67 void
68 mcnemar_execute (const struct dataset *ds,
69                   struct casereader *input,
70                   enum mv_class exclude,
71                   const struct npar_test *test,
72                   bool exact UNUSED,
73                   double timer UNUSED)
74 {
75   int i;
76   bool warn = true;
77
78   const struct dictionary *dict = dataset_dict (ds);
79
80   const struct two_sample_test *t2s = UP_CAST (test, const struct two_sample_test, parent);
81   struct ccase *c;
82   
83   struct casereader *r = input;
84
85   struct mcnemar *mc = xcalloc (t2s->n_pairs, sizeof *mc);
86
87   for (i = 0 ; i < t2s->n_pairs; ++i )
88     {
89       mc[i].val0.f = mc[i].val1.f = SYSMIS;
90     }
91
92   for (; (c = casereader_read (r)) != NULL; case_unref (c))
93     {
94       const double weight = dict_get_case_weight (dict, c, &warn);
95
96       for (i = 0 ; i < t2s->n_pairs; ++i )
97         {
98           variable_pair *vp = &t2s->pairs[i];
99           const union value *value0 = case_data (c, (*vp)[0]);
100           const union value *value1 = case_data (c, (*vp)[1]);
101
102           if (var_is_value_missing ((*vp)[0], value0, exclude))
103             continue;
104
105           if (var_is_value_missing ((*vp)[1], value1, exclude))
106             continue;
107
108
109           if ( mc[i].val0.f == SYSMIS)
110             {
111               if (mc[i].val1.f != value0->f)
112                 mc[i].val0.f = value0->f;
113               else if (mc[i].val1.f != value1->f)
114                 mc[i].val0.f = value1->f;
115             }
116
117           if ( mc[i].val1.f == SYSMIS)
118             {
119               if (mc[i].val0.f != value1->f)
120                 mc[i].val1.f = value1->f;
121               else if (mc[i].val0.f != value0->f)
122                 mc[i].val1.f = value0->f;
123             }
124
125           if (mc[i].val0.f == value0->f && mc[i].val0.f == value1->f)
126             {
127               mc[i].n00 += weight;
128             }
129           else if ( mc[i].val0.f == value0->f && mc[i].val1.f == value1->f)
130             {
131               mc[i].n10 += weight;
132             }
133           else if ( mc[i].val1.f == value0->f && mc[i].val0.f == value1->f)
134             {
135               mc[i].n01 += weight;
136             }
137           else if ( mc[i].val1.f == value0->f && mc[i].val1.f == value1->f)
138             {
139               mc[i].n11 += weight;
140             }
141           else
142             {
143               msg (ME, _("The McNemar test is appropriate only for dichotomous variables"));
144             }
145         }
146     }
147
148   casereader_destroy (r);
149
150   for (i = 0 ; i < t2s->n_pairs ; ++i)
151     output_freq_table (&t2s->pairs[i], mc + i, dict);
152
153   output_statistics_table (t2s, mc, dict);
154
155   free (mc);
156 }
157
158
159 static void
160 output_freq_table (variable_pair *vp,
161                    const struct mcnemar *param,
162                    const struct dictionary *dict)
163 {
164   const int header_rows = 2;
165   const int header_cols = 1;
166
167   struct tab_table *table = tab_create (header_cols + 2, header_rows + 2);
168
169   const struct variable *wv = dict_get_weight (dict);
170   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : & F_8_0;
171
172
173   struct string pair_name;
174   struct string val1str ;
175   struct string val0str ;
176
177   tab_set_format (table, RC_WEIGHT, wfmt);
178
179   ds_init_empty (&val0str);
180   ds_init_empty (&val1str);
181   
182   var_append_value_name ((*vp)[0], &param->val0, &val0str);
183   var_append_value_name ((*vp)[1], &param->val1, &val1str);
184
185   ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
186   ds_put_cstr (&pair_name, " & ");
187   ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
188
189   tab_title (table, "%s", ds_cstr (&pair_name));
190
191   ds_destroy (&pair_name);
192
193   tab_headers (table, header_cols, 0, header_rows, 0);
194
195   /* Vertical lines inside the box */
196   tab_box (table, 0, 0, -1, TAL_1,
197            1, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
198
199   /* Box around entire table */
200   tab_box (table, TAL_2, TAL_2, -1, -1,
201            0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
202
203   tab_vline (table, TAL_2, header_cols, 0, tab_nr (table) - 1);
204   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, header_rows);
205
206   tab_text (table,  0, 0,  TAB_CENTER, var_to_string ((*vp)[0]));
207
208   tab_joint_text (table,  1, 0,  2, 0, TAB_CENTER, var_to_string ((*vp)[1]));
209   tab_hline (table, TAL_1, 1, tab_nc (table) - 1, 1);
210
211
212   tab_text (table, 0, header_rows + 0, TAB_LEFT, ds_cstr (&val0str));
213   tab_text (table, 0, header_rows + 1, TAB_LEFT, ds_cstr (&val1str));
214
215   tab_text (table, header_cols + 0, 1, TAB_LEFT, ds_cstr (&val0str));
216   tab_text (table, header_cols + 1, 1, TAB_LEFT, ds_cstr (&val1str));
217
218   tab_double (table, header_cols + 0, header_rows + 0, TAB_RIGHT, param->n00, NULL, RC_WEIGHT);
219   tab_double (table, header_cols + 0, header_rows + 1, TAB_RIGHT, param->n01, NULL, RC_WEIGHT);
220   tab_double (table, header_cols + 1, header_rows + 0, TAB_RIGHT, param->n10, NULL, RC_WEIGHT);
221   tab_double (table, header_cols + 1, header_rows + 1, TAB_RIGHT, param->n11, NULL, RC_WEIGHT);
222
223   tab_submit (table);
224
225   ds_destroy (&val0str);
226   ds_destroy (&val1str);
227 }
228
229
230 static void
231 output_statistics_table (const struct two_sample_test *t2s,
232                          const struct mcnemar *mc,
233                          const struct dictionary *dict)
234 {
235   int i;
236
237   struct tab_table *table = tab_create (5, t2s->n_pairs + 1);
238
239   const struct variable *wv = dict_get_weight (dict);
240   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : & F_8_0;
241
242   tab_title (table, _("Test Statistics"));
243   tab_set_format (table, RC_WEIGHT, wfmt);
244
245   tab_headers (table, 0, 1,  0, 1);
246
247   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, 1);
248   tab_vline (table, TAL_2, 1, 0, tab_nr (table) - 1);
249
250
251   /* Vertical lines inside the box */
252   tab_box (table, -1, -1, -1, TAL_1,
253            0, 0,
254            tab_nc (table) - 1, tab_nr (table) - 1);
255
256   /* Box around entire table */
257   tab_box (table, TAL_2, TAL_2, -1, -1,
258            0, 0, tab_nc (table) - 1,
259            tab_nr (table) - 1);
260
261   tab_text (table,  1, 0, TAT_TITLE | TAB_CENTER,
262             _("N"));
263
264   tab_text (table,  2, 0, TAT_TITLE | TAB_CENTER,
265             _("Exact Sig. (2-tailed)"));
266
267   tab_text (table,  3, 0, TAT_TITLE | TAB_CENTER,
268             _("Exact Sig. (1-tailed)"));
269
270   tab_text (table,  4, 0, TAT_TITLE | TAB_CENTER,
271             _("Point Probability"));
272
273   for (i = 0 ; i < t2s->n_pairs; ++i)
274     {
275       double sig;
276       variable_pair *vp = &t2s->pairs[i];
277
278       struct string pair_name;
279       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
280       ds_put_cstr (&pair_name, " & ");
281       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
282
283       tab_text (table,  0, 1 + i, TAB_LEFT, ds_cstr (&pair_name));
284       ds_destroy (&pair_name);
285
286       tab_double (table, 1, 1 + i, TAB_RIGHT, mc[i].n00 + mc[i].n01 + mc[i].n10 + mc[i].n11, NULL, RC_WEIGHT);
287
288       sig = gsl_cdf_binomial_P (mc[i].n01,  0.5,  mc[i].n01 + mc[i].n10);
289
290       tab_double (table, 2, 1 + i, TAB_RIGHT, 2 * sig, NULL, RC_PVALUE);
291       tab_double (table, 3, 1 + i, TAB_RIGHT, sig, NULL, RC_PVALUE);
292
293       tab_double (table, 4, 1 + i, TAB_RIGHT, gsl_ran_binomial_pdf (mc[i].n01, 0.5, mc[i].n01 + mc[i].n10),
294                   NULL, RC_OTHER);
295     }
296
297   tab_submit (table);
298 }