src/language/stats/cochran.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2010, 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 "language/stats/cochran.h"
  20
  21 #include <gsl/gsl_cdf.h>
  22 #include <stdbool.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/val-type.h"
  29 #include "data/variable.h"
  30 #include "language/stats/npar.h"
  31 #include "libpspp/cast.h"
  32 #include "libpspp/message.h"
  33 #include "libpspp/misc.h"
  34 #include "output/tab.h"
  35
  36 #include "gettext.h"
  37 #define _(msgid) gettext (msgid)
  38
  39 struct cochran
  40 {
  41   double success;
  42   double failure;
  43
  44   double *hits;
  45   double *misses;
  46
  47   const struct dictionary *dict;
  48   double cc;
  49   double df;
  50   double q;
  51 };
  52
  53 static void show_freqs_box (const struct one_sample_test *ost, const struct cochran *ch);
  54 static void show_sig_box (const struct cochran *ch);
  55
  56 void
  57 cochran_execute (const struct dataset *ds,
  58               struct casereader *input,
  59               enum mv_class exclude,
  60               const struct npar_test *test,
  61               bool exact UNUSED, double timer UNUSED)
  62 {
  63   struct one_sample_test *ct = UP_CAST (test, struct one_sample_test, parent);
  64   int v;
  65   struct cochran ch;
  66   const struct dictionary *dict = dataset_dict (ds);
  67   const struct variable *weight = dict_get_weight (dict);
  68
  69   struct ccase *c;
  70   double rowsq = 0;
  71   ch.cc = 0.0;
  72   ch.dict = dict;
  73   ch.success = SYSMIS;
  74   ch.failure = SYSMIS;
  75   ch.hits = xcalloc (ct->n_vars, sizeof *ch.hits);
  76   ch.misses = xcalloc (ct->n_vars, sizeof *ch.misses);
  77
  78   for (; (c = casereader_read (input)); case_unref (c))
  79     {
  80       double case_hits = 0.0;
  81       const double w = weight ? case_data (c, weight)->f: 1.0;
  82       for (v = 0; v < ct->n_vars; ++v)
  83         {
  84           const struct variable *var = ct->vars[v];
  85           const union value *val = case_data (c, var);
  86
  87           if ( var_is_value_missing (var, val, exclude))
  88             continue;
  89
  90           if ( ch.success == SYSMIS)
  91             {
  92               ch.success = val->f;
  93             }
  94           else if (ch.failure == SYSMIS && val->f != ch.success)
  95             {
  96               ch.failure = val->f;
  97             }
  98           if ( ch.success == val->f)
  99             {
 100               ch.hits[v] += w;
 101               case_hits += w;
 102             }
 103           else if ( ch.failure == val->f)
 104             {
 105               ch.misses[v] += w;
 106             }
 107           else
 108             {
 109               msg (MW, _("More than two values encountered.  Cochran Q test will not be run."));
 110               goto finish;
 111             }
 112         }
 113       ch.cc += w;
 114       rowsq += pow2 (case_hits);
 115     }
 116   casereader_destroy (input);
 117
 118   {
 119     double c_l = 0;
 120     double c_l2 = 0;
 121     for (v = 0; v < ct->n_vars; ++v)
 122       {
 123         c_l += ch.hits[v];
 124         c_l2 += pow2 (ch.hits[v]);
 125       }
 126
 127     ch.q = ct->n_vars * c_l2;
 128     ch.q -= pow2 (c_l);
 129     ch.q *= ct->n_vars - 1;
 130
 131     ch.q /= ct->n_vars * c_l - rowsq;
 132
 133     ch.df = ct->n_vars - 1;
 134   }
 135
 136   show_freqs_box (ct, &ch);
 137   show_sig_box (&ch);
 138
 139  finish:
 140
 141   free (ch.hits);
 142   free (ch.misses);
 143 }
 144
 145 static void
 146 show_freqs_box (const struct one_sample_test *ost, const struct cochran *ct)
 147 {
 148   int i;
 149   const struct variable *weight = dict_get_weight (ct->dict);
 150   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 151
 152   const int row_headers = 1;
 153   const int column_headers = 2;
 154   struct tab_table *table =
 155     tab_create (row_headers + 2, column_headers + ost->n_vars);
 156
 157   tab_headers (table, row_headers, 0, column_headers, 0);
 158
 159   tab_title (table, _("Frequencies"));
 160
 161   /* Vertical lines inside the box */
 162   tab_box (table, 1, 0, -1, TAL_1,
 163            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
 164
 165   /* Box around the table */
 166   tab_box (table, TAL_2, TAL_2, -1, -1,
 167            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 168
 169   tab_joint_text (table, 1, 0, 2, 0,
 170                   TAT_TITLE | TAB_CENTER, _("Value"));
 171
 172   tab_text_format (table, 1, 1, 0, _("Success (%g)"), ct->success);
 173   tab_text_format (table, 2, 1, 0, _("Failure (%g)"), ct->failure);
 174
 175   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, column_headers);
 176   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 177
 178   for (i = 0 ; i < ost->n_vars ; ++i)
 179     {
 180       tab_text (table, 0, column_headers + i,
 181                 TAB_LEFT, var_to_string (ost->vars[i]));
 182
 183       tab_double (table, 1, column_headers + i, 0,
 184                   ct->hits[i], wfmt);
 185
 186       tab_double (table, 2, column_headers + i, 0,
 187                   ct->misses[i], wfmt);
 188     }
 189
 190   tab_submit (table);
 191 }
 192
 193
 194
 195 static void
 196 show_sig_box (const struct cochran *ch)
 197 {
 198   const struct variable *weight = dict_get_weight (ch->dict);
 199   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 200
 201   const int row_headers = 1;
 202   const int column_headers = 0;
 203   struct tab_table *table =
 204     tab_create (row_headers + 1, column_headers + 4);
 205
 206   tab_headers (table, row_headers, 0, column_headers, 0);
 207
 208   tab_title (table, _("Test Statistics"));
 209
 210   tab_text (table,  0, column_headers,
 211             TAT_TITLE | TAB_LEFT , _("N"));
 212
 213   tab_text (table,  0, 1 + column_headers,
 214             TAT_TITLE | TAB_LEFT , _("Cochran's Q"));
 215
 216   tab_text (table,  0, 2 + column_headers,
 217             TAT_TITLE | TAB_LEFT, _("df"));
 218
 219   tab_text (table,  0, 3 + column_headers,
 220             TAT_TITLE | TAB_LEFT, _("Asymp. Sig."));
 221
 222   /* Box around the table */
 223   tab_box (table, TAL_2, TAL_2, -1, -1,
 224            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 225
 226   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 227   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 228
 229   tab_double (table, 1, column_headers,
 230               0, ch->cc, wfmt);
 231
 232   tab_double (table, 1, column_headers + 1,
 233               0, ch->q, 0);
 234
 235   tab_double (table, 1, column_headers + 2,
 236               0, ch->df, &F_8_0);
 237
 238   tab_double (table, 1, column_headers + 3,
 239               0, gsl_cdf_chisq_Q (ch->q, ch->df),
 240               0);
 241
 242   tab_submit (table);
 243 }