pintos-os.org Git - pspp/blob - src/language/stats/t-test-indep.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2011 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
  19 #include "t-test.h"
  20
  21 #include <gsl/gsl_cdf.h>
  22 #include <math.h>
  23 #include "libpspp/misc.h"
  24
  25 #include "libpspp/str.h"
  26 #include "data/casereader.h"
  27 #include "data/dictionary.h"
  28 #include "data/format.h"
  29 #include "data/variable.h"
  30
  31 #include "math/moments.h"
  32 #include "math/levene.h"
  33
  34 #include <output/tab.h>
  35 #include "gettext.h"
  36 #define _(msgid) gettext (msgid)
  37
  38
  39 struct indep_samples
  40 {
  41   const struct variable *gvar;
  42   bool cut;
  43   const union value *gval0;
  44   const union value *gval1;
  45 };
  46
  47 struct pair_stats
  48 {
  49   struct moments *mom[2];
  50   double lev ;
  51   struct levene *nl;
  52 };
  53
  54
  55 static void indep_summary (const struct tt *tt, struct indep_samples *is, const struct pair_stats *ps);
  56 static void indep_test (const struct tt *tt, const struct pair_stats *ps);
  57
  58 static int
  59 which_group (const union value *v, const struct indep_samples *is)
  60 {
  61   int width = var_get_width (is->gvar);
  62   int cmp = value_compare_3way (v, is->gval0, width);
  63   if ( is->cut )
  64     return  (cmp < 0);
  65
  66   if (cmp == 0)
  67     return 0;
  68
  69   if (0 == value_compare_3way (v, is->gval1, width))
  70     return 1;
  71
  72   return -1;
  73 }
  74
  75 void
  76 indep_run (struct tt *tt, const struct variable *gvar,
  77            bool cut,
  78            const union value *gval0, const union value *gval1,
  79            struct casereader *reader)
  80 {
  81   struct indep_samples is;
  82   struct ccase *c;
  83   struct casereader *r;
  84
  85   struct pair_stats *ps = xcalloc (sizeof (*ps), tt->n_vars);
  86
  87   int v;
  88
  89   for (v = 0; v < tt->n_vars; ++v)
  90     {
  91       ps[v].mom[0] = moments_create (MOMENT_VARIANCE);
  92       ps[v].mom[1] = moments_create (MOMENT_VARIANCE);
  93       ps[v].nl = levene_create (var_get_width (gvar), cut ? gval0: NULL);
  94     }
  95
  96   is.gvar = gvar;
  97   is.gval0 = gval0;
  98   is.gval1 = gval1;
  99   is.cut = cut;
 100
 101   r = casereader_clone (reader);
 102   for ( ; (c = casereader_read (r) ); case_unref (c))
 103     {
 104       double w = dict_get_case_weight (tt->dict, c, NULL);
 105
 106       const union value *gv = case_data (c, gvar);
 107
 108       int grp = which_group (gv, &is);
 109       if ( grp < 0)
 110         continue;
 111
 112       for (v = 0; v < tt->n_vars; ++v)
 113         {
 114           const union value *val = case_data (c, tt->vars[v]);
 115           if (var_is_value_missing (tt->vars[v], val, tt->exclude))
 116             continue;
 117
 118           moments_pass_one (ps[v].mom[grp], val->f, w);
 119           levene_pass_one (ps[v].nl, val->f, w, gv);
 120         }
 121     }
 122   casereader_destroy (r);
 123
 124   r = casereader_clone (reader);
 125   for ( ; (c = casereader_read (r) ); case_unref (c))
 126     {
 127       double w = dict_get_case_weight (tt->dict, c, NULL);
 128
 129       const union value *gv = case_data (c, gvar);
 130
 131       int grp = which_group (gv, &is);
 132       if ( grp < 0)
 133         continue;
 134
 135       for (v = 0; v < tt->n_vars; ++v)
 136         {
 137           const union value *val = case_data (c, tt->vars[v]);
 138           if (var_is_value_missing (tt->vars[v], val, tt->exclude))
 139             continue;
 140
 141           moments_pass_two (ps[v].mom[grp], val->f, w);
 142           levene_pass_two (ps[v].nl, val->f, w, gv);
 143         }
 144     }
 145   casereader_destroy (r);
 146
 147   r = reader;
 148   for ( ; (c = casereader_read (r) ); case_unref (c))
 149     {
 150       double w = dict_get_case_weight (tt->dict, c, NULL);
 151
 152       const union value *gv = case_data (c, gvar);
 153
 154       int grp = which_group (gv, &is);
 155       if ( grp < 0)
 156         continue;
 157
 158       for (v = 0; v < tt->n_vars; ++v)
 159         {
 160           const union value *val = case_data (c, tt->vars[v]);
 161           if (var_is_value_missing (tt->vars[v], val, tt->exclude))
 162             continue;
 163
 164           levene_pass_three (ps[v].nl, val->f, w, gv);
 165         }
 166     }
 167   casereader_destroy (r);
 168
 169
 170   for (v = 0; v < tt->n_vars; ++v)
 171     ps[v].lev = levene_calculate (ps[v].nl);
 172
 173   indep_summary (tt, &is, ps);
 174   indep_test (tt, ps);
 175
 176
 177   for (v = 0; v < tt->n_vars; ++v)
 178     {
 179       moments_destroy (ps[v].mom[0]);
 180       moments_destroy (ps[v].mom[1]);
 181       levene_destroy (ps[v].nl);
 182     }
 183   free (ps);
 184 }
 185
 186
 187 static void
 188 indep_summary (const struct tt *tt, struct indep_samples *is, const struct pair_stats *ps)
 189 {
 190   const struct fmt_spec *wfmt = tt->wv ? var_get_print_format (tt->wv) : & F_8_0;
 191
 192   int v;
 193   int cols = 6;
 194   const int heading_rows = 1;
 195   int rows = tt->n_vars * 2 + heading_rows;
 196
 197   struct string vallab0 ;
 198   struct string vallab1 ;
 199   struct tab_table *t = tab_create (cols, rows);
 200
 201   ds_init_empty (&vallab0);
 202   ds_init_empty (&vallab1);
 203
 204   tab_headers (t, 0, 0, 1, 0);
 205   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols - 1, rows - 1);
 206   tab_hline (t, TAL_2, 0, cols - 1, 1);
 207
 208   tab_vline (t, TAL_GAP, 1, 0, rows - 1);
 209   tab_title (t, _("Group Statistics"));
 210   tab_text  (t, 1, 0, TAB_CENTER | TAT_TITLE, var_to_string (is->gvar));
 211   tab_text  (t, 2, 0, TAB_CENTER | TAT_TITLE, _("N"));
 212   tab_text  (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
 213   tab_text  (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 214   tab_text  (t, 5, 0, TAB_CENTER | TAT_TITLE, _("S.E. Mean"));
 215
 216   if (is->cut)
 217     {
 218       ds_put_cstr (&vallab0, "≥");
 219       ds_put_cstr (&vallab1, "<");
 220
 221       var_append_value_name (is->gvar, is->gval0, &vallab0);
 222       var_append_value_name (is->gvar, is->gval0, &vallab1);
 223     }
 224   else
 225     {
 226       var_append_value_name (is->gvar, is->gval0, &vallab0);
 227       var_append_value_name (is->gvar, is->gval1, &vallab1);
 228     }
 229
 230   for (v = 0; v < tt->n_vars; ++v)
 231     {
 232       int i;
 233       const struct variable *var = tt->vars[v];
 234
 235       tab_text (t, 0, v * 2 + heading_rows, TAB_LEFT,
 236                 var_to_string (var));
 237
 238       tab_text (t, 1, v * 2 + heading_rows, TAB_LEFT,
 239                        ds_cstr (&vallab0));
 240
 241       tab_text (t, 1, v * 2 + 1 + heading_rows, TAB_LEFT,
 242                        ds_cstr (&vallab1));
 243
 244       for (i = 0 ; i < 2; ++i)
 245         {
 246           double cc, mean, sigma;
 247           moments_calculate (ps[v].mom[i], &cc, &mean, &sigma, NULL, NULL);
 248
 249           tab_double (t, 2, v * 2 + i + heading_rows, TAB_RIGHT, cc, wfmt);
 250           tab_double (t, 3, v * 2 + i + heading_rows, TAB_RIGHT, mean, NULL);
 251           tab_double (t, 4, v * 2 + i + heading_rows, TAB_RIGHT, sqrt (sigma), NULL);
 252           tab_double (t, 5, v * 2 + i + heading_rows, TAB_RIGHT, sqrt (sigma / cc), NULL);
 253         }
 254     }
 255
 256   tab_submit (t);
 257
 258   ds_destroy (&vallab0);
 259   ds_destroy (&vallab1);
 260 }
 261
 262
 263 static void
 264 indep_test (const struct tt *tt, const struct pair_stats *ps)
 265 {
 266   int v;
 267   const int heading_rows = 3;
 268   const int rows= tt->n_vars * 2 + heading_rows;
 269
 270   const size_t cols = 11;
 271
 272   struct tab_table *t = tab_create (cols, rows);
 273   tab_headers (t, 0, 0, 3, 0);
 274   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols - 1, rows - 1);
 275   tab_hline (t, TAL_2, 0, cols - 1, 3);
 276
 277   tab_title (t, _("Independent Samples Test"));
 278
 279   tab_hline (t, TAL_1, 2, cols - 1, 1);
 280   tab_vline (t, TAL_2, 2, 0, rows - 1);
 281   tab_vline (t, TAL_1, 4, 0, rows - 1);
 282   tab_box (t, -1, -1, -1, TAL_1, 2, 1, cols - 2, rows - 1);
 283   tab_hline (t, TAL_1, cols - 2, cols - 1, 2);
 284   tab_box (t, -1, -1, -1, TAL_1, cols - 2, 2, cols - 1, rows - 1);
 285   tab_joint_text (t, 2, 0, 3, 0, TAB_CENTER, _("Levene's Test for Equality of Variances"));
 286   tab_joint_text (t, 4, 0, cols - 1, 0, TAB_CENTER, _("t-test for Equality of Means"));
 287
 288   tab_text (t, 2, 2, TAB_CENTER | TAT_TITLE, _("F"));
 289   tab_text (t, 3, 2, TAB_CENTER | TAT_TITLE, _("Sig."));
 290   tab_text (t, 4, 2, TAB_CENTER | TAT_TITLE, _("t"));
 291   tab_text (t, 5, 2, TAB_CENTER | TAT_TITLE, _("df"));
 292   tab_text (t, 6, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
 293   tab_text (t, 7, 2, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
 294   tab_text (t, 8, 2, TAB_CENTER | TAT_TITLE, _("Std. Error Difference"));
 295   tab_text (t, 9, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
 296   tab_text (t, 10, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
 297
 298   tab_joint_text_format (t, 9, 1, 10, 1, TAB_CENTER,
 299                          _("%g%% Confidence Interval of the Difference"),
 300                          tt->confidence * 100.0);
 301
 302   for (v = 0; v < tt->n_vars; ++v)
 303   {
 304     double df, pooled_variance, mean_diff, tval;
 305     double se2, std_err_diff;
 306     double p, q;
 307     double cc0, mean0, sigma0;
 308     double cc1, mean1, sigma1;
 309     moments_calculate (ps[v].mom[0], &cc0, &mean0, &sigma0, NULL, NULL);
 310     moments_calculate (ps[v].mom[1], &cc1, &mean1, &sigma1, NULL, NULL);
 311
 312     tab_text (t, 0, v * 2 + heading_rows, TAB_LEFT, var_to_string (tt->vars[v]));
 313     tab_text (t, 1, v * 2 + heading_rows, TAB_LEFT, _("Equal variances assumed"));
 314
 315     df = cc0 + cc1 - 2.0;
 316     tab_double (t, 5, v * 2 + heading_rows, TAB_RIGHT, df, NULL);
 317
 318     pooled_variance = ((cc0 - 1)* sigma0 + (cc1 - 1) * sigma1) / df ;
 319
 320     tval = (mean0 - mean1) / sqrt (pooled_variance);
 321     tval /= sqrt ((cc0 + cc1) / (cc0 * cc1));
 322
 323     tab_double (t, 4, v * 2 + heading_rows, TAB_RIGHT, tval, NULL);
 324
 325     p = gsl_cdf_tdist_P (tval, df);
 326     q = gsl_cdf_tdist_Q (tval, df);
 327
 328     mean_diff = mean0 - mean1;
 329
 330     tab_double (t, 6, v * 2 + heading_rows, TAB_RIGHT, 2.0 * (tval > 0 ? q : p),   NULL);
 331     tab_double (t, 7, v * 2 + heading_rows, TAB_RIGHT, mean_diff, NULL);
 332
 333     std_err_diff = sqrt ((sigma0 / cc0) + (sigma1 / cc1));
 334     tab_double (t, 8, v * 2 + heading_rows, TAB_RIGHT, std_err_diff, NULL);
 335
 336
 337     /* Now work out the confidence interval */
 338     q = (1 - tt->confidence)/2.0;  /* 2-tailed test */
 339
 340     tval = gsl_cdf_tdist_Qinv (q, df);
 341     tab_double (t,  9, v * 2 + heading_rows, TAB_RIGHT, mean_diff - tval * std_err_diff, NULL);
 342     tab_double (t, 10, v * 2 + heading_rows, TAB_RIGHT, mean_diff + tval * std_err_diff, NULL);
 343
 344     /* Equal variances not assumed */
 345     tab_text (t, 1, v * 2 + heading_rows + 1,  TAB_LEFT, _("Equal variances not assumed"));
 346
 347     se2 = sigma0 / cc0 + sigma1 / cc1;
 348     tval = mean_diff / sqrt (se2);
 349     tab_double (t, 4, v * 2 + heading_rows + 1, TAB_RIGHT, tval, NULL);
 350
 351     {
 352       double p, q;
 353       const double s0 = sigma0 / (cc0);
 354       const double s1 = sigma1 / (cc1);
 355       double df = pow2 (s0 + s1) ;
 356       df /= pow2 (s0) / (cc0 - 1) + pow2 (s1) / (cc1 - 1);
 357
 358       tab_double (t, 5, v * 2 + heading_rows + 1, TAB_RIGHT, df, NULL);
 359
 360       p = gsl_cdf_tdist_P (tval, df);
 361       q = gsl_cdf_tdist_Q (tval, df);
 362
 363       tab_double (t, 6, v * 2 + heading_rows + 1, TAB_RIGHT, 2.0 * (tval > 0 ? q : p), NULL);
 364
 365       /* Now work out the confidence interval */
 366       q = (1 - tt->confidence) / 2.0;  /* 2-tailed test */
 367
 368       tval = gsl_cdf_tdist_Qinv (q, df);
 369     }
 370
 371     tab_double (t, 7, v * 2 + heading_rows + 1, TAB_RIGHT, mean_diff, NULL);
 372     tab_double (t, 8, v * 2 + heading_rows + 1, TAB_RIGHT, std_err_diff, NULL);
 373     tab_double (t, 9, v * 2 + heading_rows + 1, TAB_RIGHT,  mean_diff - tval * std_err_diff, NULL);
 374     tab_double (t, 10, v * 2 + heading_rows + 1, TAB_RIGHT, mean_diff + tval * std_err_diff, NULL);
 375
 376     tab_double (t, 2, v * 2 + heading_rows, TAB_CENTER, ps[v].lev, NULL);
 377
 378
 379     {
 380       /* Now work out the significance of the Levene test */
 381       double df1 = 1;
 382       double df2 = cc0 + cc1 - 2;
 383       double q = gsl_cdf_fdist_Q (ps[v].lev, df1, df2);
 384       tab_double (t, 3, v * 2 + heading_rows, TAB_CENTER, q, NULL);
 385     }
 386   }
 387
 388   tab_submit (t);
 389 }