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