pintos-os.org Git - pspp/blob - src/language/stats/t-test-paired.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 <math.h>
  20 #include <gsl/gsl_cdf.h>
  21
  22 #include "t-test.h"
  23
  24 #include "math/moments.h"
  25 #include "math/correlation.h"
  26 #include "data/casereader.h"
  27 #include "data/dictionary.h"
  28 #include "data/format.h"
  29 #include "data/variable.h"
  30
  31 #include "output/pivot-table.h"
  32
  33 #include "gettext.h"
  34 #define N_(msgid) msgid
  35 #define _(msgid) gettext (msgid)
  36
  37
  38 struct pair_stats
  39 {
  40   double sum_of_prod;
  41   struct moments *mom0;
  42   const struct variable *var0;
  43
  44   struct moments *mom1;
  45   const struct variable *var1;
  46
  47   struct moments *mom_diff;
  48 };
  49
  50 struct paired_samp
  51 {
  52   struct pair_stats *ps;
  53   size_t n_ps;
  54 };
  55
  56 static void paired_summary (const struct tt *tt, struct paired_samp *os);
  57 static void paired_correlations (const struct tt *tt, struct paired_samp *os);
  58 static void paired_test (const struct tt *tt, const struct paired_samp *os);
  59
  60 void
  61 paired_run (const struct tt *tt, size_t n_pairs, vp *pairs, struct casereader *reader)
  62 {
  63   struct ccase *c;
  64   struct paired_samp ps;
  65   struct casereader *r;
  66
  67   ps.ps = xcalloc (n_pairs, sizeof *ps.ps);
  68   ps.n_ps = n_pairs;
  69   for (size_t i = 0; i < n_pairs; ++i)
  70     {
  71       vp *pair = &pairs[i];
  72       struct pair_stats *pp = &ps.ps[i];
  73       pp->var0 = (*pair)[0];
  74       pp->var1 = (*pair)[1];
  75       pp->mom0 = moments_create (MOMENT_VARIANCE);
  76       pp->mom1 = moments_create (MOMENT_VARIANCE);
  77       pp->mom_diff = moments_create (MOMENT_VARIANCE);
  78     }
  79
  80   r = casereader_clone (reader);
  81   for (; (c = casereader_read (r)); case_unref (c))
  82     {
  83       double w = dict_get_case_weight (tt->dict, c, NULL);
  84
  85       for (int i = 0; i < ps.n_ps; i++)
  86         {
  87           struct pair_stats *pp = &ps.ps[i];
  88           const union value *val0 = case_data (c, pp->var0);
  89           const union value *val1 = case_data (c, pp->var1);
  90           if (var_is_value_missing (pp->var0, val0) & tt->exclude)
  91             continue;
  92
  93           if (var_is_value_missing (pp->var1, val1) & tt->exclude)
  94             continue;
  95
  96           moments_pass_one (pp->mom0, val0->f, w);
  97           moments_pass_one (pp->mom1, val1->f, w);
  98           moments_pass_one (pp->mom_diff, val0->f - val1->f, w);
  99         }
 100     }
 101   casereader_destroy (r);
 102
 103   r = reader;
 104   for (; (c = casereader_read (r)); case_unref (c))
 105     {
 106       double w = dict_get_case_weight (tt->dict, c, NULL);
 107
 108       for (int i = 0; i < ps.n_ps; i++)
 109         {
 110           struct pair_stats *pp = &ps.ps[i];
 111           const union value *val0 = case_data (c, pp->var0);
 112           const union value *val1 = case_data (c, pp->var1);
 113           if (var_is_value_missing (pp->var0, val0) & tt->exclude)
 114             continue;
 115
 116           if (var_is_value_missing (pp->var1, val1) & tt->exclude)
 117             continue;
 118
 119           moments_pass_two (pp->mom0, val0->f, w);
 120           moments_pass_two (pp->mom1, val1->f, w);
 121           moments_pass_two (pp->mom_diff, val0->f - val1->f, w);
 122           pp->sum_of_prod += val0->f * val1->f * w;
 123         }
 124     }
 125   casereader_destroy (r);
 126
 127   paired_summary (tt, &ps);
 128   paired_correlations (tt, &ps);
 129   paired_test (tt, &ps);
 130
 131   /* Clean up */
 132
 133   for (int i = 0; i < ps.n_ps; i++)
 134     {
 135       struct pair_stats *pp = &ps.ps[i];
 136       moments_destroy (pp->mom0);
 137       moments_destroy (pp->mom1);
 138       moments_destroy (pp->mom_diff);
 139     }
 140   free (ps.ps);
 141 }
 142
 143 static void
 144 paired_summary (const struct tt *tt, struct paired_samp *os)
 145 {
 146   struct pivot_table *table = pivot_table_create (
 147     N_("Paired Sample Statistics"));
 148   pivot_table_set_weight_var (table, tt->wv);
 149
 150   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 151                           N_("N"), PIVOT_RC_COUNT,
 152                           N_("Mean"), PIVOT_RC_OTHER,
 153                           N_("Std. Deviation"), PIVOT_RC_OTHER,
 154                           N_("S.E. Mean"), PIVOT_RC_OTHER);
 155
 156   struct pivot_dimension *variables = pivot_dimension_create (
 157     table, PIVOT_AXIS_ROW, N_("Variables"));
 158
 159   for (size_t i = 0; i < os->n_ps; i++)
 160     {
 161       struct pair_stats *pp = &os->ps[i];
 162       struct pivot_category *pair = pivot_category_create_group__ (
 163         variables->root, pivot_value_new_text_format (N_("Pair %zu"), i + 1));
 164
 165       for (int j = 0; j < 2; j++)
 166         {
 167           const struct variable *var = j ? pp->var1 : pp->var0;
 168           const struct moments *mom = j ? pp->mom1 : pp->mom0;
 169           double cc, mean, sigma;
 170           moments_calculate (mom, &cc, &mean, &sigma, NULL, NULL);
 171
 172           int var_idx = pivot_category_create_leaf (
 173             pair, pivot_value_new_variable (var));
 174
 175           double entries[] = { cc, mean, sqrt (sigma), sqrt (sigma / cc) };
 176           for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 177             pivot_table_put2 (table, j, var_idx,
 178                               pivot_value_new_number (entries[j]));
 179         }
 180     }
 181
 182   pivot_table_submit (table);
 183 }
 184
 185
 186 static void
 187 paired_correlations (const struct tt *tt, struct paired_samp *os)
 188 {
 189   struct pivot_table *table = pivot_table_create (
 190     N_("Paired Samples Correlations"));
 191   pivot_table_set_weight_var (table, tt->wv);
 192
 193   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 194                           N_("N"), PIVOT_RC_COUNT,
 195                           N_("Correlation"), PIVOT_RC_CORRELATION,
 196                           N_("Sig."), PIVOT_RC_SIGNIFICANCE);
 197
 198   struct pivot_dimension *pairs = pivot_dimension_create (
 199     table, PIVOT_AXIS_ROW, N_("Pairs"));
 200
 201   for (size_t i = 0; i < os->n_ps; i++)
 202     {
 203       struct pair_stats *pp = &os->ps[i];
 204       struct pivot_category *group = pivot_category_create_group__ (
 205         pairs->root, pivot_value_new_text_format (N_("Pair %zu"), i + 1));
 206
 207       int row = pivot_category_create_leaf (
 208         group, pivot_value_new_text_format (N_("%s & %s"),
 209                                             var_to_string (pp->var0),
 210                                             var_to_string (pp->var1)));
 211
 212       double cc0, mean0, sigma0;
 213       double cc1, mean1, sigma1;
 214       moments_calculate (pp->mom0, &cc0, &mean0, &sigma0, NULL, NULL);
 215       moments_calculate (pp->mom1, &cc1, &mean1, &sigma1, NULL, NULL);
 216       /* If this fails, then we're not dealing with missing values properly */
 217       assert (cc0 == cc1);
 218
 219       double corr = ((pp->sum_of_prod / cc0 - mean0 * mean1)
 220                      / sqrt (sigma0 * sigma1) * cc0 / (cc0 - 1));
 221       double sig = 2.0 * significance_of_correlation (corr, cc0);
 222       double entries[] = { cc0, corr, sig };
 223       for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
 224         pivot_table_put2 (table, i, row, pivot_value_new_number (entries[i]));
 225     }
 226
 227   pivot_table_submit (table);
 228 }
 229
 230
 231 static void
 232 paired_test (const struct tt *tt, const struct paired_samp *os)
 233 {
 234   struct pivot_table *table = pivot_table_create (N_("Paired Samples Test"));
 235   pivot_table_set_weight_var (table, tt->wv);
 236
 237   struct pivot_dimension *statistics = pivot_dimension_create (
 238     table, PIVOT_AXIS_COLUMN, N_("Statistics"));
 239   struct pivot_category *group = pivot_category_create_group (
 240     statistics->root, N_("Paired Differences"),
 241     N_("Mean"), PIVOT_RC_OTHER,
 242     N_("Std. Deviation"), PIVOT_RC_OTHER,
 243     N_("S.E. Mean"), PIVOT_RC_OTHER);
 244   struct pivot_category *interval = pivot_category_create_group__ (
 245     group, pivot_value_new_text_format (
 246       N_("%g%% Confidence Interval of the Difference"),
 247       tt->confidence * 100.0));
 248   pivot_category_create_leaves (interval,
 249                                 N_("Lower"), PIVOT_RC_OTHER,
 250                                 N_("Upper"), PIVOT_RC_OTHER);
 251   pivot_category_create_leaves (statistics->root,
 252                                 N_("t"), PIVOT_RC_OTHER,
 253                                 N_("df"), PIVOT_RC_COUNT,
 254                                 N_("Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE);
 255
 256   struct pivot_dimension *pairs = pivot_dimension_create (
 257     table, PIVOT_AXIS_ROW, N_("Pairs"));
 258
 259   for (size_t i = 0; i < os->n_ps; i++)
 260     {
 261       struct pair_stats *pp = &os->ps[i];
 262       struct pivot_category *group = pivot_category_create_group__ (
 263         pairs->root, pivot_value_new_text_format (N_("Pair %zu"), i + 1));
 264
 265       int row = pivot_category_create_leaf (
 266         group, pivot_value_new_text_format (N_("%s - %s"),
 267                                             var_to_string (pp->var0),
 268                                             var_to_string (pp->var1)));
 269
 270       double cc, mean, sigma;
 271       moments_calculate (pp->mom_diff, &cc, &mean, &sigma, NULL, NULL);
 272
 273       double df = cc - 1.0;
 274
 275       double t = mean * sqrt (cc / sigma);
 276       double se_mean = sqrt (sigma / cc);
 277
 278       double p = gsl_cdf_tdist_P (t, df);
 279       double q = gsl_cdf_tdist_Q (t, df);
 280       double sig = 2.0 * (t > 0 ? q : p);
 281
 282       double t_qinv = gsl_cdf_tdist_Qinv ((1.0 - tt->confidence) / 2.0, df);
 283
 284       double entries[] = {
 285         mean,
 286         sqrt (sigma),
 287         se_mean,
 288         mean - t_qinv * se_mean,
 289         mean + t_qinv * se_mean,
 290         t,
 291         df,
 292         sig,
 293       };
 294       for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
 295         pivot_table_put2 (table, i, row, pivot_value_new_number (entries[i]));
 296
 297     }
 298
 299   pivot_table_submit (table);
 300 }