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