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