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