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 "error.h"
  24 #include "casefile.h"
  25 #include "levene.h"
  26 #include "hash.h"
  27 #include "str.h"
  28 #include "var.h"
  29 #include "vfm.h"
  30 #include "alloc.h"
  31 #include "misc.h"
  32
  33
  34 #include <math.h>
  35 #include <stdlib.h>
  36
  37
  38 /* This module calculates the Levene statistic for variables.
  39
  40    Just for reference, the Levene Statistic is a defines as follows:
  41
  42    W = \frac{ (n-k)\sum_{i=1}^k n_i(Z_{iL} - Z_{LL})^2}
  43             { (k-1)\sum_{i=1}^k \sum_{j=1}^{n_i} (Z_{ij} - Z_{iL})^2}
  44
  45    where:
  46         k is the number of groups
  47         n is the total number of samples
  48         n_i is the number of samples in the ith group
  49         Z_{ij} is | Y_{ij} - Y_{iL} | where Y_{iL} is the mean of the ith group
  50         Z_{iL} is the  mean of Z_{ij} over the ith group
  51         Z_{LL} is the grand mean of Z_{ij}
  52
  53    Imagine calculating that with pencil and paper!
  54
  55  */
  56
  57 static struct group_statistics *get_group(int v, struct group_statistics *key);
  58
  59
  60 struct levene_info
  61 {
  62
  63   /* Per group statistics */
  64   struct t_test_proc **group_stats;
  65
  66   /* The independent variable */
  67   struct variable *v_indep;
  68
  69   /* Number of dependent variables */
  70   int n_dep;
  71
  72   /* The dependent variables */
  73   struct variable  **v_dep;
  74
  75   /* How to treat missing values */
  76   enum lev_missing missing;
  77
  78   /* Function to test for missing values */
  79   is_missing_func is_missing;
  80
  81 };
  82
  83 /* First pass */
  84 static void  levene_precalc (const struct levene_info *l);
  85 static int levene_calc (const struct ccase *, void *);
  86 static void levene_postcalc (void *);
  87
  88
  89 /* Second pass */
  90 static void levene2_precalc (void *);
  91 static int levene2_calc (const struct ccase *, void *);
  92 static void levene2_postcalc (void *);
  93
  94
  95 void
  96 levene(const struct casefile *cf,
  97        struct variable *v_indep, int n_dep, struct variable **v_dep,
  98              enum lev_missing missing,   is_missing_func value_is_missing)
  99 {
 100   struct casereader *r;
 101   const struct ccase *c;
 102   struct levene_info l;
 103
 104   l.n_dep      = n_dep;
 105   l.v_indep    = v_indep;
 106   l.v_dep      = v_dep;
 107   l.missing    = missing;
 108   l.is_missing = value_is_missing;
 109
 110
 111
 112   levene_precalc(&l);
 113   for(r = casefile_get_reader (cf);
 114       casereader_read (r, &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     {
 125       levene2_calc(c,&l);
 126     }
 127   casereader_destroy (r);
 128   levene2_postcalc(&l);
 129
 130 }
 131
 132 static struct hsh_table **hash;
 133
 134
 135 /* Return -1 if the id of a is less than b; +1 if greater than and
 136    0 if equal */
 137 static int
 138 compare_group(const struct group_statistics *a,
 139                  const struct group_statistics *b,
 140                  int width)
 141 {
 142   int id_cmp = compare_values(&a->id, &b->id, width);
 143
 144   if (id_cmp == 0 )
 145     {
 146       int c;
 147       c= memcmp(&a->criterion,&b->criterion,sizeof(enum comparison));
 148       return c;
 149     }
 150   else
 151     return id_cmp;
 152 }
 153
 154
 155 static unsigned
 156 hash_group(const struct group_statistics *g, int width)
 157 {
 158   unsigned id_hash;
 159
 160   if ( 0 == width )
 161     id_hash = hsh_hash_double (g->id.f);
 162   else
 163     id_hash = hsh_hash_bytes (g->id.s, width);
 164
 165   return id_hash;
 166 }
 167
 168 /* Internal variables used in calculating the Levene statistic */
 169
 170 /* Per variable statistics */
 171 struct lz_stats
 172 {
 173   /* Total of all lz */
 174   double grand_total;
 175
 176   /* Mean of all lz */
 177   double grand_mean;
 178
 179   /* The total number of cases */
 180   double total_n ;
 181
 182   /* Number of groups */
 183   int n_groups;
 184 };
 185
 186 /* An array of lz_stats for each variable */
 187 static struct lz_stats *lz;
 188
 189 /* Set to 1 if the groups require inequality comparisions */
 190 static int inequality_compare;
 191
 192
 193 static void
 194 levene_precalc (const struct levene_info *l)
 195 {
 196   int i;
 197
 198   lz  = xmalloc (sizeof (struct lz_stats ) * l->n_dep ) ;
 199
 200   hash = xmalloc (sizeof ( struct hsh_table *) * l->n_dep );
 201
 202   for(i=0; i < l->n_dep ; ++i )
 203     {
 204       struct variable *v = l->v_dep[i];
 205       int g;
 206       int number_of_groups = v->p.t_t.n_groups ;
 207
 208       hash[i] = hsh_create (l->n_dep * number_of_groups,
 209                             (hsh_compare_func *) compare_group,
 210                             (hsh_hash_func *) hash_group,
 211                             0,(void *) l->v_indep->width);
 212
 213       lz[i].grand_total = 0;
 214       lz[i].total_n = 0;
 215       lz[i].n_groups = number_of_groups;
 216
 217       for (g = 0 ; g < v->p.t_t.n_groups ; ++g )
 218         {
 219           struct group_statistics *gs = &v->p.t_t.gs[g];
 220           gs->lz_total = 0;
 221           hsh_insert(hash[i], gs);
 222           if ( gs->criterion != CMP_EQ )
 223             {
 224               inequality_compare = 1;
 225             }
 226         }
 227     }
 228
 229 }
 230
 231 static int
 232 levene_calc (const struct ccase *c, void *_l)
 233 {
 234   int i;
 235   int warn = 0;
 236   struct levene_info *l = (struct levene_info *) _l;
 237   const union value *gv = &c->data[l->v_indep->fv];
 238   struct group_statistics key;
 239   double weight = dict_get_case_weight(default_dict,c,&warn);
 240
 241
 242   /* Skip the entire case if /MISSING=LISTWISE is set */
 243   if ( l->missing == LEV_LISTWISE )
 244     {
 245       for (i = 0; i < l->n_dep; ++i)
 246         {
 247           struct variable *v = l->v_dep[i];
 248           const union value *val = &c->data[v->fv];
 249
 250           if (l->is_missing(val,v) )
 251             {
 252               return 0;
 253             }
 254         }
 255     }
 256
 257
 258   key.id = *gv;
 259   key.criterion = CMP_EQ;
 260
 261   for (i = 0; i < l->n_dep; ++i)
 262     {
 263       struct variable *var = l->v_dep[i];
 264       double levene_z;
 265       const union value *v = &c->data[var->fv];
 266       struct group_statistics *gs;
 267       gs = get_group(i,&key);
 268       if ( 0 == gs )
 269         continue ;
 270
 271       if ( ! l->is_missing(v,var))
 272         {
 273           levene_z= fabs(v->f - gs->mean);
 274           lz[i].grand_total += levene_z * weight;
 275           lz[i].total_n += weight;
 276
 277           gs->lz_total += levene_z * weight;
 278         }
 279     }
 280   return 0;
 281 }
 282
 283
 284 static void
 285 levene_postcalc (void *_l)
 286 {
 287   int v;
 288
 289   struct levene_info *l = (struct levene_info *) _l;
 290
 291   for (v = 0; v < l->n_dep; ++v)
 292     {
 293       lz[v].grand_mean = lz[v].grand_total / lz[v].total_n ;
 294
 295     }
 296
 297 }
 298
 299
 300 /* The denominator for the expression for the Levene */
 301 static double *lz_denominator;
 302
 303 static void
 304 levene2_precalc (void *_l)
 305 {
 306   int v;
 307
 308   struct levene_info *l = (struct levene_info *) _l;
 309
 310   lz_denominator = (double *) xmalloc(sizeof(double) * l->n_dep);
 311
 312   /* This stuff could go in the first post calc . . . */
 313   for (v = 0; v < l->n_dep; ++v)
 314     {
 315       struct hsh_iterator hi;
 316       struct group_statistics *g;
 317       for(g = (struct group_statistics *) hsh_first(hash[v],&hi);
 318           g != 0 ;
 319           g = (struct group_statistics *) hsh_next(hash[v],&hi) )
 320         {
 321           g->lz_mean = g->lz_total/g->n ;
 322         }
 323       lz_denominator[v] = 0;
 324   }
 325 }
 326
 327 static int
 328 levene2_calc (const struct ccase *c, void *_l)
 329 {
 330   int i;
 331   int warn = 0;
 332
 333   struct levene_info *l = (struct levene_info *) _l;
 334
 335   double weight = dict_get_case_weight(default_dict,c,&warn);
 336
 337   const union value *gv = &c->data[l->v_indep->fv];
 338   struct group_statistics key;
 339
 340   /* Skip the entire case if /MISSING=LISTWISE is set */
 341   if ( l->missing == LEV_LISTWISE )
 342     {
 343       for (i = 0; i < l->n_dep; ++i)
 344         {
 345           struct variable *v = l->v_dep[i];
 346           const union value *val = &c->data[v->fv];
 347
 348           if (l->is_missing(val,v) )
 349             {
 350               return 0;
 351             }
 352         }
 353     }
 354
 355   key.id = *gv;
 356   key.criterion = CMP_EQ;
 357
 358   for (i = 0; i < l->n_dep; ++i)
 359     {
 360       double levene_z;
 361       struct variable *var = l->v_dep[i] ;
 362       const union value *v = &c->data[var->fv];
 363       struct group_statistics *gs;
 364       gs = get_group(i,&key);
 365       if ( 0 == gs )
 366         continue;
 367
 368       if ( ! l->is_missing(v,var) )
 369         {
 370           levene_z = fabs(v->f - gs->mean);
 371           lz_denominator[i] += weight * pow2(levene_z - gs->lz_mean);
 372         }
 373     }
 374
 375   return 0;
 376 }
 377
 378
 379 static void
 380 levene2_postcalc (void *_l)
 381 {
 382   int v;
 383
 384   struct levene_info *l = (struct levene_info *) _l;
 385
 386   for (v = 0; v < l->n_dep; ++v)
 387     {
 388       double lz_numerator = 0;
 389       struct hsh_iterator hi;
 390       struct group_statistics *g;
 391       for(g = (struct group_statistics *) hsh_first(hash[v],&hi);
 392           g != 0 ;
 393           g = (struct group_statistics *) hsh_next(hash[v],&hi) )
 394         {
 395
 396           lz_numerator += g->n * pow2(g->lz_mean - lz[v].grand_mean );
 397
 398
 399         }
 400       lz_numerator *= ( l->v_dep[v]->p.t_t.ugs.n -
 401                         l->v_dep[v]->p.t_t.n_groups );
 402
 403       lz_denominator[v] /= (l->v_dep[v]->p.t_t.n_groups - 1);
 404
 405       l->v_dep[v]->p.t_t.levene = lz_numerator/lz_denominator[v] ;
 406     }
 407
 408   /* Now clear up after ourselves */
 409   free(lz_denominator);
 410   for (v = 0; v < l->n_dep; ++v)
 411     {
 412       hsh_destroy(hash[v]);
 413     }
 414
 415   free(hash);
 416   free(lz);
 417 }
 418
 419
 420 /* Return the group belonging to the v_th dependent variable
 421    which matches the key */
 422 static struct group_statistics *
 423 get_group(int v, struct group_statistics *key)
 424 {
 425   struct group_statistics *gs;
 426   gs = hsh_find(hash[v],key);
 427
 428
 429   if ( ( !gs )  && inequality_compare)
 430     {
 431       /* Here we degrade to a linear search.
 432          This would seem inefficient.  However, it should only ever happen
 433          with the T-TEST, for which there are exactly two groups */
 434
 435       struct hsh_iterator hi;
 436
 437       assert( hsh_count(hash[v]) == 2 ) ;
 438       for(gs = (struct group_statistics *) hsh_first(hash[v],&hi);
 439           gs != 0 ;
 440           gs = (struct group_statistics *) hsh_next(hash[v],&hi) )
 441         {
 442           int cmp;
 443
 444           cmp = compare_values(&gs->id, &key->id, 0);
 445
 446           assert( cmp != 0 ); /* or else the hash would have found something */
 447
 448           if ( cmp == -1 &&
 449                ( gs->criterion == CMP_GT || gs->criterion == CMP_GE )
 450              )
 451             break;
 452
 453           if ( cmp == 1 &&
 454                ( gs->criterion == CMP_LT || gs->criterion == CMP_LE )
 455              )
 456             break;
 457         }
 458     }
 459
 460   return gs;
 461 }