Fixed bug #22073
[pspp] / src / math / levene.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2004 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 #include "levene.h"
19 #include <libpspp/message.h>
20 #include <data/case.h>
21 #include <data/casereader.h>
22 #include <data/dictionary.h>
23 #include "group-proc.h"
24 #include <libpspp/hash.h>
25 #include <libpspp/str.h>
26 #include <data/variable.h>
27 #include <data/procedure.h>
28 #include <libpspp/misc.h>
29 #include "group.h"
30
31 #include <math.h>
32 #include <stdlib.h>
33
34 #include "xalloc.h"
35
36
37 /* This module calculates the Levene statistic for variables.
38
39    Just for reference, the Levene Statistic is a defines as follows:
40
41    W = \frac{ (n-k)\sum_{i=1}^k n_i(Z_{iL} - Z_{LL})^2}
42             { (k-1)\sum_{i=1}^k \sum_{j=1}^{n_i} (Z_{ij} - Z_{iL})^2}
43
44    where:
45         k is the number of groups
46         n is the total number of samples
47         n_i is the number of samples in the ith group
48         Z_{ij} is | Y_{ij} - Y_{iL} | where Y_{iL} is the mean of the ith group
49         Z_{iL} is the  mean of Z_{ij} over the ith group
50         Z_{LL} is the grand mean of Z_{ij}
51
52    Imagine calculating that with pencil and paper!
53
54  */
55
56
57 struct levene_info
58 {
59
60   /* Per group statistics */
61   struct t_test_proc **group_stats;
62
63   /* The independent variable */
64   const struct variable *v_indep;
65
66   /* Number of dependent variables */
67   size_t n_dep;
68
69   /* The dependent variables */
70   const struct variable  **v_dep;
71
72   /* Filter for missing values */
73   enum mv_class exclude;
74
75   /* An array of lz_stats for each variable */
76   struct lz_stats *lz;
77
78   /* The denominator for the expression for the Levene */
79   double *lz_denominator;
80
81 };
82
83 /* Per variable statistics */
84 struct lz_stats
85 {
86   /* Total of all lz */
87   double grand_total;
88
89   /* Mean of all lz */
90   double grand_mean;
91
92   /* The total number of cases */
93   double total_n ;
94
95   /* Number of groups */
96   int n_groups;
97 };
98
99 /* First pass */
100 static void  levene_precalc (const struct levene_info *l);
101 static int levene_calc (const struct dictionary *dict, const struct ccase *,
102                         const struct levene_info *l);
103 static void levene_postcalc (struct levene_info *);
104
105
106 /* Second pass */
107 static void levene2_precalc (struct levene_info *l);
108 static int levene2_calc (const struct dictionary *, const struct ccase *,
109                          struct levene_info *l);
110 static void levene2_postcalc (struct levene_info *);
111
112
113 void
114 levene(const struct dictionary *dict,
115        struct casereader *reader,
116        const struct variable *v_indep, size_t n_dep,
117        const struct variable **v_dep,
118        enum mv_class exclude)
119 {
120   struct casereader *pass1, *pass2;
121   struct ccase c;
122   struct levene_info l;
123
124   l.n_dep      = n_dep;
125   l.v_indep    = v_indep;
126   l.v_dep      = v_dep;
127   l.exclude    = exclude;
128   l.lz         = xnmalloc (l.n_dep, sizeof *l.lz);
129   l.lz_denominator = xnmalloc (l.n_dep, sizeof *l.lz_denominator);
130
131   casereader_split (reader, &pass1, &pass2);
132
133   levene_precalc (&l);
134   for (; casereader_read (pass1, &c); case_destroy (&c))
135     levene_calc (dict, &c, &l);
136   casereader_destroy (pass1);
137   levene_postcalc (&l);
138
139   levene2_precalc(&l);
140   for (; casereader_read (pass2, &c); case_destroy (&c))
141     levene2_calc (dict, &c, &l);
142   casereader_destroy (pass2);
143   levene2_postcalc (&l);
144
145   free (l.lz_denominator);
146   free (l.lz);
147 }
148
149 static void
150 levene_precalc (const struct levene_info *l)
151 {
152   size_t i;
153
154   for(i = 0; i < l->n_dep ; ++i )
155     {
156       const struct variable *var = l->v_dep[i];
157       struct group_proc *gp = group_proc_get (var);
158       struct group_statistics *gs;
159       struct hsh_iterator hi;
160
161       l->lz[i].grand_total = 0;
162       l->lz[i].total_n = 0;
163       l->lz[i].n_groups = gp->n_groups ;
164
165
166       for ( gs = hsh_first(gp->group_hash, &hi);
167             gs != 0;
168             gs = hsh_next(gp->group_hash, &hi))
169         {
170           gs->lz_total = 0;
171         }
172
173     }
174
175 }
176
177 static int
178 levene_calc (const struct dictionary *dict, const struct ccase *c,
179              const struct levene_info *l)
180 {
181   size_t i;
182   bool warn = false;
183   const union value *gv = case_data (c, l->v_indep);
184   struct group_statistics key;
185   double weight = dict_get_case_weight (dict, c, &warn);
186
187   key.id = *gv;
188
189   for (i = 0; i < l->n_dep; ++i)
190     {
191       const struct variable *var = l->v_dep[i];
192       struct group_proc *gp = group_proc_get (var);
193       double levene_z;
194       const union value *v = case_data (c, var);
195       struct group_statistics *gs;
196
197       gs = hsh_find(gp->group_hash,(void *) &key );
198
199       if ( 0 == gs )
200         continue ;
201
202       if ( !var_is_value_missing (var, v, l->exclude))
203         {
204           levene_z= fabs(v->f - gs->mean);
205           l->lz[i].grand_total += levene_z * weight;
206           l->lz[i].total_n += weight;
207
208           gs->lz_total += levene_z * weight;
209         }
210     }
211   return 0;
212 }
213
214
215 static void
216 levene_postcalc (struct levene_info *l)
217 {
218   size_t v;
219
220   for (v = 0; v < l->n_dep; ++v)
221     {
222       /* This is Z_LL */
223       l->lz[v].grand_mean = l->lz[v].grand_total / l->lz[v].total_n ;
224     }
225
226
227 }
228
229
230
231 static void
232 levene2_precalc (struct levene_info *l)
233 {
234   size_t v;
235
236
237   /* This stuff could go in the first post calc . . . */
238   for (v = 0;
239        v < l->n_dep;
240        ++v)
241     {
242       struct hsh_iterator hi;
243       struct group_statistics *g;
244
245       const struct variable *var = l->v_dep[v] ;
246       struct hsh_table *hash = group_proc_get (var)->group_hash;
247
248
249       for(g = (struct group_statistics *) hsh_first(hash,&hi);
250           g != 0 ;
251           g = (struct group_statistics *) hsh_next(hash,&hi) )
252         {
253           g->lz_mean = g->lz_total / g->n ;
254         }
255       l->lz_denominator[v] = 0;
256   }
257 }
258
259 static int
260 levene2_calc (const struct dictionary *dict, const struct ccase *c,
261               struct levene_info *l)
262 {
263   size_t i;
264   bool warn = false;
265
266   double weight = dict_get_case_weight (dict, c, &warn);
267
268   const union value *gv = case_data (c, l->v_indep);
269   struct group_statistics key;
270
271   key.id = *gv;
272
273   for (i = 0; i < l->n_dep; ++i)
274     {
275       double levene_z;
276       const struct variable *var = l->v_dep[i] ;
277       const union value *v = case_data (c, var);
278       struct group_statistics *gs;
279
280       gs = hsh_find(group_proc_get (var)->group_hash,(void *) &key );
281
282       if ( 0 == gs )
283         continue;
284
285       if ( !var_is_value_missing (var, v, l->exclude))
286         {
287           levene_z = fabs(v->f - gs->mean);
288           l->lz_denominator[i] += weight * pow2 (levene_z - gs->lz_mean);
289         }
290     }
291
292   return 0;
293 }
294
295
296 static void
297 levene2_postcalc (struct levene_info *l)
298 {
299   size_t v;
300
301   for (v = 0; v < l->n_dep; ++v)
302     {
303       double lz_numerator = 0;
304       struct hsh_iterator hi;
305       struct group_statistics *g;
306
307       const struct variable *var = l->v_dep[v] ;
308       struct group_proc *gp = group_proc_get (var);
309       struct hsh_table *hash = gp->group_hash;
310
311       for(g = (struct group_statistics *) hsh_first(hash,&hi);
312           g != 0 ;
313           g = (struct group_statistics *) hsh_next(hash,&hi) )
314         {
315           lz_numerator += g->n * pow2(g->lz_mean - l->lz[v].grand_mean );
316         }
317       lz_numerator *= ( gp->ugs.n - gp->n_groups );
318
319       l->lz_denominator[v] *= (gp->n_groups - 1);
320
321       gp->levene = lz_numerator / l->lz_denominator[v] ;
322
323     }
324 }
325