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