1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2004, 2009 Free Software Foundation, Inc.
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.
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.
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/>. */
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>
37 /* This module calculates the Levene statistic for variables.
39 Just for reference, the Levene Statistic is a defines as follows:
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}
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}
52 Imagine calculating that with pencil and paper!
60 /* Per group statistics */
61 struct t_test_proc **group_stats;
63 /* The independent variable */
64 const struct variable *v_indep;
66 /* Number of dependent variables */
69 /* The dependent variables */
70 const struct variable **v_dep;
72 /* Filter for missing values */
73 enum mv_class exclude;
75 /* An array of lz_stats for each variable */
78 /* The denominator for the expression for the Levene */
79 double *lz_denominator;
83 /* Per variable statistics */
92 /* The total number of cases */
95 /* Number of groups */
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 *);
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 *);
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)
120 struct casereader *pass1, *pass2;
122 struct levene_info l;
128 l.lz = xnmalloc (l.n_dep, sizeof *l.lz);
129 l.lz_denominator = xnmalloc (l.n_dep, sizeof *l.lz_denominator);
131 casereader_split (reader, &pass1, &pass2);
134 for (; (c = casereader_read (pass1)) != NULL; case_unref (c))
135 levene_calc (dict, c, &l);
136 casereader_destroy (pass1);
137 levene_postcalc (&l);
140 for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
141 levene2_calc (dict, c, &l);
142 casereader_destroy (pass2);
143 levene2_postcalc (&l);
145 free (l.lz_denominator);
150 levene_precalc (const struct levene_info *l)
154 for(i = 0; i < l->n_dep ; ++i )
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;
161 l->lz[i].grand_total = 0;
162 l->lz[i].total_n = 0;
163 l->lz[i].n_groups = gp->n_groups ;
166 for ( gs = hsh_first(gp->group_hash, &hi);
168 gs = hsh_next(gp->group_hash, &hi))
178 levene_calc (const struct dictionary *dict, const struct ccase *c,
179 const struct levene_info *l)
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);
189 for (i = 0; i < l->n_dep; ++i)
191 const struct variable *var = l->v_dep[i];
192 struct group_proc *gp = group_proc_get (var);
194 const union value *v = case_data (c, var);
195 struct group_statistics *gs;
197 gs = hsh_find(gp->group_hash,(void *) &key );
202 if ( !var_is_value_missing (var, v, l->exclude))
204 levene_z= fabs(v->f - gs->mean);
205 l->lz[i].grand_total += levene_z * weight;
206 l->lz[i].total_n += weight;
208 gs->lz_total += levene_z * weight;
216 levene_postcalc (struct levene_info *l)
220 for (v = 0; v < l->n_dep; ++v)
223 l->lz[v].grand_mean = l->lz[v].grand_total / l->lz[v].total_n ;
232 levene2_precalc (struct levene_info *l)
237 /* This stuff could go in the first post calc . . . */
242 struct hsh_iterator hi;
243 struct group_statistics *g;
245 const struct variable *var = l->v_dep[v] ;
246 struct hsh_table *hash = group_proc_get (var)->group_hash;
249 for(g = (struct group_statistics *) hsh_first(hash,&hi);
251 g = (struct group_statistics *) hsh_next(hash,&hi) )
253 g->lz_mean = g->lz_total / g->n ;
255 l->lz_denominator[v] = 0;
260 levene2_calc (const struct dictionary *dict, const struct ccase *c,
261 struct levene_info *l)
266 double weight = dict_get_case_weight (dict, c, &warn);
268 const union value *gv = case_data (c, l->v_indep);
269 struct group_statistics key;
273 for (i = 0; i < l->n_dep; ++i)
276 const struct variable *var = l->v_dep[i] ;
277 const union value *v = case_data (c, var);
278 struct group_statistics *gs;
280 gs = hsh_find(group_proc_get (var)->group_hash,(void *) &key );
285 if ( !var_is_value_missing (var, v, l->exclude))
287 levene_z = fabs(v->f - gs->mean);
288 l->lz_denominator[i] += weight * pow2 (levene_z - gs->lz_mean);
297 levene2_postcalc (struct levene_info *l)
301 for (v = 0; v < l->n_dep; ++v)
303 double lz_numerator = 0;
304 struct hsh_iterator hi;
305 struct group_statistics *g;
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;
311 for(g = (struct group_statistics *) hsh_first(hash,&hi);
313 g = (struct group_statistics *) hsh_next(hash,&hi) )
315 lz_numerator += g->n * pow2(g->lz_mean - l->lz[v].grand_mean );
317 lz_numerator *= ( gp->ugs.n - gp->n_groups );
319 l->lz_denominator[v] *= (gp->n_groups - 1);
321 gp->levene = lz_numerator / l->lz_denominator[v] ;