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