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
  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 static struct group_statistics *get_group(int v, struct group_statistics *key);
  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 static struct hsh_table **hash;
 136
 137
 138 /* Return -1 if the id of a is less than b; +1 if greater than and
 139    0 if equal */
 140 static int
 141 compare_group(const struct group_statistics *a,
 142                  const struct group_statistics *b,
 143                  int width)
 144 {
 145   int id_cmp = compare_values(&a->id, &b->id, width);
 146
 147   if (id_cmp == 0 )
 148     {
 149       int c;
 150       c= memcmp(&a->criterion,&b->criterion,sizeof(enum comparison));
 151       return c;
 152     }
 153   else
 154     return id_cmp;
 155 }
 156
 157
 158 static unsigned
 159 hash_group(const struct group_statistics *g, int width)
 160 {
 161   unsigned id_hash;
 162
 163   if ( 0 == width )
 164     id_hash = hsh_hash_double (g->id.f);
 165   else
 166     id_hash = hsh_hash_bytes (g->id.s, width);
 167
 168   return id_hash;
 169 }
 170
 171 /* Internal variables used in calculating the Levene statistic */
 172
 173 /* Per variable statistics */
 174 struct lz_stats
 175 {
 176   /* Total of all lz */
 177   double grand_total;
 178
 179   /* Mean of all lz */
 180   double grand_mean;
 181
 182   /* The total number of cases */
 183   double total_n ;
 184
 185   /* Number of groups */
 186   int n_groups;
 187 };
 188
 189 /* An array of lz_stats for each variable */
 190 static struct lz_stats *lz;
 191
 192 /* Set to 1 if the groups require inequality comparisions */
 193 static int inequality_compare;
 194
 195
 196 static void
 197 levene_precalc (const struct levene_info *l)
 198 {
 199   int i;
 200
 201   lz  = xmalloc (sizeof (struct lz_stats ) * l->n_dep ) ;
 202
 203   hash = xmalloc (sizeof ( struct hsh_table *) * l->n_dep );
 204
 205   for(i=0; i < l->n_dep ; ++i )
 206     {
 207       struct variable *v = l->v_dep[i];
 208       int g;
 209       int number_of_groups = v->p.t_t.n_groups ;
 210
 211       hash[i] = hsh_create (l->n_dep * number_of_groups,
 212                             (hsh_compare_func *) compare_group,
 213                             (hsh_hash_func *) hash_group,
 214                             0,(void *) l->v_indep->width);
 215
 216       lz[i].grand_total = 0;
 217       lz[i].total_n = 0;
 218       lz[i].n_groups = number_of_groups;
 219
 220       for (g = 0 ; g < v->p.t_t.n_groups ; ++g )
 221         {
 222           struct group_statistics *gs = &v->p.t_t.gs[g];
 223           gs->lz_total = 0;
 224           hsh_insert(hash[i], gs);
 225           if ( gs->criterion != CMP_EQ )
 226             {
 227               inequality_compare = 1;
 228             }
 229         }
 230     }
 231
 232 }
 233
 234 static int
 235 levene_calc (const struct ccase *c, void *_l)
 236 {
 237   int i;
 238   int warn = 0;
 239   struct levene_info *l = (struct levene_info *) _l;
 240   const union value *gv = case_data (c, l->v_indep->fv);
 241   struct group_statistics key;
 242   double weight = dict_get_case_weight(default_dict,c,&warn);
 243
 244
 245   /* Skip the entire case if /MISSING=LISTWISE is set */
 246   if ( l->missing == LEV_LISTWISE )
 247     {
 248       for (i = 0; i < l->n_dep; ++i)
 249         {
 250           struct variable *v = l->v_dep[i];
 251           const union value *val = case_data (c, v->fv);
 252
 253           if (l->is_missing(val,v) )
 254             {
 255               return 0;
 256             }
 257         }
 258     }
 259
 260
 261   key.id = *gv;
 262   key.criterion = CMP_EQ;
 263
 264   for (i = 0; i < l->n_dep; ++i)
 265     {
 266       struct variable *var = l->v_dep[i];
 267       double levene_z;
 268       const union value *v = case_data (c, var->fv);
 269       struct group_statistics *gs;
 270       gs = get_group(i,&key);
 271       if ( 0 == gs )
 272         continue ;
 273
 274       if ( ! l->is_missing(v,var))
 275         {
 276           levene_z= fabs(v->f - gs->mean);
 277           lz[i].grand_total += levene_z * weight;
 278           lz[i].total_n += weight;
 279
 280           gs->lz_total += levene_z * weight;
 281         }
 282     }
 283   return 0;
 284 }
 285
 286
 287 static void
 288 levene_postcalc (void *_l)
 289 {
 290   int v;
 291
 292   struct levene_info *l = (struct levene_info *) _l;
 293
 294   for (v = 0; v < l->n_dep; ++v)
 295     {
 296       lz[v].grand_mean = lz[v].grand_total / lz[v].total_n ;
 297
 298     }
 299
 300 }
 301
 302
 303 /* The denominator for the expression for the Levene */
 304 static double *lz_denominator;
 305
 306 static void
 307 levene2_precalc (void *_l)
 308 {
 309   int v;
 310
 311   struct levene_info *l = (struct levene_info *) _l;
 312
 313   lz_denominator = (double *) xmalloc(sizeof(double) * l->n_dep);
 314
 315   /* This stuff could go in the first post calc . . . */
 316   for (v = 0; v < l->n_dep; ++v)
 317     {
 318       struct hsh_iterator hi;
 319       struct group_statistics *g;
 320       for(g = (struct group_statistics *) hsh_first(hash[v],&hi);
 321           g != 0 ;
 322           g = (struct group_statistics *) hsh_next(hash[v],&hi) )
 323         {
 324           g->lz_mean = g->lz_total/g->n ;
 325         }
 326       lz_denominator[v] = 0;
 327   }
 328 }
 329
 330 static int
 331 levene2_calc (const struct ccase *c, void *_l)
 332 {
 333   int i;
 334   int warn = 0;
 335
 336   struct levene_info *l = (struct levene_info *) _l;
 337
 338   double weight = dict_get_case_weight(default_dict,c,&warn);
 339
 340   const union value *gv = case_data (c, l->v_indep->fv);
 341   struct group_statistics key;
 342
 343   /* Skip the entire case if /MISSING=LISTWISE is set */
 344   if ( l->missing == LEV_LISTWISE )
 345     {
 346       for (i = 0; i < l->n_dep; ++i)
 347         {
 348           struct variable *v = l->v_dep[i];
 349           const union value *val = case_data (c, v->fv);
 350
 351           if (l->is_missing(val,v) )
 352             {
 353               return 0;
 354             }
 355         }
 356     }
 357
 358   key.id = *gv;
 359   key.criterion = CMP_EQ;
 360
 361   for (i = 0; i < l->n_dep; ++i)
 362     {
 363       double levene_z;
 364       struct variable *var = l->v_dep[i] ;
 365       const union value *v = case_data (c, var->fv);
 366       struct group_statistics *gs;
 367       gs = get_group(i,&key);
 368       if ( 0 == gs )
 369         continue;
 370
 371       if ( ! l->is_missing(v,var) )
 372         {
 373           levene_z = fabs(v->f - gs->mean);
 374           lz_denominator[i] += weight * pow2(levene_z - gs->lz_mean);
 375         }
 376     }
 377
 378   return 0;
 379 }
 380
 381
 382 static void
 383 levene2_postcalc (void *_l)
 384 {
 385   int v;
 386
 387   struct levene_info *l = (struct levene_info *) _l;
 388
 389   for (v = 0; v < l->n_dep; ++v)
 390     {
 391       double lz_numerator = 0;
 392       struct hsh_iterator hi;
 393       struct group_statistics *g;
 394       for(g = (struct group_statistics *) hsh_first(hash[v],&hi);
 395           g != 0 ;
 396           g = (struct group_statistics *) hsh_next(hash[v],&hi) )
 397         {
 398
 399           lz_numerator += g->n * pow2(g->lz_mean - lz[v].grand_mean );
 400
 401
 402         }
 403       lz_numerator *= ( l->v_dep[v]->p.t_t.ugs.n -
 404                         l->v_dep[v]->p.t_t.n_groups );
 405
 406       lz_denominator[v] /= (l->v_dep[v]->p.t_t.n_groups - 1);
 407
 408       l->v_dep[v]->p.t_t.levene = lz_numerator/lz_denominator[v] ;
 409     }
 410
 411   /* Now clear up after ourselves */
 412   free(lz_denominator);
 413   for (v = 0; v < l->n_dep; ++v)
 414     {
 415       hsh_destroy(hash[v]);
 416     }
 417
 418   free(hash);
 419   free(lz);
 420 }
 421
 422
 423 /* Return the group belonging to the v_th dependent variable
 424    which matches the key */
 425 static struct group_statistics *
 426 get_group(int v, struct group_statistics *key)
 427 {
 428   struct group_statistics *gs;
 429   gs = hsh_find(hash[v],key);
 430
 431
 432   if ( ( !gs )  && inequality_compare)
 433     {
 434       /* Here we degrade to a linear search.
 435          This would seem inefficient.  However, it should only ever happen
 436          with the T-TEST, for which there are exactly two groups */
 437
 438       struct hsh_iterator hi;
 439
 440       assert( hsh_count(hash[v]) == 2 ) ;
 441       for(gs = (struct group_statistics *) hsh_first(hash[v],&hi);
 442           gs != 0 ;
 443           gs = (struct group_statistics *) hsh_next(hash[v],&hi) )
 444         {
 445           int cmp;
 446
 447           cmp = compare_values(&gs->id, &key->id, 0);
 448
 449           assert( cmp != 0 ); /* or else the hash would have found something */
 450
 451           if ( cmp == -1 &&
 452                ( gs->criterion == CMP_GT || gs->criterion == CMP_GE )
 453              )
 454             break;
 455
 456           if ( cmp == 1 &&
 457                ( gs->criterion == CMP_LT || gs->criterion == CMP_LE )
 458              )
 459             break;
 460         }
 461     }
 462
 463   return gs;
 464 }