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