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