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/casereader.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 <libpspp/alloc.h>
  33 #include <libpspp/misc.h>
  34 #include "group.h"
  35
  36 #include <math.h>
  37 #include <stdlib.h>
  38
  39
  40 /* This module calculates the Levene statistic for variables.
  41
  42    Just for reference, the Levene Statistic is a defines as follows:
  43
  44    W = \frac{ (n-k)\sum_{i=1}^k n_i(Z_{iL} - Z_{LL})^2}
  45             { (k-1)\sum_{i=1}^k \sum_{j=1}^{n_i} (Z_{ij} - Z_{iL})^2}
  46
  47    where:
  48         k is the number of groups
  49         n is the total number of samples
  50         n_i is the number of samples in the ith group
  51         Z_{ij} is | Y_{ij} - Y_{iL} | where Y_{iL} is the mean of the ith group
  52         Z_{iL} is the  mean of Z_{ij} over the ith group
  53         Z_{LL} is the grand mean of Z_{ij}
  54
  55    Imagine calculating that with pencil and paper!
  56
  57  */
  58
  59
  60 struct levene_info
  61 {
  62
  63   /* Per group statistics */
  64   struct t_test_proc **group_stats;
  65
  66   /* The independent variable */
  67   const struct variable *v_indep;
  68
  69   /* Number of dependent variables */
  70   size_t n_dep;
  71
  72   /* The dependent variables */
  73   const struct variable  **v_dep;
  74
  75   /* Filter for missing values */
  76   enum mv_class exclude;
  77
  78   /* An array of lz_stats for each variable */
  79   struct lz_stats *lz;
  80
  81   /* The denominator for the expression for the Levene */
  82   double *lz_denominator;
  83
  84 };
  85
  86 /* Per variable statistics */
  87 struct lz_stats
  88 {
  89   /* Total of all lz */
  90   double grand_total;
  91
  92   /* Mean of all lz */
  93   double grand_mean;
  94
  95   /* The total number of cases */
  96   double total_n ;
  97
  98   /* Number of groups */
  99   int n_groups;
 100 };
 101
 102 /* First pass */
 103 static void  levene_precalc (const struct levene_info *l);
 104 static int levene_calc (const struct dictionary *dict, const struct ccase *,
 105                         const struct levene_info *l);
 106 static void levene_postcalc (struct levene_info *);
 107
 108
 109 /* Second pass */
 110 static void levene2_precalc (struct levene_info *l);
 111 static int levene2_calc (const struct dictionary *, const struct ccase *,
 112                          struct levene_info *l);
 113 static void levene2_postcalc (struct levene_info *);
 114
 115
 116 void
 117 levene(const struct dictionary *dict,
 118        struct casereader *reader,
 119        const struct variable *v_indep, size_t n_dep,
 120        const struct variable **v_dep,
 121        enum mv_class exclude)
 122 {
 123   struct casereader *pass1, *pass2;
 124   struct ccase c;
 125   struct levene_info l;
 126
 127   l.n_dep      = n_dep;
 128   l.v_indep    = v_indep;
 129   l.v_dep      = v_dep;
 130   l.exclude    = exclude;
 131   l.lz         = xnmalloc (l.n_dep, sizeof *l.lz);
 132   l.lz_denominator = xnmalloc (l.n_dep, sizeof *l.lz_denominator);
 133
 134   casereader_split (reader, &pass1, &pass2);
 135
 136   levene_precalc (&l);
 137   for (; casereader_read (pass1, &c); case_destroy (&c))
 138     levene_calc (dict, &c, &l);
 139   casereader_destroy (pass1);
 140   levene_postcalc (&l);
 141
 142   levene2_precalc(&l);
 143   for (; casereader_read (pass2, &c); case_destroy (&c))
 144     levene2_calc (dict, &c, &l);
 145   casereader_destroy (pass2);
 146   levene2_postcalc (&l);
 147
 148   free (l.lz_denominator);
 149   free (l.lz);
 150 }
 151
 152 static void
 153 levene_precalc (const struct levene_info *l)
 154 {
 155   size_t i;
 156
 157   for(i = 0; i < l->n_dep ; ++i )
 158     {
 159       const struct variable *var = l->v_dep[i];
 160       struct group_proc *gp = group_proc_get (var);
 161       struct group_statistics *gs;
 162       struct hsh_iterator hi;
 163
 164       l->lz[i].grand_total = 0;
 165       l->lz[i].total_n = 0;
 166       l->lz[i].n_groups = gp->n_groups ;
 167
 168
 169       for ( gs = hsh_first(gp->group_hash, &hi);
 170             gs != 0;
 171             gs = hsh_next(gp->group_hash, &hi))
 172         {
 173           gs->lz_total = 0;
 174         }
 175
 176     }
 177
 178 }
 179
 180 static int
 181 levene_calc (const struct dictionary *dict, const struct ccase *c,
 182              const struct levene_info *l)
 183 {
 184   size_t i;
 185   bool warn = false;
 186   const union value *gv = case_data (c, l->v_indep);
 187   struct group_statistics key;
 188   double weight = dict_get_case_weight (dict, c, &warn);
 189
 190   key.id = *gv;
 191
 192   for (i = 0; i < l->n_dep; ++i)
 193     {
 194       const struct variable *var = l->v_dep[i];
 195       struct group_proc *gp = group_proc_get (var);
 196       double levene_z;
 197       const union value *v = case_data (c, var);
 198       struct group_statistics *gs;
 199
 200       gs = hsh_find(gp->group_hash,(void *) &key );
 201
 202       if ( 0 == gs )
 203         continue ;
 204
 205       if ( !var_is_value_missing (var, v, l->exclude))
 206         {
 207           levene_z= fabs(v->f - gs->mean);
 208           l->lz[i].grand_total += levene_z * weight;
 209           l->lz[i].total_n += weight;
 210
 211           gs->lz_total += levene_z * weight;
 212         }
 213     }
 214   return 0;
 215 }
 216
 217
 218 static void
 219 levene_postcalc (struct levene_info *l)
 220 {
 221   size_t v;
 222
 223   for (v = 0; v < l->n_dep; ++v)
 224     {
 225       /* This is Z_LL */
 226       l->lz[v].grand_mean = l->lz[v].grand_total / l->lz[v].total_n ;
 227     }
 228
 229
 230 }
 231
 232
 233
 234 static void
 235 levene2_precalc (struct levene_info *l)
 236 {
 237   size_t v;
 238
 239
 240   /* This stuff could go in the first post calc . . . */
 241   for (v = 0;
 242        v < l->n_dep;
 243        ++v)
 244     {
 245       struct hsh_iterator hi;
 246       struct group_statistics *g;
 247
 248       const struct variable *var = l->v_dep[v] ;
 249       struct hsh_table *hash = group_proc_get (var)->group_hash;
 250
 251
 252       for(g = (struct group_statistics *) hsh_first(hash,&hi);
 253           g != 0 ;
 254           g = (struct group_statistics *) hsh_next(hash,&hi) )
 255         {
 256           g->lz_mean = g->lz_total / g->n ;
 257         }
 258       l->lz_denominator[v] = 0;
 259   }
 260 }
 261
 262 static int
 263 levene2_calc (const struct dictionary *dict, const struct ccase *c,
 264               struct levene_info *l)
 265 {
 266   size_t i;
 267   bool warn = false;
 268
 269   double weight = dict_get_case_weight (dict, c, &warn);
 270
 271   const union value *gv = case_data (c, l->v_indep);
 272   struct group_statistics key;
 273
 274   key.id = *gv;
 275
 276   for (i = 0; i < l->n_dep; ++i)
 277     {
 278       double levene_z;
 279       const struct variable *var = l->v_dep[i] ;
 280       const union value *v = case_data (c, var);
 281       struct group_statistics *gs;
 282
 283       gs = hsh_find(group_proc_get (var)->group_hash,(void *) &key );
 284
 285       if ( 0 == gs )
 286         continue;
 287
 288       if ( !var_is_value_missing (var, v, l->exclude))
 289         {
 290           levene_z = fabs(v->f - gs->mean);
 291           l->lz_denominator[i] += weight * pow2 (levene_z - gs->lz_mean);
 292         }
 293     }
 294
 295   return 0;
 296 }
 297
 298
 299 static void
 300 levene2_postcalc (struct levene_info *l)
 301 {
 302   size_t v;
 303
 304   for (v = 0; v < l->n_dep; ++v)
 305     {
 306       double lz_numerator = 0;
 307       struct hsh_iterator hi;
 308       struct group_statistics *g;
 309
 310       const struct variable *var = l->v_dep[v] ;
 311       struct group_proc *gp = group_proc_get (var);
 312       struct hsh_table *hash = gp->group_hash;
 313
 314       for(g = (struct group_statistics *) hsh_first(hash,&hi);
 315           g != 0 ;
 316           g = (struct group_statistics *) hsh_next(hash,&hi) )
 317         {
 318           lz_numerator += g->n * pow2(g->lz_mean - l->lz[v].grand_mean );
 319         }
 320       lz_numerator *= ( gp->ugs.n - gp->n_groups );
 321
 322       l->lz_denominator[v] *= (gp->n_groups - 1);
 323
 324       gp->levene = lz_numerator / l->lz_denominator[v] ;
 325
 326     }
 327 }
 328