src/math/levene.c

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