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