QUICK CLUSTER: Implement pairwise missing option
[pspp] / src / language / stats / t-test-paired.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 <math.h>
20 #include <gsl/gsl_cdf.h>
21
22 #include "t-test.h"
23
24 #include "math/moments.h"
25 #include "math/correlation.h"
26 #include "data/casereader.h"
27 #include "data/dictionary.h"
28 #include "data/format.h"
29 #include "data/variable.h"
30 #include "libpspp/hmapx.h"
31 #include "libpspp/hash-functions.h"
32
33 #include "output/tab.h"
34
35 #include "gettext.h"
36 #define _(msgid) gettext (msgid)
37
38
39 struct pair_stats
40 {
41   int posn;
42   double sum_of_prod;
43   struct moments *mom0;
44   const struct variable *var0;
45
46   struct moments *mom1;
47   const struct variable *var1;
48
49   struct moments *mom_diff;
50 };
51
52 struct paired_samp
53 {
54   struct hmapx hmap;
55 };
56
57 static void paired_summary (const struct tt *tt, struct paired_samp *os);
58 static void paired_correlations (const struct tt *tt, struct paired_samp *os);
59 static void paired_test (const struct tt *tt, const struct paired_samp *os);
60
61 void
62 paired_run (const struct tt *tt, size_t n_pairs, vp *pairs, struct casereader *reader)
63 {
64   int i;
65   struct ccase *c;
66   struct paired_samp ps;
67   struct casereader *r;
68   struct hmapx_node *node;
69   struct pair_stats *pp = NULL;
70
71   hmapx_init (&ps.hmap);
72
73   for (i = 0; i < n_pairs; ++i)
74     {
75       vp *pair = &pairs[i];
76       unsigned int hash;
77       struct pair_stats *pp = xzalloc (sizeof *pp);
78       pp->posn = i;
79       pp->var0 = (*pair)[0];
80       pp->var1 = (*pair)[1];
81       pp->mom0 = moments_create (MOMENT_VARIANCE);
82       pp->mom1 = moments_create (MOMENT_VARIANCE);
83       pp->mom_diff = moments_create (MOMENT_VARIANCE);
84
85       hash = hash_pointer ((*pair)[0], 0);
86       hash = hash_pointer ((*pair)[1], hash);
87
88       hmapx_insert (&ps.hmap, pp, hash);
89     }
90
91   r = casereader_clone (reader);
92   for ( ; (c = casereader_read (r) ); case_unref (c))
93     {
94       double w = dict_get_case_weight (tt->dict, c, NULL);
95
96       struct hmapx_node *node;
97       struct pair_stats *pp = NULL;
98       HMAPX_FOR_EACH (pp, node, &ps.hmap)
99         {
100           const union value *val0 = case_data (c, pp->var0);
101           const union value *val1 = case_data (c, pp->var1);
102           if (var_is_value_missing (pp->var0, val0, tt->exclude))
103             continue;
104
105           if (var_is_value_missing (pp->var1, val1, tt->exclude))
106             continue;
107
108           moments_pass_one (pp->mom0, val0->f, w);
109           moments_pass_one (pp->mom1, val1->f, w);
110           moments_pass_one (pp->mom_diff, val0->f - val1->f, w);
111         }
112     }
113   casereader_destroy (r);
114
115   r = reader;
116   for ( ; (c = casereader_read (r) ); case_unref (c))
117     {
118       double w = dict_get_case_weight (tt->dict, c, NULL);
119
120       struct hmapx_node *node;
121       struct pair_stats *pp = NULL;
122       HMAPX_FOR_EACH (pp, node, &ps.hmap)
123         {
124           const union value *val0 = case_data (c, pp->var0);
125           const union value *val1 = case_data (c, pp->var1);
126           if (var_is_value_missing (pp->var0, val0, tt->exclude))
127             continue;
128
129           if (var_is_value_missing (pp->var1, val1, tt->exclude))
130             continue;
131
132           moments_pass_two (pp->mom0, val0->f, w);
133           moments_pass_two (pp->mom1, val1->f, w);
134           moments_pass_two (pp->mom_diff, val0->f - val1->f, w);
135           pp->sum_of_prod += val0->f * val1->f;
136         }
137     }
138   casereader_destroy (r);
139
140   paired_summary (tt, &ps);
141   paired_correlations (tt, &ps);
142   paired_test (tt, &ps);
143
144   /* Clean up */
145   HMAPX_FOR_EACH (pp, node, &ps.hmap)
146     {
147       moments_destroy (pp->mom0);
148       moments_destroy (pp->mom1);
149       moments_destroy (pp->mom_diff);
150       free (pp);
151     }
152
153   hmapx_destroy (&ps.hmap);
154 }
155
156 static void
157 paired_summary (const struct tt *tt, struct paired_samp *os)
158 {
159   size_t n_pairs = hmapx_count (&os->hmap);
160   struct hmapx_node *node;
161   struct pair_stats *pp = NULL;
162
163   const int heading_rows = 1;
164   const int heading_cols = 2;
165
166   const int cols = 4 + heading_cols;
167   const int rows = n_pairs * 2 + heading_rows;
168   struct tab_table *t = tab_create (cols, rows);
169   const struct fmt_spec *wfmt = tt->wv ? var_get_print_format (tt->wv) : & F_8_0;
170
171   tab_headers (t, 0, 0, heading_rows, 0);
172   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols - 1, rows - 1);
173   tab_box (t, -1, -1,       TAL_0, TAL_1, heading_cols, 0, cols - 1, rows - 1);
174
175   tab_hline (t, TAL_2, 0, cols - 1, 1);
176
177   tab_title (t, _("Paired Sample Statistics"));
178   tab_vline (t, TAL_2, heading_cols, 0, rows - 1);
179   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("N"));
180   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
181   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
182   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("S.E. Mean"));
183
184   HMAPX_FOR_EACH (pp, node, &os->hmap)
185     {
186       int v = pp->posn;
187       double cc, mean, sigma;
188
189       tab_text_format (t, 0, v * 2 + heading_rows, TAB_LEFT, _("Pair %d"), pp->posn);
190
191       /* first var */
192       moments_calculate (pp->mom0, &cc, &mean, &sigma, NULL, NULL);
193       tab_text (t, 1, v * 2 + heading_rows, TAB_LEFT, var_to_string (pp->var0));
194       tab_double (t, 3, v * 2 + heading_rows, TAB_RIGHT, cc, wfmt);
195       tab_double (t, 2, v * 2 + heading_rows, TAB_RIGHT, mean, NULL);
196       tab_double (t, 4, v * 2 + heading_rows, TAB_RIGHT, sqrt (sigma), NULL);
197       tab_double (t, 5, v * 2 + heading_rows, TAB_RIGHT, sqrt (sigma / cc), NULL);
198
199       /* second var */
200       moments_calculate (pp->mom1, &cc, &mean, &sigma, NULL, NULL);
201       tab_text (t, 1, v * 2 + 1 + heading_rows, TAB_LEFT, var_to_string (pp->var1));      
202       tab_double (t, 3, v * 2 + 1 + heading_rows, TAB_RIGHT, cc, wfmt);
203       tab_double (t, 2, v * 2 + 1 + heading_rows, TAB_RIGHT, mean, NULL);
204       tab_double (t, 4, v * 2 + 1 + heading_rows, TAB_RIGHT, sqrt (sigma), NULL);
205       tab_double (t, 5, v * 2 + 1 + heading_rows, TAB_RIGHT, sqrt (sigma / cc), NULL);
206     }
207
208   tab_submit (t);
209 }
210
211
212 static void
213 paired_correlations (const struct tt *tt, struct paired_samp *os)
214 {
215   size_t n_pairs = hmapx_count (&os->hmap);
216   struct hmapx_node *node;
217   struct pair_stats *pp = NULL;
218   const int heading_rows = 1;
219   const int heading_cols = 2;
220
221   const int cols = 5;
222   const int rows = n_pairs + heading_rows;
223   struct tab_table *t = tab_create (cols, rows);
224   const struct fmt_spec *wfmt = tt->wv ? var_get_print_format (tt->wv) : & F_8_0;
225
226   tab_headers (t, 0, 0, heading_rows, 0);
227   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols - 1, rows - 1);
228
229   tab_hline (t, TAL_2, 0, cols - 1, 1);
230
231   tab_title (t, _("Paired Samples Correlations"));
232   tab_vline (t, TAL_2, heading_cols, 0, rows - 1);
233   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("N"));
234   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Correlation"));
235   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
236
237   HMAPX_FOR_EACH (pp, node, &os->hmap)
238     {
239       double corr;
240       double cc0, mean0, sigma0;
241       double cc1, mean1, sigma1;
242       int v = pp->posn;
243
244       tab_text_format (t, 0, v + heading_rows, TAB_LEFT, _("Pair %d"), pp->posn);
245
246       tab_text_format (t, 1, v + heading_rows, TAB_LEFT, _("%s & %s"), 
247                        var_to_string (pp->var0),
248                        var_to_string (pp->var1));
249
250       moments_calculate (pp->mom0, &cc0, &mean0, &sigma0, NULL, NULL);
251       moments_calculate (pp->mom1, &cc1, &mean1, &sigma1, NULL, NULL);
252
253       /* If this fails, then we're not dealing with missing values properly */
254       assert (cc0 == cc1);
255
256       tab_double (t, 2, v + heading_rows, TAB_RIGHT, cc0, wfmt);
257
258       corr = pp->sum_of_prod / cc0  - (mean0 * mean1);
259       corr /= sqrt (sigma0 * sigma1);
260       corr *= cc0 / (cc0 - 1);
261
262       tab_double (t, 3, v + heading_rows, TAB_RIGHT, corr, NULL);
263       tab_double (t, 4, v + heading_rows, TAB_RIGHT, 
264                   2.0 * significance_of_correlation (corr, cc0), NULL);
265     }
266
267   tab_submit (t);
268 }
269
270
271 static void
272 paired_test (const struct tt *tt, const struct paired_samp *os)
273 {
274   size_t n_pairs = hmapx_count (&os->hmap);
275   struct hmapx_node *node;
276   struct pair_stats *pp = NULL;
277
278   const int heading_rows = 3;
279   const int heading_cols = 2;
280   const size_t rows = heading_rows + n_pairs;
281   const size_t cols = 10;
282   const struct fmt_spec *wfmt = tt->wv ? var_get_print_format (tt->wv) : & F_8_0;
283
284   struct tab_table *t = tab_create (cols, rows);
285
286   tab_headers (t, 0, 0, heading_rows, 0);
287   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols - 1, rows - 1);
288   tab_hline (t, TAL_2, 0, cols - 1, 3);
289
290   tab_title (t, _("Paired Samples Test"));
291   tab_hline (t, TAL_1, heading_cols, 6, 1);
292   tab_vline (t, TAL_2, heading_cols, 0, rows - 1);
293
294   tab_box (t, -1, -1, -1, TAL_1, heading_cols, 0, cols - 1, rows - 1);
295
296   tab_joint_text (t, 2, 0, 6, 0, TAB_CENTER,
297                   _("Paired Differences"));
298
299   tab_joint_text_format (t, 5, 1, 6, 1, TAB_CENTER,
300                          _("%g%% Confidence Interval of the Difference"),
301                          tt->confidence * 100.0);
302
303   tab_vline (t, TAL_GAP, 6, 1, 1);
304   tab_hline (t, TAL_1, 5, 6, 2);
305   tab_text (t, 7, 2, TAB_CENTER | TAT_TITLE, _("t"));
306   tab_text (t, 8, 2, TAB_CENTER | TAT_TITLE, _("df"));
307   tab_text (t, 9, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
308   tab_text (t, 4, 2, TAB_CENTER | TAT_TITLE, _("Std. Error Mean"));
309   tab_text (t, 3, 2, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
310   tab_text (t, 2, 2, TAB_CENTER | TAT_TITLE, _("Mean"));
311
312   tab_text (t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
313   tab_text (t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
314
315   HMAPX_FOR_EACH (pp, node, &os->hmap)
316     {
317       int v = pp->posn;
318       double cc, mean, sigma;
319       double df ;
320       double tval;
321       double p, q;
322       double se_mean;
323
324       moments_calculate (pp->mom_diff, &cc, &mean, &sigma, NULL, NULL);
325
326       df = cc - 1.0;
327       tab_text_format (t, 0, v + heading_rows, TAB_LEFT, _("Pair %d"), v);
328
329       tab_text_format (t, 1, v + heading_rows, TAB_LEFT, _("%s - %s"), 
330                        var_to_string (pp->var0),
331                        var_to_string (pp->var1));
332
333       tval = mean * sqrt (cc / sigma);
334       se_mean = sqrt (sigma / cc);
335
336       tab_double (t, 2, v + heading_rows, TAB_RIGHT, mean, NULL);
337       tab_double (t, 3, v + heading_rows, TAB_RIGHT, sqrt (sigma), NULL);
338       tab_double (t, 4, v + heading_rows, TAB_RIGHT, se_mean, NULL);
339
340       tab_double (t, 7, v + heading_rows, TAB_RIGHT, tval, NULL);
341       tab_double (t, 8, v + heading_rows, TAB_RIGHT, df, wfmt);
342
343
344       p = gsl_cdf_tdist_P (tval, df);
345       q = gsl_cdf_tdist_Q (tval, df);
346
347       tab_double (t, 9, v + heading_rows, TAB_RIGHT, 2.0 * (tval > 0 ? q : p), NULL);
348
349       tval = gsl_cdf_tdist_Qinv ( (1.0 - tt->confidence) / 2.0, df);
350
351       tab_double (t, 5, v + heading_rows, TAB_RIGHT, mean - tval * se_mean, NULL);
352       tab_double (t, 6, v + heading_rows, TAB_RIGHT, mean + tval * se_mean, NULL);
353     }
354
355   tab_submit (t);
356 }