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 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 }