0e36b87a8c12b59817c944be0e14bf14430cb2f6
[pspp] / src / levene.c
1 /* This file is part of GNU PSPP 
2    Computes Levene test  statistic.
3
4    Copyright (C) 2004 Free Software Foundation, Inc.
5    Written by John Darrington <john@darrington.wattle.id.au>
6
7    This program is free software; you can redistribute it and/or
8    modify it under the terms of the GNU General Public License as
9    published by the Free Software Foundation; either version 2 of the
10    License, or (at your option) any later version.
11
12    This program is distributed in the hope that it will be useful, but
13    WITHOUT ANY WARRANTY; without even the implied warranty of
14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15    General Public License for more details.
16
17    You should have received a copy of the GNU General Public License
18    along with this program; if not, write to the Free Software
19    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
20    02111-1307, USA. */
21
22 #include <config.h>
23 #include "levene.h"
24 #include "error.h"
25 #include "case.h"
26 #include "casefile.h"
27 #include "hash.h"
28 #include "str.h"
29 #include "var.h"
30 #include "vfm.h"
31 #include "alloc.h"
32 #include "misc.h"
33 #include "group.h"
34
35 #include <math.h>
36 #include <stdlib.h>
37
38
39 /* This module calculates the Levene statistic for variables.
40
41    Just for reference, the Levene Statistic is a defines as follows:
42
43    W = \frac{ (n-k)\sum_{i=1}^k n_i(Z_{iL} - Z_{LL})^2}
44             { (k-1)\sum_{i=1}^k \sum_{j=1}^{n_i} (Z_{ij} - Z_{iL})^2}
45
46    where:
47         k is the number of groups
48         n is the total number of samples
49         n_i is the number of samples in the ith group
50         Z_{ij} is | Y_{ij} - Y_{iL} | where Y_{iL} is the mean of the ith group
51         Z_{iL} is the  mean of Z_{ij} over the ith group
52         Z_{LL} is the grand mean of Z_{ij}
53
54    Imagine calculating that with pencil and paper!
55
56  */
57
58
59 struct levene_info
60 {
61
62   /* Per group statistics */
63   struct t_test_proc **group_stats;
64
65   /* The independent variable */
66   struct variable *v_indep; 
67
68   /* Number of dependent variables */
69   int n_dep;
70
71   /* The dependent variables */
72   struct variable  **v_dep;
73
74   /* How to treat missing values */
75   enum lev_missing missing;
76
77   /* Function to test for missing values */
78   is_missing_func is_missing;
79
80 };
81
82 /* First pass */
83 static void  levene_precalc (const struct levene_info *l);
84 static int levene_calc (const struct ccase *, void *);
85 static void levene_postcalc (void *);
86
87
88 /* Second pass */
89 static void levene2_precalc (void *);
90 static int levene2_calc (const struct ccase *, void *);
91 static void levene2_postcalc (void *);
92
93
94 void  
95 levene(const struct casefile *cf,
96        struct variable *v_indep, int n_dep, struct variable **v_dep,
97              enum lev_missing missing,   is_missing_func value_is_missing)
98 {
99   struct casereader *r;
100   struct ccase c;
101   struct levene_info l;
102
103   l.n_dep      = n_dep;
104   l.v_indep    = v_indep;
105   l.v_dep      = v_dep;
106   l.missing    = missing;
107   l.is_missing = value_is_missing;
108
109
110
111   levene_precalc(&l);
112   for(r = casefile_get_reader (cf);
113       casereader_read (r, &c) ;
114       case_destroy (&c)) 
115     {
116       levene_calc(&c,&l);
117     }
118   casereader_destroy (r);
119   levene_postcalc(&l);
120
121   levene2_precalc(&l);
122   for(r = casefile_get_reader (cf);
123       casereader_read (r, &c) ;
124       case_destroy (&c)) 
125     {
126       levene2_calc(&c,&l);
127     }
128   casereader_destroy (r);
129   levene2_postcalc(&l);
130
131 }
132
133 /* Internal variables used in calculating the Levene statistic */
134
135 /* Per variable statistics */
136 struct lz_stats
137 {
138   /* Total of all lz */
139   double grand_total;
140
141   /* Mean of all lz */
142   double grand_mean;
143
144   /* The total number of cases */
145   double total_n ; 
146
147   /* Number of groups */
148   int n_groups;
149 };
150
151 /* An array of lz_stats for each variable */
152 static struct lz_stats *lz;
153
154
155 static void 
156 levene_precalc (const struct levene_info *l)
157 {
158   int i;
159
160   lz  = xmalloc (sizeof (struct lz_stats ) * l->n_dep ) ;
161
162   for(i = 0; i < l->n_dep ; ++i ) 
163     {
164       struct variable *var = l->v_dep[i];
165       struct group_statistics *gs;
166       struct hsh_iterator hi;
167
168       lz[i].grand_total = 0;
169       lz[i].total_n = 0;
170       lz[i].n_groups = var->p.grp_data.n_groups ; 
171
172       
173       for ( gs = hsh_first(var->p.grp_data.group_hash, &hi);
174             gs != 0;
175             gs = hsh_next(var->p.grp_data.group_hash, &hi))
176         {
177           gs->lz_total = 0;
178         }
179             
180     }
181
182 }
183
184 static int 
185 levene_calc (const struct ccase *c, void *_l)
186 {
187   int i;
188   int warn = 0;
189   struct levene_info *l = (struct levene_info *) _l;
190   const union value *gv = case_data (c, l->v_indep->fv);
191   struct group_statistics key;
192   double weight = dict_get_case_weight(default_dict,c,&warn); 
193
194   /* Skip the entire case if /MISSING=LISTWISE is set */
195   if ( l->missing == LEV_LISTWISE ) 
196     {
197       for (i = 0; i < l->n_dep; ++i) 
198         {
199           struct variable *v = l->v_dep[i];
200           const union value *val = case_data (c, v->fv);
201
202           if (l->is_missing(val,v) )
203             {
204               return 0;
205             }
206         }
207     }
208
209   
210   key.id = *gv;
211
212   for (i = 0; i < l->n_dep; ++i) 
213     {
214       struct variable *var = l->v_dep[i];
215       double levene_z;
216       const union value *v = case_data (c, var->fv);
217       struct group_statistics *gs;
218
219       gs = hsh_find(var->p.grp_data.group_hash,(void *) &key );
220
221       if ( 0 == gs ) 
222         continue ;
223
224       if ( ! l->is_missing(v,var))
225         {
226           levene_z= fabs(v->f - gs->mean);
227           lz[i].grand_total += levene_z * weight;
228           lz[i].total_n += weight; 
229
230           gs->lz_total += levene_z * weight;
231         }
232
233     }
234   return 0;
235 }
236
237
238 static void 
239 levene_postcalc (void *_l)
240 {
241   int v;
242
243   struct levene_info *l = (struct levene_info *) _l;
244
245   for (v = 0; v < l->n_dep; ++v) 
246     {
247       /* This is Z_LL */
248       lz[v].grand_mean = lz[v].grand_total / lz[v].total_n ;
249     }
250
251   
252 }
253
254
255 /* The denominator for the expression for the Levene */
256 static double *lz_denominator;
257
258 static void 
259 levene2_precalc (void *_l)
260 {
261   int v;
262
263   struct levene_info *l = (struct levene_info *) _l;
264
265   lz_denominator = (double *) xmalloc(sizeof(double) * l->n_dep);
266
267   /* This stuff could go in the first post calc . . . */
268   for (v = 0; v < l->n_dep; ++v) 
269     {
270       struct hsh_iterator hi;
271       struct group_statistics *g;
272
273       struct variable *var = l->v_dep[v] ;
274       struct hsh_table *hash = var->p.grp_data.group_hash;
275
276
277       for(g = (struct group_statistics *) hsh_first(hash,&hi);
278           g != 0 ;
279           g = (struct group_statistics *) hsh_next(hash,&hi) )
280         {
281           g->lz_mean = g->lz_total / g->n ;
282         }
283       lz_denominator[v] = 0;
284   }
285 }
286
287 static int 
288 levene2_calc (const struct ccase *c, void *_l)
289 {
290   int i;
291   int warn = 0;
292
293   struct levene_info *l = (struct levene_info *) _l;
294
295   double weight = dict_get_case_weight(default_dict,c,&warn); 
296
297   const union value *gv = case_data (c, l->v_indep->fv);
298   struct group_statistics key;
299
300   /* Skip the entire case if /MISSING=LISTWISE is set */
301   if ( l->missing == LEV_LISTWISE ) 
302     {
303       for (i = 0; i < l->n_dep; ++i) 
304         {
305           struct variable *v = l->v_dep[i];
306           const union value *val = case_data (c, v->fv);
307
308           if (l->is_missing(val,v) )
309             {
310               return 0;
311             }
312         }
313     }
314
315   key.id = *gv;
316
317   for (i = 0; i < l->n_dep; ++i) 
318     {
319       double levene_z;
320       struct variable *var = l->v_dep[i] ;
321       const union value *v = case_data (c, var->fv);
322       struct group_statistics *gs;
323
324       gs = hsh_find(var->p.grp_data.group_hash,(void *) &key );
325
326       if ( 0 == gs ) 
327         continue;
328
329       if ( ! l->is_missing(v,var) )
330         {
331           levene_z = fabs(v->f - gs->mean); 
332           lz_denominator[i] += weight * pow2(levene_z - gs->lz_mean);
333         }
334     }
335
336   return 0;
337 }
338
339
340 static void 
341 levene2_postcalc (void *_l)
342 {
343   int v;
344
345   struct levene_info *l = (struct levene_info *) _l;
346
347   for (v = 0; v < l->n_dep; ++v) 
348     {
349       double lz_numerator = 0;
350       struct hsh_iterator hi;
351       struct group_statistics *g;
352
353       struct variable *var = l->v_dep[v] ;
354       struct hsh_table *hash = var->p.grp_data.group_hash;
355
356       for(g = (struct group_statistics *) hsh_first(hash,&hi);
357           g != 0 ;
358           g = (struct group_statistics *) hsh_next(hash,&hi) )
359         {
360           lz_numerator += g->n * pow2(g->lz_mean - lz[v].grand_mean );
361         }
362       lz_numerator *= ( l->v_dep[v]->p.grp_data.ugs.n - 
363                         l->v_dep[v]->p.grp_data.n_groups );
364
365       lz_denominator[v] *= (l->v_dep[v]->p.grp_data.n_groups - 1);
366
367       l->v_dep[v]->p.grp_data.levene = lz_numerator / lz_denominator[v] ;
368
369     }
370
371   /* Now clear up after ourselves */
372   free(lz_denominator);
373   free(lz);
374 }
375