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   const struct variable *v_indep;
  69
  70   /* Number of dependent variables */
  71   size_t n_dep;
  72
  73   /* The dependent variables */
  74   const 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        const struct variable *v_indep, size_t n_dep,
  98        const struct variable **v_dep,
  99        struct casefilter *filter)
 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.filter = filter;
 109
 110
 111   levene_precalc (&l);
 112   for(r = casefile_get_reader (cf, filter);
 113       casereader_read (r, &c) ;
 114       case_destroy (&c))
 115     {
 116       levene_calc (dict, &c, &l);
 117     }
 118   casereader_destroy (r);
 119   levene_postcalc (&l);
 120
 121   levene2_precalc(&l);
 122   for(r = casefile_get_reader (cf, filter);
 123       casereader_read (r, &c) ;
 124       case_destroy (&c))
 125     {
 126       levene2_calc (dict, &c,&l);
 127     }
 128   casereader_destroy (r);
 129   levene2_postcalc (&l);
 130 }
 131
 132 /* Internal variables used in calculating the Levene statistic */
 133
 134 /* Per variable statistics */
 135 struct lz_stats
 136 {
 137   /* Total of all lz */
 138   double grand_total;
 139
 140   /* Mean of all lz */
 141   double grand_mean;
 142
 143   /* The total number of cases */
 144   double total_n ;
 145
 146   /* Number of groups */
 147   int n_groups;
 148 };
 149
 150 /* An array of lz_stats for each variable */
 151 static struct lz_stats *lz;
 152
 153
 154 static void
 155 levene_precalc (const struct levene_info *l)
 156 {
 157   size_t i;
 158
 159   lz = xnmalloc (l->n_dep, sizeof *lz);
 160
 161   for(i = 0; i < l->n_dep ; ++i )
 162     {
 163       const struct variable *var = l->v_dep[i];
 164       struct group_proc *gp = group_proc_get (var);
 165       struct group_statistics *gs;
 166       struct hsh_iterator hi;
 167
 168       lz[i].grand_total = 0;
 169       lz[i].total_n = 0;
 170       lz[i].n_groups = gp->n_groups ;
 171
 172
 173       for ( gs = hsh_first(gp->group_hash, &hi);
 174             gs != 0;
 175             gs = hsh_next(gp->group_hash, &hi))
 176         {
 177           gs->lz_total = 0;
 178         }
 179
 180     }
 181
 182 }
 183
 184 static int
 185 levene_calc (const struct dictionary *dict, const struct ccase *c,
 186              const struct levene_info *l)
 187 {
 188   size_t i;
 189   bool warn = false;
 190   const union value *gv = case_data (c, l->v_indep);
 191   struct group_statistics key;
 192   double weight = dict_get_case_weight (dict, c, &warn);
 193
 194   key.id = *gv;
 195
 196   for (i = 0; i < l->n_dep; ++i)
 197     {
 198       const struct variable *var = l->v_dep[i];
 199       struct group_proc *gp = group_proc_get (var);
 200       double levene_z;
 201       const union value *v = case_data (c, var);
 202       struct group_statistics *gs;
 203
 204       gs = hsh_find(gp->group_hash,(void *) &key );
 205
 206       if ( 0 == gs )
 207         continue ;
 208
 209       if ( ! casefilter_variable_missing (l->filter, c, var))
 210         {
 211           levene_z= fabs(v->f - gs->mean);
 212           lz[i].grand_total += levene_z * weight;
 213           lz[i].total_n += weight;
 214
 215           gs->lz_total += levene_z * weight;
 216         }
 217     }
 218   return 0;
 219 }
 220
 221
 222 static void
 223 levene_postcalc (void *_l)
 224 {
 225   size_t v;
 226
 227   struct levene_info *l = (struct levene_info *) _l;
 228
 229   for (v = 0; v < l->n_dep; ++v)
 230     {
 231       /* This is Z_LL */
 232       lz[v].grand_mean = lz[v].grand_total / lz[v].total_n ;
 233     }
 234
 235
 236 }
 237
 238
 239
 240 /* The denominator for the expression for the Levene */
 241 static double *lz_denominator = 0;
 242
 243 static void
 244 levene2_precalc (struct levene_info *l)
 245 {
 246   size_t v;
 247
 248   lz_denominator = xnmalloc (l->n_dep, sizeof *lz_denominator);
 249
 250   /* This stuff could go in the first post calc . . . */
 251   for (v = 0;
 252        v < l->n_dep;
 253        ++v)
 254     {
 255       struct hsh_iterator hi;
 256       struct group_statistics *g;
 257
 258       const struct variable *var = l->v_dep[v] ;
 259       struct hsh_table *hash = group_proc_get (var)->group_hash;
 260
 261
 262       for(g = (struct group_statistics *) hsh_first(hash,&hi);
 263           g != 0 ;
 264           g = (struct group_statistics *) hsh_next(hash,&hi) )
 265         {
 266           g->lz_mean = g->lz_total / g->n ;
 267         }
 268       lz_denominator[v] = 0;
 269   }
 270 }
 271
 272 static int
 273 levene2_calc (const struct dictionary *dict, const struct ccase *c,
 274               struct levene_info *l)
 275 {
 276   size_t i;
 277   bool warn = false;
 278
 279   double weight = dict_get_case_weight (dict, c, &warn);
 280
 281   const union value *gv = case_data (c, l->v_indep);
 282   struct group_statistics key;
 283
 284   key.id = *gv;
 285
 286   for (i = 0; i < l->n_dep; ++i)
 287     {
 288       double levene_z;
 289       const struct variable *var = l->v_dep[i] ;
 290       const union value *v = case_data (c, var);
 291       struct group_statistics *gs;
 292
 293       gs = hsh_find(group_proc_get (var)->group_hash,(void *) &key );
 294
 295       if ( 0 == gs )
 296         continue;
 297
 298       if ( ! casefilter_variable_missing (l->filter, c, var))
 299
 300         {
 301           levene_z = fabs(v->f - gs->mean);
 302           lz_denominator[i] += weight * pow2 (levene_z - gs->lz_mean);
 303         }
 304     }
 305
 306   return 0;
 307 }
 308
 309
 310 static void
 311 levene2_postcalc (void *_l)
 312 {
 313   size_t v;
 314
 315   struct levene_info *l = (struct levene_info *) _l;
 316
 317   for (v = 0; v < l->n_dep; ++v)
 318     {
 319       double lz_numerator = 0;
 320       struct hsh_iterator hi;
 321       struct group_statistics *g;
 322
 323       const struct variable *var = l->v_dep[v] ;
 324       struct group_proc *gp = group_proc_get (var);
 325       struct hsh_table *hash = gp->group_hash;
 326
 327       for(g = (struct group_statistics *) hsh_first(hash,&hi);
 328           g != 0 ;
 329           g = (struct group_statistics *) hsh_next(hash,&hi) )
 330         {
 331           lz_numerator += g->n * pow2(g->lz_mean - lz[v].grand_mean );
 332         }
 333       lz_numerator *= ( gp->ugs.n - gp->n_groups );
 334
 335       lz_denominator[v] *= (gp->n_groups - 1);
 336
 337       gp->levene = lz_numerator / lz_denominator[v] ;
 338
 339     }
 340
 341   /* Now clear up after ourselves */
 342   free(lz_denominator);
 343   free(lz);
 344 }
 345