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