GLM: Disable the interactions option until it actually gives the right results
[pspp] / src / language / stats / friedman.c
1 /* PSPP - a program for statistical analysis. -*-c-*-
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
18 #include <config.h>
19
20 #include "language/stats/friedman.h"
21
22 #include <gsl/gsl_cdf.h>
23 #include <math.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/variable.h"
30 #include "libpspp/message.h"
31 #include "libpspp/misc.h"
32 #include "output/tab.h"
33
34 #include "gettext.h"
35 #define _(msgid) gettext (msgid)
36
37
38 struct friedman
39 {
40   double *rank_sum;
41   double cc;
42   double chi_sq;
43   double w;
44   const struct dictionary *dict;
45 };
46
47 static void show_ranks_box (const struct one_sample_test *ost, 
48                             const struct friedman *fr);
49
50 static void show_sig_box (const struct one_sample_test *ost,
51                           const struct friedman *fr);
52
53 struct datum
54 {
55   long posn;
56   double x;
57 };
58
59 static int
60 cmp_x (const void *a_, const void *b_)
61 {
62   const struct datum *a = a_;
63   const struct datum *b = b_;
64
65   if (a->x < b->x)
66     return -1;
67   
68   return (a->x > b->x);
69 }
70
71 static int
72 cmp_posn (const void *a_, const void *b_)
73 {
74   const struct datum *a = a_;
75   const struct datum *b = b_;
76
77   if (a->posn < b->posn)
78     return -1;
79   
80   return (a->posn > b->posn);
81 }
82
83 void
84 friedman_execute (const struct dataset *ds,
85                   struct casereader *input,
86                   enum mv_class exclude,
87                   const struct npar_test *test,
88                   bool exact UNUSED,
89                   double timer UNUSED)
90 {
91   double numerator = 0.0;
92   double denominator = 0.0;
93   int v;
94   struct ccase *c;
95   const struct dictionary *dict = dataset_dict (ds);
96   const struct variable *weight = dict_get_weight (dict);
97
98   struct one_sample_test *ost = UP_CAST (test, struct one_sample_test, parent);
99   struct friedman_test *ft = UP_CAST (ost, struct friedman_test, parent);
100   bool warn = true;
101
102   double sigma_t = 0.0; 
103   struct datum *row = xcalloc (ost->n_vars, sizeof *row);
104   double rsq;
105   struct friedman fr;
106   fr.rank_sum = xcalloc (ost->n_vars, sizeof *fr.rank_sum);
107   fr.cc = 0.0;
108   fr.dict = dict;
109   for (v = 0; v < ost->n_vars; ++v)
110     {
111       row[v].posn = v;
112       fr.rank_sum[v] = 0.0;
113     }
114
115   input = casereader_create_filter_weight (input, dict, &warn, NULL);
116   input = casereader_create_filter_missing (input,
117                                             ost->vars, ost->n_vars,
118                                             exclude, 0, 0);
119
120   for (; (c = casereader_read (input)); case_unref (c))
121     {
122       double prev_x = SYSMIS;
123       int run_length = 0;
124
125       const double w = weight ? case_data (c, weight)->f: 1.0;
126
127       fr.cc += w;
128
129       for (v = 0; v < ost->n_vars; ++v)
130         {
131           const struct variable *var = ost->vars[v];
132           const union value *val = case_data (c, var);
133           row[v].x = val->f;
134         }
135
136       qsort (row, ost->n_vars, sizeof *row, cmp_x);
137       for (v = 0; v < ost->n_vars; ++v)
138         {
139           double x = row[v].x;
140           /* Replace value by the Rank */
141           if ( prev_x == x)
142             {
143               /* Deal with ties */
144               int i;
145               run_length++;
146               for (i = v - run_length; i < v; ++i)
147                 {
148                   row[i].x *= run_length ;
149                   row[i].x += v + 1;
150                   row[i].x /= run_length + 1;
151                 }
152               row[v].x = row[v-1].x;
153             }
154           else
155             {
156               row[v].x = v + 1;
157               if ( run_length > 0)
158                 {
159                   double t = run_length + 1;
160                   sigma_t += w * (pow3 (t) - t);
161                 }
162               run_length = 0;
163             }
164           prev_x = x;
165         }
166       if ( run_length > 0)
167         {
168           double t = run_length + 1;
169           sigma_t += w * (pow3 (t) - t );
170         }
171
172       qsort (row, ost->n_vars, sizeof *row, cmp_posn);
173
174       for (v = 0; v < ost->n_vars; ++v)
175         fr.rank_sum[v] += row[v].x * w;
176     }
177   casereader_destroy (input);
178   free (row);
179
180
181   for (v = 0; v < ost->n_vars; ++v)
182     {
183       numerator += pow2 (fr.rank_sum[v]);
184     }
185
186   rsq = numerator;
187
188   numerator *= 12.0 / (fr.cc * ost->n_vars * ( ost->n_vars + 1));
189   numerator -= 3 * fr.cc * ( ost->n_vars + 1);
190
191   denominator = 1 - sigma_t / ( fr.cc * ost->n_vars * ( pow2 (ost->n_vars) - 1));
192
193   fr.chi_sq = numerator / denominator;
194
195   if ( ft->kendalls_w)
196     {
197       fr.w = 12 * rsq ;
198       fr.w -= 3 * pow2 (fr.cc) *
199         ost->n_vars * pow2 (ost->n_vars + 1);
200
201       fr.w /= pow2 (fr.cc) * (pow3 (ost->n_vars) - ost->n_vars)
202         - fr.cc * sigma_t;
203     }
204   else
205     fr.w = SYSMIS;
206
207   show_ranks_box (ost, &fr);
208   show_sig_box (ost, &fr);
209
210   free (fr.rank_sum);
211 }
212
213 \f
214
215
216 static void
217 show_ranks_box (const struct one_sample_test *ost, const struct friedman *fr)
218 {
219   int i;
220   const int row_headers = 1;
221   const int column_headers = 1;
222   struct tab_table *table =
223     tab_create (row_headers + 1, column_headers + ost->n_vars);
224
225   tab_headers (table, row_headers, 0, column_headers, 0);
226
227   tab_title (table, _("Ranks"));
228
229   /* Vertical lines inside the box */
230   tab_box (table, 1, 0, -1, TAL_1,
231            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
232
233   /* Box around the table */
234   tab_box (table, TAL_2, TAL_2, -1, -1,
235            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
236
237
238   tab_text (table, 1, 0, 0, _("Mean Rank"));
239
240   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, column_headers);
241   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
242
243   for (i = 0 ; i < ost->n_vars ; ++i)
244     {
245       tab_text (table, 0, row_headers + i,
246                 TAB_LEFT, var_to_string (ost->vars[i]));
247
248       tab_double (table, 1, row_headers + i,
249                   0, fr->rank_sum[i] / fr->cc, 0);
250     }
251
252   tab_submit (table);
253 }
254
255
256 static void
257 show_sig_box (const struct one_sample_test *ost, const struct friedman *fr)
258 {
259   const struct friedman_test *ft = UP_CAST (ost, const struct friedman_test, parent);
260   
261   int row = 0;
262   const struct variable *weight = dict_get_weight (fr->dict);
263   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
264
265   const int row_headers = 1;
266   const int column_headers = 0;
267   struct tab_table *table =
268     tab_create (row_headers + 1, column_headers + (ft->kendalls_w ? 5 : 4));
269
270   tab_headers (table, row_headers, 0, column_headers, 0);
271
272   tab_title (table, _("Test Statistics"));
273
274   tab_text (table,  0, column_headers + row++,
275             TAT_TITLE | TAB_LEFT , _("N"));
276
277   if ( ft->kendalls_w)
278     tab_text (table,  0, column_headers + row++,
279               TAT_TITLE | TAB_LEFT , _("Kendall's W"));
280
281   tab_text (table,  0, column_headers + row++,
282             TAT_TITLE | TAB_LEFT , _("Chi-Square"));
283
284   tab_text (table,  0, column_headers + row++,
285             TAT_TITLE | TAB_LEFT, _("df"));
286
287   tab_text (table,  0, column_headers + row++,
288             TAT_TITLE | TAB_LEFT, _("Asymp. Sig."));
289
290   /* Box around the table */
291   tab_box (table, TAL_2, TAL_2, -1, -1,
292            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
293
294
295   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
296   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
297
298   row = 0;
299   tab_double (table, 1, column_headers + row++, 
300               0, fr->cc, wfmt);
301
302   if (ft->kendalls_w)
303     tab_double (table, 1, column_headers + row++, 
304                 0, fr->w, 0);
305
306   tab_double (table, 1, column_headers + row++, 
307               0, fr->chi_sq, 0);
308
309   tab_double (table, 1, column_headers + row++, 
310               0, ost->n_vars - 1, &F_8_0);
311
312   tab_double (table, 1, column_headers + row++, 
313               0, gsl_cdf_chisq_Q (fr->chi_sq, ost->n_vars - 1), 
314               0);
315
316   tab_submit (table);
317 }