src/math/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., 51 Franklin Street, Fifth Floor, Boston, MA
  20    02110-1301, USA. */
  21
  22 #include <config.h>
  23 #include "levene.h"
  24 #include <libpspp/message.h>
  25 #include <data/case.h>
  26 #include <data/casefile.h>
  27 #include <data/dictionary.h>
  28 #include "group-proc.h"
  29 #include <libpspp/hash.h>
  30 #include <libpspp/str.h>
  31 #include <data/variable.h>
  32 #include <data/procedure.h>
  33 #include <libpspp/alloc.h>
  34 #include <libpspp/misc.h>
  35 #include "group.h"
  36
  37 #include <math.h>
  38 #include <stdlib.h>
  39
  40
  41 /* This module calculates the Levene statistic for variables.
  42
  43    Just for reference, the Levene Statistic is a defines as follows:
  44
  45    W = \frac{ (n-k)\sum_{i=1}^k n_i(Z_{iL} - Z_{LL})^2}
  46             { (k-1)\sum_{i=1}^k \sum_{j=1}^{n_i} (Z_{ij} - Z_{iL})^2}
  47
  48    where:
  49         k is the number of groups
  50         n is the total number of samples
  51         n_i is the number of samples in the ith group
  52         Z_{ij} is | Y_{ij} - Y_{iL} | where Y_{iL} is the mean of the ith group
  53         Z_{iL} is the  mean of Z_{ij} over the ith group
  54         Z_{LL} is the grand mean of Z_{ij}
  55
  56    Imagine calculating that with pencil and paper!
  57
  58  */
  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   size_t 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 /* First pass */
  84 static void  levene_precalc (const struct levene_info *l);
  85 static int levene_calc (const struct ccase *, void *);
  86 static void levene_postcalc (void *);
  87
  88
  89 /* Second pass */
  90 static void levene2_precalc (void *);
  91 static int levene2_calc (const struct ccase *, void *);
  92 static void levene2_postcalc (void *);
  93
  94
  95 void
  96 levene(const struct casefile *cf,
  97        struct variable *v_indep, size_t n_dep, struct variable **v_dep,
  98              enum lev_missing missing,   is_missing_func value_is_missing)
  99 {
 100   struct casereader *r;
 101   struct ccase c;
 102   struct levene_info l;
 103
 104   l.n_dep      = n_dep;
 105   l.v_indep    = v_indep;
 106   l.v_dep      = v_dep;
 107   l.missing    = missing;
 108   l.is_missing = value_is_missing;
 109
 110
 111
 112   levene_precalc(&l);
 113   for(r = casefile_get_reader (cf);
 114       casereader_read (r, &c) ;
 115       case_destroy (&c))
 116     {
 117       levene_calc(&c,&l);
 118     }
 119   casereader_destroy (r);
 120   levene_postcalc(&l);
 121
 122   levene2_precalc(&l);
 123   for(r = casefile_get_reader (cf);
 124       casereader_read (r, &c) ;
 125       case_destroy (&c))
 126     {
 127       levene2_calc(&c,&l);
 128     }
 129   casereader_destroy (r);
 130   levene2_postcalc(&l);
 131
 132 }
 133
 134 /* Internal variables used in calculating the Levene statistic */
 135
 136 /* Per variable statistics */
 137 struct lz_stats
 138 {
 139   /* Total of all lz */
 140   double grand_total;
 141
 142   /* Mean of all lz */
 143   double grand_mean;
 144
 145   /* The total number of cases */
 146   double total_n ;
 147
 148   /* Number of groups */
 149   int n_groups;
 150 };
 151
 152 /* An array of lz_stats for each variable */
 153 static struct lz_stats *lz;
 154
 155
 156 static void
 157 levene_precalc (const struct levene_info *l)
 158 {
 159   size_t i;
 160
 161   lz = xnmalloc (l->n_dep, sizeof *lz);
 162
 163   for(i = 0; i < l->n_dep ; ++i )
 164     {
 165       struct variable *var = l->v_dep[i];
 166       struct group_proc *gp = group_proc_get (var);
 167       struct group_statistics *gs;
 168       struct hsh_iterator hi;
 169
 170       lz[i].grand_total = 0;
 171       lz[i].total_n = 0;
 172       lz[i].n_groups = gp->n_groups ;
 173
 174
 175       for ( gs = hsh_first(gp->group_hash, &hi);
 176             gs != 0;
 177             gs = hsh_next(gp->group_hash, &hi))
 178         {
 179           gs->lz_total = 0;
 180         }
 181
 182     }
 183
 184 }
 185
 186 static int
 187 levene_calc (const struct ccase *c, void *_l)
 188 {
 189   size_t i;
 190   bool warn = false;
 191   struct levene_info *l = (struct levene_info *) _l;
 192   const union value *gv = case_data (c, l->v_indep->fv);
 193   struct group_statistics key;
 194   double weight = dict_get_case_weight (dataset_dict (current_dataset), c, &warn);
 195
 196   /* Skip the entire case if /MISSING=LISTWISE is set */
 197   if ( l->missing == LEV_LISTWISE )
 198     {
 199       for (i = 0; i < l->n_dep; ++i)
 200         {
 201           struct variable *v = l->v_dep[i];
 202           const union value *val = case_data (c, v->fv);
 203
 204           if (l->is_missing (&v->miss, val) )
 205             {
 206               return 0;
 207             }
 208         }
 209     }
 210
 211
 212   key.id = *gv;
 213
 214   for (i = 0; i < l->n_dep; ++i)
 215     {
 216       struct variable *var = l->v_dep[i];
 217       struct group_proc *gp = group_proc_get (var);
 218       double levene_z;
 219       const union value *v = case_data (c, var->fv);
 220       struct group_statistics *gs;
 221
 222       gs = hsh_find(gp->group_hash,(void *) &key );
 223
 224       if ( 0 == gs )
 225         continue ;
 226
 227       if ( ! l->is_missing(&var->miss, v))
 228         {
 229           levene_z= fabs(v->f - gs->mean);
 230           lz[i].grand_total += levene_z * weight;
 231           lz[i].total_n += weight;
 232
 233           gs->lz_total += levene_z * weight;
 234         }
 235
 236     }
 237   return 0;
 238 }
 239
 240
 241 static void
 242 levene_postcalc (void *_l)
 243 {
 244   size_t v;
 245
 246   struct levene_info *l = (struct levene_info *) _l;
 247
 248   for (v = 0; v < l->n_dep; ++v)
 249     {
 250       /* This is Z_LL */
 251       lz[v].grand_mean = lz[v].grand_total / lz[v].total_n ;
 252     }
 253
 254
 255 }
 256
 257
 258 /* The denominator for the expression for the Levene */
 259 static double *lz_denominator;
 260
 261 static void
 262 levene2_precalc (void *_l)
 263 {
 264   size_t v;
 265
 266   struct levene_info *l = (struct levene_info *) _l;
 267
 268   lz_denominator = xnmalloc (l->n_dep, sizeof *lz_denominator);
 269
 270   /* This stuff could go in the first post calc . . . */
 271   for (v = 0; v < l->n_dep; ++v)
 272     {
 273       struct hsh_iterator hi;
 274       struct group_statistics *g;
 275
 276       struct variable *var = l->v_dep[v] ;
 277       struct hsh_table *hash = group_proc_get (var)->group_hash;
 278
 279
 280       for(g = (struct group_statistics *) hsh_first(hash,&hi);
 281           g != 0 ;
 282           g = (struct group_statistics *) hsh_next(hash,&hi) )
 283         {
 284           g->lz_mean = g->lz_total / g->n ;
 285         }
 286       lz_denominator[v] = 0;
 287   }
 288 }
 289
 290 static int
 291 levene2_calc (const struct ccase *c, void *_l)
 292 {
 293   size_t i;
 294   bool warn = false;
 295
 296   struct levene_info *l = (struct levene_info *) _l;
 297
 298   double weight = dict_get_case_weight (dataset_dict (current_dataset), c, &warn);
 299
 300   const union value *gv = case_data (c, l->v_indep->fv);
 301   struct group_statistics key;
 302
 303   /* Skip the entire case if /MISSING=LISTWISE is set */
 304   if ( l->missing == LEV_LISTWISE )
 305     {
 306       for (i = 0; i < l->n_dep; ++i)
 307         {
 308           struct variable *v = l->v_dep[i];
 309           const union value *val = case_data (c, v->fv);
 310
 311           if (l->is_missing(&v->miss, val) )
 312             {
 313               return 0;
 314             }
 315         }
 316     }
 317
 318   key.id = *gv;
 319
 320   for (i = 0; i < l->n_dep; ++i)
 321     {
 322       double levene_z;
 323       struct variable *var = l->v_dep[i] ;
 324       const union value *v = case_data (c, var->fv);
 325       struct group_statistics *gs;
 326
 327       gs = hsh_find(group_proc_get (var)->group_hash,(void *) &key );
 328
 329       if ( 0 == gs )
 330         continue;
 331
 332       if ( ! l->is_missing (&var->miss, v) )
 333         {
 334           levene_z = fabs(v->f - gs->mean);
 335           lz_denominator[i] += weight * pow2(levene_z - gs->lz_mean);
 336         }
 337     }
 338
 339   return 0;
 340 }
 341
 342
 343 static void
 344 levene2_postcalc (void *_l)
 345 {
 346   size_t v;
 347
 348   struct levene_info *l = (struct levene_info *) _l;
 349
 350   for (v = 0; v < l->n_dep; ++v)
 351     {
 352       double lz_numerator = 0;
 353       struct hsh_iterator hi;
 354       struct group_statistics *g;
 355
 356       struct variable *var = l->v_dep[v] ;
 357       struct group_proc *gp = group_proc_get (var);
 358       struct hsh_table *hash = gp->group_hash;
 359
 360       for(g = (struct group_statistics *) hsh_first(hash,&hi);
 361           g != 0 ;
 362           g = (struct group_statistics *) hsh_next(hash,&hi) )
 363         {
 364           lz_numerator += g->n * pow2(g->lz_mean - lz[v].grand_mean );
 365         }
 366       lz_numerator *= ( gp->ugs.n - gp->n_groups );
 367
 368       lz_denominator[v] *= (gp->n_groups - 1);
 369
 370       gp->levene = lz_numerator / lz_denominator[v] ;
 371
 372     }
 373
 374   /* Now clear up after ourselves */
 375   free(lz_denominator);
 376   free(lz);
 377 }
 378