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
cochran_execute(const struct dataset * ds,struct casereader * input,enum mv_class exclude,const struct npar_test * test,bool exact UNUSED,double timer UNUSED)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_data (c, weight)->f: 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
show_freqs_box(const struct one_sample_test * ost,const struct cochran * ct)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
show_sig_box(const struct cochran * ch)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 }
200