Constness
[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 (sizeof *mc, t2s->n_pairs);
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   ds_init_empty (&val0str);
178   ds_init_empty (&val1str);
179   
180   var_append_value_name ((*vp)[0], &param->val0, &val0str);
181   var_append_value_name ((*vp)[1], &param->val1, &val1str);
182
183   ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
184   ds_put_cstr (&pair_name, " & ");
185   ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
186
187   tab_title (table, ds_cstr (&pair_name));
188
189   ds_destroy (&pair_name);
190
191   tab_headers (table, header_cols, 0, header_rows, 0);
192
193   /* Vertical lines inside the box */
194   tab_box (table, 0, 0, -1, TAL_1,
195            1, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
196
197   /* Box around entire table */
198   tab_box (table, TAL_2, TAL_2, -1, -1,
199            0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
200
201   tab_vline (table, TAL_2, header_cols, 0, tab_nr (table) - 1);
202   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, header_rows);
203
204   tab_text (table,  0, 0,  TAB_CENTER, var_to_string ((*vp)[0]));
205
206   tab_joint_text (table,  1, 0,  2, 0, TAB_CENTER, var_to_string ((*vp)[1]));
207   tab_hline (table, TAL_1, 1, tab_nc (table) - 1, 1);
208
209
210   tab_text (table, 0, header_rows + 0, TAB_LEFT, ds_cstr (&val0str));
211   tab_text (table, 0, header_rows + 1, TAB_LEFT, ds_cstr (&val1str));
212
213   tab_text (table, header_cols + 0, 1, TAB_LEFT, ds_cstr (&val0str));
214   tab_text (table, header_cols + 1, 1, TAB_LEFT, ds_cstr (&val1str));
215
216   tab_double (table, header_cols + 0, header_rows + 0, TAB_RIGHT, param->n00, wfmt);
217   tab_double (table, header_cols + 0, header_rows + 1, TAB_RIGHT, param->n01, wfmt);
218   tab_double (table, header_cols + 1, header_rows + 0, TAB_RIGHT, param->n10, wfmt);
219   tab_double (table, header_cols + 1, header_rows + 1, TAB_RIGHT, param->n11, wfmt);
220
221   tab_submit (table);
222
223   ds_destroy (&val0str);
224   ds_destroy (&val1str);
225 }
226
227
228 static void
229 output_statistics_table (const struct two_sample_test *t2s,
230                          const struct mcnemar *mc,
231                          const struct dictionary *dict)
232 {
233   int i;
234
235   struct tab_table *table = tab_create (5, t2s->n_pairs + 1);
236
237   const struct variable *wv = dict_get_weight (dict);
238   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : & F_8_0;
239
240   tab_title (table, _("Test Statistics"));
241
242   tab_headers (table, 0, 1,  0, 1);
243
244   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, 1);
245   tab_vline (table, TAL_2, 1, 0, tab_nr (table) - 1);
246
247
248   /* Vertical lines inside the box */
249   tab_box (table, -1, -1, -1, TAL_1,
250            0, 0,
251            tab_nc (table) - 1, tab_nr (table) - 1);
252
253   /* Box around entire table */
254   tab_box (table, TAL_2, TAL_2, -1, -1,
255            0, 0, tab_nc (table) - 1,
256            tab_nr (table) - 1);
257
258   tab_text (table,  1, 0, TAT_TITLE | TAB_CENTER,
259             _("N"));
260
261   tab_text (table,  2, 0, TAT_TITLE | TAB_CENTER,
262             _("Exact Sig. (2-tailed)"));
263
264   tab_text (table,  3, 0, TAT_TITLE | TAB_CENTER,
265             _("Exact Sig. (1-tailed)"));
266
267   tab_text (table,  4, 0, TAT_TITLE | TAB_CENTER,
268             _("Point Probability"));
269
270   for (i = 0 ; i < t2s->n_pairs; ++i)
271     {
272       double sig;
273       variable_pair *vp = &t2s->pairs[i];
274
275       struct string pair_name;
276       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
277       ds_put_cstr (&pair_name, " & ");
278       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
279
280       tab_text (table,  0, 1 + i, TAB_LEFT, ds_cstr (&pair_name));
281       ds_destroy (&pair_name);
282
283       tab_double (table, 1, 1 + i, TAB_RIGHT, mc[i].n00 + mc[i].n01 + mc[i].n10 + mc[i].n11, wfmt);
284
285       sig = gsl_cdf_binomial_P (mc[i].n01,  0.5,  mc[i].n01 + mc[i].n10);
286
287       tab_double (table, 2, 1 + i, TAB_RIGHT, 2 * sig, NULL);
288       tab_double (table, 3, 1 + i, TAB_RIGHT, sig, NULL);
289
290       tab_double (table, 4, 1 + i, TAB_RIGHT, gsl_ran_binomial_pdf (mc[i].n01, 0.5, mc[i].n01 + mc[i].n10),
291                   NULL);
292     }
293
294   tab_submit (table);
295 }