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