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/pivot-table.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 N_(msgid) msgid
  42 #define _(msgid) gettext (msgid)
  43
  44
  45 struct mcnemar
  46 {
  47   union value val0;
  48   union value val1;
  49
  50   double n00;
  51   double n01;
  52   double n10;
  53   double n11;
  54 };
  55
  56 static void
  57 output_freq_table (variable_pair *vp,
  58                    const struct mcnemar *param,
  59                    const struct dictionary *dict);
  60
  61
  62 static void
  63 output_statistics_table (const struct two_sample_test *t2s,
  64                          const struct mcnemar *param,
  65                          const struct dictionary *dict);
  66
  67
  68 void
  69 mcnemar_execute (const struct dataset *ds,
  70                   struct casereader *input,
  71                   enum mv_class exclude,
  72                   const struct npar_test *test,
  73                   bool exact UNUSED,
  74                   double timer UNUSED)
  75 {
  76   int i;
  77   bool warn = true;
  78
  79   const struct dictionary *dict = dataset_dict (ds);
  80
  81   const struct two_sample_test *t2s = UP_CAST (test, const struct two_sample_test, parent);
  82   struct ccase *c;
  83
  84   struct casereader *r = input;
  85
  86   struct mcnemar *mc = XCALLOC (t2s->n_pairs,  struct mcnemar);
  87
  88   for (i = 0 ; i < t2s->n_pairs; ++i)
  89     {
  90       mc[i].val0.f = mc[i].val1.f = SYSMIS;
  91     }
  92
  93   for (; (c = casereader_read (r)) != NULL; case_unref (c))
  94     {
  95       const double weight = dict_get_case_weight (dict, c, &warn);
  96
  97       for (i = 0 ; i < t2s->n_pairs; ++i)
  98         {
  99           variable_pair *vp = &t2s->pairs[i];
 100           const union value *value0 = case_data (c, (*vp)[0]);
 101           const union value *value1 = case_data (c, (*vp)[1]);
 102
 103           if (var_is_value_missing ((*vp)[0], value0) & exclude)
 104             continue;
 105
 106           if (var_is_value_missing ((*vp)[1], value1) & exclude)
 107             continue;
 108
 109
 110           if (mc[i].val0.f == SYSMIS)
 111             {
 112               if (mc[i].val1.f != value0->f)
 113                 mc[i].val0.f = value0->f;
 114               else if (mc[i].val1.f != value1->f)
 115                 mc[i].val0.f = value1->f;
 116             }
 117
 118           if (mc[i].val1.f == SYSMIS)
 119             {
 120               if (mc[i].val0.f != value1->f)
 121                 mc[i].val1.f = value1->f;
 122               else if (mc[i].val0.f != value0->f)
 123                 mc[i].val1.f = value0->f;
 124             }
 125
 126           if (mc[i].val0.f == value0->f && mc[i].val0.f == value1->f)
 127             {
 128               mc[i].n00 += weight;
 129             }
 130           else if (mc[i].val0.f == value0->f && mc[i].val1.f == value1->f)
 131             {
 132               mc[i].n10 += weight;
 133             }
 134           else if (mc[i].val1.f == value0->f && mc[i].val0.f == value1->f)
 135             {
 136               mc[i].n01 += weight;
 137             }
 138           else if (mc[i].val1.f == value0->f && mc[i].val1.f == value1->f)
 139             {
 140               mc[i].n11 += weight;
 141             }
 142           else
 143             {
 144               msg (ME, _("The McNemar test is appropriate only for dichotomous variables"));
 145             }
 146         }
 147     }
 148
 149   casereader_destroy (r);
 150
 151   for (i = 0 ; i < t2s->n_pairs ; ++i)
 152     output_freq_table (&t2s->pairs[i], mc + i, dict);
 153
 154   output_statistics_table (t2s, mc, dict);
 155
 156   free (mc);
 157 }
 158
 159 static char *
 160 make_pair_name (const variable_pair *pair)
 161 {
 162   return xasprintf ("%s & %s", var_to_string ((*pair)[0]),
 163                     var_to_string ((*pair)[1]));
 164 }
 165
 166 static void
 167 output_freq_table (variable_pair *vp,
 168                    const struct mcnemar *param,
 169                    const struct dictionary *dict)
 170 {
 171   struct pivot_table *table = pivot_table_create__ (
 172     pivot_value_new_user_text_nocopy (make_pair_name (vp)), "Frequencies");
 173   pivot_table_set_weight_var (table, dict_get_weight (dict));
 174
 175   struct pivot_dimension *vars[2];
 176   for (int i = 0; i < 2; i++)
 177     {
 178       vars[i] = pivot_dimension_create__ (
 179         table, i ? PIVOT_AXIS_COLUMN : PIVOT_AXIS_ROW,
 180         pivot_value_new_variable ((*vp)[i]));
 181       vars[i]->root->show_label = true;
 182
 183       for (int j = 0; j < 2; j++)
 184         {
 185           const union value *val = j ? &param->val1 : &param->val0;
 186           pivot_category_create_leaf_rc (
 187             vars[i]->root, pivot_value_new_var_value ((*vp)[0], val),
 188             PIVOT_RC_COUNT);
 189         }
 190     }
 191
 192   struct entry
 193     {
 194       int idx0;
 195       int idx1;
 196       double x;
 197     }
 198   entries[] = {
 199     { 0, 0, param->n00 },
 200     { 1, 0, param->n01 },
 201     { 0, 1, param->n10 },
 202     { 1, 1, param->n11 },
 203   };
 204   for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
 205     {
 206       const struct entry *e = &entries[i];
 207       pivot_table_put2 (table, e->idx0, e->idx1,
 208                         pivot_value_new_number (e->x));
 209     }
 210
 211   pivot_table_submit (table);
 212 }
 213
 214 static void
 215 output_statistics_table (const struct two_sample_test *t2s,
 216                          const struct mcnemar *mc,
 217                          const struct dictionary *dict)
 218 {
 219   struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
 220   pivot_table_set_weight_var (table, dict_get_weight (dict));
 221
 222   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 223                           N_("N"), PIVOT_RC_COUNT,
 224                           N_("Exact Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE,
 225                           N_("Exact Sig. (1-tailed)"), PIVOT_RC_SIGNIFICANCE,
 226                           N_("Point Probability"), PIVOT_RC_OTHER);
 227
 228   struct pivot_dimension *pairs = pivot_dimension_create (
 229     table, PIVOT_AXIS_ROW, N_("Pairs"));
 230
 231   for (size_t i = 0 ; i < t2s->n_pairs; ++i)
 232     {
 233       variable_pair *vp = &t2s->pairs[i];
 234       int pair_idx = pivot_category_create_leaf (
 235         pairs->root, pivot_value_new_user_text_nocopy (make_pair_name (vp)));
 236
 237       double n = mc[i].n00 + mc[i].n01 + mc[i].n10 + mc[i].n11;
 238       double sig = gsl_cdf_binomial_P ((mc[i].n01 > mc[i].n10) ? mc[i].n10: mc[i].n01,
 239                                        0.5, mc[i].n01 + mc[i].n10);
 240
 241       double point = gsl_ran_binomial_pdf (mc[i].n01, 0.5,
 242                                            mc[i].n01 + mc[i].n10);
 243       double entries[] = { n, 2.0 * sig, sig, point };
 244       for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 245         pivot_table_put2 (table, j, pair_idx,
 246                           pivot_value_new_number (entries[j]));
 247     }
 248
 249   pivot_table_submit (table);
 250 }