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
   6    This program is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU General Public License as
   8    published by the Free Software Foundation; either version 2 of the
   9    License, or (at your option) any later version.
  10
  11    This program is distributed in the hope that it will be useful, but
  12    WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    General Public License for more details.
  15
  16    You should have received a copy of the GNU General Public License
  17    along with this program; if not, write to the Free Software
  18    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  19    02110-1301, USA. */
  20
  21 #include <config.h>
  22 #include "levene.h"
  23 #include <libpspp/message.h>
  24 #include <data/case.h>
  25 #include <data/casefile.h>
  26 #include <data/dictionary.h>
  27 #include "group-proc.h"
  28 #include <libpspp/hash.h>
  29 #include <libpspp/str.h>
  30 #include <data/variable.h>
  31 #include <data/procedure.h>
  32 #include <data/casefilter.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   /* Filter for missing values */
  77   struct casefilter *filter;
  78 };
  79
  80 /* First pass */
  81 static void  levene_precalc (const struct levene_info *l);
  82 static int levene_calc (const struct dictionary *dict, const struct ccase *,
  83                         const struct levene_info *l);
  84 static void levene_postcalc (void *);
  85
  86
  87 /* Second pass */
  88 static void levene2_precalc (struct levene_info *l);
  89 static int levene2_calc (const struct dictionary *, const struct ccase *,
  90                          struct levene_info *l);
  91 static void levene2_postcalc (void *);
  92
  93
  94 void
  95 levene(const struct dictionary *dict,
  96        const struct casefile *cf,
  97        struct variable *v_indep, size_t n_dep, struct variable **v_dep,
  98        struct casefilter *filter)
  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.filter = filter;
 108
 109
 110   levene_precalc (&l);
 111   for(r = casefile_get_reader (cf, filter);
 112       casereader_read (r, &c) ;
 113       case_destroy (&c))
 114     {
 115       levene_calc (dict, &c, &l);
 116     }
 117   casereader_destroy (r);
 118   levene_postcalc (&l);
 119
 120   levene2_precalc(&l);
 121   for(r = casefile_get_reader (cf, filter);
 122       casereader_read (r, &c) ;
 123       case_destroy (&c))
 124     {
 125       levene2_calc (dict, &c,&l);
 126     }
 127   casereader_destroy (r);
 128   levene2_postcalc (&l);
 129 }
 130
 131 /* Internal variables used in calculating the Levene statistic */
 132
 133 /* Per variable statistics */
 134 struct lz_stats
 135 {
 136   /* Total of all lz */
 137   double grand_total;
 138
 139   /* Mean of all lz */
 140   double grand_mean;
 141
 142   /* The total number of cases */
 143   double total_n ;
 144
 145   /* Number of groups */
 146   int n_groups;
 147 };
 148
 149 /* An array of lz_stats for each variable */
 150 static struct lz_stats *lz;
 151
 152
 153 static void
 154 levene_precalc (const struct levene_info *l)
 155 {
 156   size_t i;
 157
 158   lz = xnmalloc (l->n_dep, sizeof *lz);
 159
 160   for(i = 0; i < l->n_dep ; ++i )
 161     {
 162       struct variable *var = l->v_dep[i];
 163       struct group_proc *gp = group_proc_get (var);
 164       struct group_statistics *gs;
 165       struct hsh_iterator hi;
 166
 167       lz[i].grand_total = 0;
 168       lz[i].total_n = 0;
 169       lz[i].n_groups = gp->n_groups ;
 170
 171
 172       for ( gs = hsh_first(gp->group_hash, &hi);
 173             gs != 0;
 174             gs = hsh_next(gp->group_hash, &hi))
 175         {
 176           gs->lz_total = 0;
 177         }
 178
 179     }
 180
 181 }
 182
 183 static int
 184 levene_calc (const struct dictionary *dict, const struct ccase *c,
 185              const struct levene_info *l)
 186 {
 187   size_t i;
 188   bool warn = false;
 189   const union value *gv = case_data (c, l->v_indep);
 190   struct group_statistics key;
 191   double weight = dict_get_case_weight (dict, c, &warn);
 192
 193   key.id = *gv;
 194
 195   for (i = 0; i < l->n_dep; ++i)
 196     {
 197       struct variable *var = l->v_dep[i];
 198       struct group_proc *gp = group_proc_get (var);
 199       double levene_z;
 200       const union value *v = case_data (c, var);
 201       struct group_statistics *gs;
 202
 203       gs = hsh_find(gp->group_hash,(void *) &key );
 204
 205       if ( 0 == gs )
 206         continue ;
 207
 208       if ( ! casefilter_variable_missing (l->filter, c, var))
 209         {
 210           levene_z= fabs(v->f - gs->mean);
 211           lz[i].grand_total += levene_z * weight;
 212           lz[i].total_n += weight;
 213
 214           gs->lz_total += levene_z * weight;
 215         }
 216     }
 217   return 0;
 218 }
 219
 220
 221 static void
 222 levene_postcalc (void *_l)
 223 {
 224   size_t v;
 225
 226   struct levene_info *l = (struct levene_info *) _l;
 227
 228   for (v = 0; v < l->n_dep; ++v)
 229     {
 230       /* This is Z_LL */
 231       lz[v].grand_mean = lz[v].grand_total / lz[v].total_n ;
 232     }
 233
 234
 235 }
 236
 237
 238
 239 /* The denominator for the expression for the Levene */
 240 static double *lz_denominator = 0;
 241
 242 static void
 243 levene2_precalc (struct levene_info *l)
 244 {
 245   size_t v;
 246
 247   lz_denominator = xnmalloc (l->n_dep, sizeof *lz_denominator);
 248
 249   /* This stuff could go in the first post calc . . . */
 250   for (v = 0;
 251        v < l->n_dep;
 252        ++v)
 253     {
 254       struct hsh_iterator hi;
 255       struct group_statistics *g;
 256
 257       struct variable *var = l->v_dep[v] ;
 258       struct hsh_table *hash = group_proc_get (var)->group_hash;
 259
 260
 261       for(g = (struct group_statistics *) hsh_first(hash,&hi);
 262           g != 0 ;
 263           g = (struct group_statistics *) hsh_next(hash,&hi) )
 264         {
 265           g->lz_mean = g->lz_total / g->n ;
 266         }
 267       lz_denominator[v] = 0;
 268   }
 269 }
 270
 271 static int
 272 levene2_calc (const struct dictionary *dict, const struct ccase *c,
 273               struct levene_info *l)
 274 {
 275   size_t i;
 276   bool warn = false;
 277
 278   double weight = dict_get_case_weight (dict, c, &warn);
 279
 280   const union value *gv = case_data (c, l->v_indep);
 281   struct group_statistics key;
 282
 283   key.id = *gv;
 284
 285   for (i = 0; i < l->n_dep; ++i)
 286     {
 287       double levene_z;
 288       struct variable *var = l->v_dep[i] ;
 289       const union value *v = case_data (c, var);
 290       struct group_statistics *gs;
 291
 292       gs = hsh_find(group_proc_get (var)->group_hash,(void *) &key );
 293
 294       if ( 0 == gs )
 295         continue;
 296
 297       if ( ! casefilter_variable_missing (l->filter, c, var))
 298
 299         {
 300           levene_z = fabs(v->f - gs->mean);
 301           lz_denominator[i] += weight * pow2 (levene_z - gs->lz_mean);
 302         }
 303     }
 304
 305   return 0;
 306 }
 307
 308
 309 static void
 310 levene2_postcalc (void *_l)
 311 {
 312   size_t v;
 313
 314   struct levene_info *l = (struct levene_info *) _l;
 315
 316   for (v = 0; v < l->n_dep; ++v)
 317     {
 318       double lz_numerator = 0;
 319       struct hsh_iterator hi;
 320       struct group_statistics *g;
 321
 322       struct variable *var = l->v_dep[v] ;
 323       struct group_proc *gp = group_proc_get (var);
 324       struct hsh_table *hash = gp->group_hash;
 325
 326       for(g = (struct group_statistics *) hsh_first(hash,&hi);
 327           g != 0 ;
 328           g = (struct group_statistics *) hsh_next(hash,&hi) )
 329         {
 330           lz_numerator += g->n * pow2(g->lz_mean - lz[v].grand_mean );
 331         }
 332       lz_numerator *= ( gp->ugs.n - gp->n_groups );
 333
 334       lz_denominator[v] *= (gp->n_groups - 1);
 335
 336       gp->levene = lz_numerator / lz_denominator[v] ;
 337
 338     }
 339
 340   /* Now clear up after ourselves */
 341   free(lz_denominator);
 342   free(lz);
 343 }
 344