src/language/stats/cochran.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2010, 2011, 2014 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 <float.h>
  22 #include <gsl/gsl_cdf.h>
  23 #include <stdbool.h>
  24
  25 #include "data/casereader.h"
  26 #include "data/dataset.h"
  27 #include "data/dictionary.h"
  28 #include "data/format.h"
  29 #include "data/val-type.h"
  30 #include "data/variable.h"
  31 #include "language/stats/npar.h"
  32 #include "libpspp/cast.h"
  33 #include "libpspp/message.h"
  34 #include "libpspp/misc.h"
  35 #include "output/pivot-table.h"
  36
  37 #include "gettext.h"
  38 #define N_(msgid) msgid
  39 #define _(msgid) gettext (msgid)
  40
  41 struct cochran
  42 {
  43   double success;
  44   double failure;
  45
  46   double *hits;
  47   double *misses;
  48
  49   const struct dictionary *dict;
  50   double cc;
  51   double df;
  52   double q;
  53 };
  54
  55 static void show_freqs_box (const struct one_sample_test *ost, const struct cochran *ch);
  56 static void show_sig_box (const struct cochran *ch);
  57
  58 void
  59 cochran_execute (const struct dataset *ds,
  60               struct casereader *input,
  61               enum mv_class exclude,
  62               const struct npar_test *test,
  63               bool exact UNUSED, double timer UNUSED)
  64 {
  65   struct one_sample_test *ct = UP_CAST (test, struct one_sample_test, parent);
  66   int v;
  67   struct cochran ch;
  68   const struct dictionary *dict = dataset_dict (ds);
  69   const struct variable *weight = dict_get_weight (dict);
  70
  71   struct ccase *c;
  72   double rowsq = 0;
  73   ch.cc = 0.0;
  74   ch.dict = dict;
  75   ch.success = SYSMIS;
  76   ch.failure = SYSMIS;
  77   ch.hits = xcalloc (ct->n_vars, sizeof *ch.hits);
  78   ch.misses = xcalloc (ct->n_vars, sizeof *ch.misses);
  79
  80   for (; (c = casereader_read (input)); case_unref (c))
  81     {
  82       double case_hits = 0.0;
  83       const double w = weight ? case_num (c, weight) : 1.0;
  84       for (v = 0; v < ct->n_vars; ++v)
  85         {
  86           const struct variable *var = ct->vars[v];
  87           const union value *val = case_data (c, var);
  88
  89           if (var_is_value_missing (var, val) & exclude)
  90             continue;
  91
  92           if (ch.success == SYSMIS)
  93             {
  94               ch.success = val->f;
  95             }
  96           else if (ch.failure == SYSMIS && val->f != ch.success)
  97             {
  98               ch.failure = val->f;
  99             }
 100           if (ch.success == val->f)
 101             {
 102               ch.hits[v] += w;
 103               case_hits += w;
 104             }
 105           else if (ch.failure == val->f)
 106             {
 107               ch.misses[v] += w;
 108             }
 109           else
 110             {
 111               msg (MW, _("More than two values encountered.  Cochran Q test will not be run."));
 112               goto finish;
 113             }
 114         }
 115       ch.cc += w;
 116       rowsq += pow2 (case_hits);
 117     }
 118   casereader_destroy (input);
 119
 120   {
 121     double c_l = 0;
 122     double c_l2 = 0;
 123     for (v = 0; v < ct->n_vars; ++v)
 124       {
 125         c_l += ch.hits[v];
 126         c_l2 += pow2 (ch.hits[v]);
 127       }
 128
 129     ch.q = ct->n_vars * c_l2;
 130     ch.q -= pow2 (c_l);
 131     ch.q *= ct->n_vars - 1;
 132
 133     ch.q /= ct->n_vars * c_l - rowsq;
 134
 135     ch.df = ct->n_vars - 1;
 136   }
 137
 138   show_freqs_box (ct, &ch);
 139   show_sig_box (&ch);
 140
 141  finish:
 142
 143   free (ch.hits);
 144   free (ch.misses);
 145 }
 146
 147 static void
 148 show_freqs_box (const struct one_sample_test *ost, const struct cochran *ct)
 149 {
 150   struct pivot_table *table = pivot_table_create (N_("Frequencies"));
 151   pivot_table_set_weight_var (table, dict_get_weight (ct->dict));
 152
 153   char *success = xasprintf (_("Success (%.*g)"), DBL_DIG + 1, ct->success);
 154   char *failure = xasprintf (_("Failure (%.*g)"), DBL_DIG + 1, ct->failure);
 155   struct pivot_dimension *values = pivot_dimension_create (
 156     table, PIVOT_AXIS_COLUMN, N_("Value"),
 157     success, PIVOT_RC_COUNT,
 158     failure, PIVOT_RC_COUNT);
 159   values->root->show_label = true;
 160   free (failure);
 161   free (success);
 162
 163   struct pivot_dimension *variables = pivot_dimension_create (
 164     table, PIVOT_AXIS_ROW, N_("Variable"));
 165
 166   for (size_t i = 0 ; i < ost->n_vars ; ++i)
 167     {
 168       int row = pivot_category_create_leaf (
 169         variables->root, pivot_value_new_variable (ost->vars[i]));
 170
 171       pivot_table_put2 (table, 0, row, pivot_value_new_number (ct->hits[i]));
 172       pivot_table_put2 (table, 1, row, pivot_value_new_number (ct->misses[i]));
 173     }
 174
 175   pivot_table_submit (table);
 176 }
 177
 178 static void
 179 show_sig_box (const struct cochran *ch)
 180 {
 181   struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
 182
 183   pivot_table_set_weight_format (table, dict_get_weight_format (ch->dict));
 184
 185   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Value"), N_("Value"));
 186
 187   pivot_dimension_create (
 188     table, PIVOT_AXIS_ROW, N_("Statistics"),
 189     N_("N"), PIVOT_RC_COUNT,
 190     N_("Cochran's Q"), PIVOT_RC_SIGNIFICANCE,
 191     N_("df"), PIVOT_RC_INTEGER,
 192     N_("Asymp. Sig."), PIVOT_RC_SIGNIFICANCE);
 193
 194   double sig = gsl_cdf_chisq_Q (ch->q, ch->df);
 195   double entries[] = { ch->cc, ch->q, ch->df, sig };
 196   for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
 197     pivot_table_put2 (table, 0, i, pivot_value_new_number (entries[i]));
 198   pivot_table_submit (table);
 199 }