pintos-os.org Git - pspp/blob - src/language/stats/correlations.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2009, 2010, 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 <gsl/gsl_cdf.h>
  20 #include <gsl/gsl_matrix.h>
  21 #include <math.h>
  22
  23 #include "data/casegrouper.h"
  24 #include "data/casereader.h"
  25 #include "data/dataset.h"
  26 #include "data/dictionary.h"
  27 #include "data/format.h"
  28 #include "data/variable.h"
  29 #include "language/command.h"
  30 #include "language/dictionary/split-file.h"
  31 #include "language/lexer/lexer.h"
  32 #include "language/lexer/variable-parser.h"
  33 #include "libpspp/assertion.h"
  34 #include "libpspp/message.h"
  35 #include "libpspp/misc.h"
  36 #include "math/correlation.h"
  37 #include "math/covariance.h"
  38 #include "math/moments.h"
  39 #include "output/pivot-table.h"
  40
  41 #include "gl/xalloc.h"
  42 #include "gl/minmax.h"
  43
  44 #include "gettext.h"
  45 #define _(msgid) gettext (msgid)
  46 #define N_(msgid) msgid
  47
  48
  49 struct corr
  50   {
  51     size_t n_vars_total;
  52     size_t n_vars1;
  53
  54     const struct variable **vars;
  55   };
  56
  57
  58 /* Handling of missing values. */
  59 enum corr_missing_type
  60   {
  61     CORR_PAIRWISE,    /* Handle missing values on a per-variable-pair basis. */
  62     CORR_LISTWISE     /* Discard entire case if any variable is missing. */
  63   };
  64
  65 struct corr_opts
  66 {
  67   enum corr_missing_type missing_type;
  68   enum mv_class exclude;      /* Classes of missing values to exclude. */
  69
  70   bool sig;   /* Flag significant values or not */
  71   int tails;  /* Report significance with how many tails ? */
  72   bool descriptive_stats;
  73   bool xprod_stats;
  74
  75   const struct variable *wv;  /* The weight variable (if any) */
  76 };
  77
  78
  79 static void
  80 output_descriptives (const struct corr *corr, const struct corr_opts *opts,
  81                      const gsl_matrix *means,
  82                      const gsl_matrix *vars, const gsl_matrix *ns)
  83 {
  84   struct pivot_table *table = pivot_table_create (
  85     N_("Descriptive Statistics"));
  86   pivot_table_set_weight_var (table, opts->wv);
  87
  88   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
  89                           N_("Mean"), PIVOT_RC_OTHER,
  90                           N_("Std. Deviation"), PIVOT_RC_OTHER,
  91                           N_("N"), PIVOT_RC_COUNT);
  92
  93   struct pivot_dimension *variables = pivot_dimension_create (
  94     table, PIVOT_AXIS_ROW, N_("Variable"));
  95
  96   for (size_t r = 0; r < corr->n_vars_total; ++r)
  97     {
  98       const struct variable *v = corr->vars[r];
  99
 100       int row = pivot_category_create_leaf (variables->root,
 101                                             pivot_value_new_variable (v));
 102
 103       double mean = gsl_matrix_get (means, r, 0);
 104       /* Here we want to display the non-biased estimator */
 105       double n = gsl_matrix_get (ns, r, 0);
 106       double stddev = sqrt (gsl_matrix_get (vars, r, 0) * n / (n - 1));
 107       double entries[] = { mean, stddev, n };
 108       for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
 109         pivot_table_put2 (table, i, row, pivot_value_new_number (entries[i]));
 110     }
 111
 112   pivot_table_submit (table);
 113 }
 114
 115 static void
 116 output_correlation (const struct corr *corr, const struct corr_opts *opts,
 117                     const gsl_matrix *cm, const gsl_matrix *samples,
 118                     const gsl_matrix *cv)
 119 {
 120   struct pivot_table *table = pivot_table_create (N_("Correlations"));
 121   pivot_table_set_weight_var (table, opts->wv);
 122
 123   /* Column variable dimension. */
 124   struct pivot_dimension *columns = pivot_dimension_create (
 125     table, PIVOT_AXIS_COLUMN, N_("Variables"));
 126
 127   size_t matrix_cols = (corr->n_vars_total > corr->n_vars1
 128                         ? corr->n_vars_total - corr->n_vars1
 129                         : corr->n_vars1);
 130   for (size_t c = 0; c < matrix_cols; c++)
 131     {
 132       const struct variable *v = corr->n_vars_total > corr->n_vars1 ?
 133         corr->vars[corr->n_vars1 + c] : corr->vars[c];
 134       pivot_category_create_leaf (columns->root, pivot_value_new_variable (v));
 135     }
 136
 137   /* Statistics dimension. */
 138   struct pivot_dimension *statistics = pivot_dimension_create (
 139     table, PIVOT_AXIS_ROW, N_("Statistics"),
 140     N_("Pearson Correlation"), PIVOT_RC_CORRELATION,
 141     opts->tails == 2 ? N_("Sig. (2-tailed)") : N_("Sig. (1-tailed)"),
 142     PIVOT_RC_SIGNIFICANCE);
 143
 144   if (opts->xprod_stats)
 145     pivot_category_create_leaves (statistics->root, N_("Cross-products"),
 146                                   N_("Covariance"));
 147
 148   if (opts->missing_type != CORR_LISTWISE)
 149     pivot_category_create_leaves (statistics->root, N_("N"), PIVOT_RC_COUNT);
 150
 151   /* Row variable dimension. */
 152   struct pivot_dimension *rows = pivot_dimension_create (
 153     table, PIVOT_AXIS_ROW, N_("Variables"));
 154   for (size_t r = 0; r < corr->n_vars1; r++)
 155     pivot_category_create_leaf (rows->root,
 156                                 pivot_value_new_variable (corr->vars[r]));
 157
 158   struct pivot_footnote *sig_footnote = pivot_table_create_footnote (
 159     table, pivot_value_new_text (N_("Significant at .05 level")));
 160
 161   for (size_t r = 0; r < corr->n_vars1; r++)
 162     for (size_t c = 0; c < matrix_cols; c++)
 163       {
 164         const int col_index = (corr->n_vars_total > corr->n_vars1
 165                                ? corr->n_vars1 + c
 166                                : c);
 167         double pearson = gsl_matrix_get (cm, r, col_index);
 168         double w = gsl_matrix_get (samples, r, col_index);
 169         double sig = opts->tails * significance_of_correlation (pearson, w);
 170
 171         double entries[5];
 172         int n = 0;
 173         entries[n++] = pearson;
 174         entries[n++] = col_index != r ? sig : SYSMIS;
 175         if (opts->xprod_stats)
 176           {
 177             double cov = gsl_matrix_get (cv, r, col_index);
 178             const double xprod_dev = cov * w;
 179             cov *= w / (w - 1.0);
 180
 181             entries[n++] = xprod_dev;
 182             entries[n++] = cov;
 183           }
 184         if (opts->missing_type != CORR_LISTWISE)
 185           entries[n++] = w;
 186
 187         for (int i = 0; i < n; i++)
 188           if (entries[i] != SYSMIS)
 189             {
 190               struct pivot_value *v = pivot_value_new_number (entries[i]);
 191               if (!i && opts->sig && col_index != r && sig < 0.05)
 192                 pivot_value_add_footnote (v, sig_footnote);
 193               pivot_table_put3 (table, c, i, r, v);
 194             }
 195       }
 196
 197   pivot_table_submit (table);
 198 }
 199
 200
 201 static void
 202 run_corr (struct casereader *r, const struct corr_opts *opts, const struct corr *corr)
 203 {
 204   struct covariance *cov = covariance_2pass_create (
 205     corr->n_vars_total, corr->vars, NULL,opts->wv, opts->exclude, true);
 206
 207   struct casereader *rc = casereader_clone (r);
 208   struct ccase *c;
 209   for (; (c = casereader_read (r)); case_unref (c))
 210     covariance_accumulate_pass1 (cov, c);
 211   for (; (c = casereader_read (rc)); case_unref (c))
 212     covariance_accumulate_pass2 (cov, c);
 213   casereader_destroy (rc);
 214
 215   gsl_matrix *cov_matrix = covariance_calculate (cov);
 216   if (!cov_matrix)
 217     {
 218       msg (SE, _("The data for the chosen variables are all missing or empty."));
 219       covariance_destroy (cov);
 220       return;
 221     }
 222
 223   const gsl_matrix *samples_matrix = covariance_moments (cov, MOMENT_NONE);
 224   const gsl_matrix *var_matrix = covariance_moments (cov, MOMENT_VARIANCE);
 225   const gsl_matrix *mean_matrix = covariance_moments (cov, MOMENT_MEAN);
 226
 227   gsl_matrix *corr_matrix = correlation_from_covariance (cov_matrix, var_matrix);
 228
 229   if (opts->descriptive_stats)
 230     output_descriptives (corr, opts, mean_matrix, var_matrix, samples_matrix);
 231
 232   output_correlation (corr, opts, corr_matrix, samples_matrix, cov_matrix);
 233
 234   covariance_destroy (cov);
 235   gsl_matrix_free (corr_matrix);
 236   gsl_matrix_free (cov_matrix);
 237 }
 238
 239 int
 240 cmd_correlations (struct lexer *lexer, struct dataset *ds)
 241 {
 242   size_t n_all_vars = 0; /* Total number of variables involved in this command */
 243   const struct dictionary *dict = dataset_dict (ds);
 244
 245   struct corr *corrs = NULL;
 246   size_t n_corrs = 0;
 247   size_t allocated_corrs = 0;
 248
 249   struct corr_opts opts = {
 250     .missing_type = CORR_PAIRWISE,
 251     .wv = dict_get_weight (dict),
 252     .tails = 2,
 253     .exclude = MV_ANY,
 254   };
 255
 256   /* Parse CORRELATIONS. */
 257   while (lex_token (lexer) != T_ENDCMD)
 258     {
 259       lex_match (lexer, T_SLASH);
 260       if (lex_match_id (lexer, "MISSING"))
 261         {
 262           lex_match (lexer, T_EQUALS);
 263           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 264             {
 265               if (lex_match_id (lexer, "PAIRWISE"))
 266                 opts.missing_type = CORR_PAIRWISE;
 267               else if (lex_match_id (lexer, "LISTWISE"))
 268                 opts.missing_type = CORR_LISTWISE;
 269               else if (lex_match_id (lexer, "INCLUDE"))
 270                 opts.exclude = MV_SYSTEM;
 271               else if (lex_match_id (lexer, "EXCLUDE"))
 272                 opts.exclude = MV_ANY;
 273               else
 274                 {
 275                   lex_error_expecting (lexer, "PAIRWISE", "LISTWISE",
 276                                        "INCLUDE", "EXCLUDE");
 277                   goto error;
 278                 }
 279               lex_match (lexer, T_COMMA);
 280             }
 281         }
 282       else if (lex_match_id (lexer, "PRINT"))
 283         {
 284           lex_match (lexer, T_EQUALS);
 285           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 286             {
 287               if (lex_match_id (lexer, "TWOTAIL"))
 288                 opts.tails = 2;
 289               else if (lex_match_id (lexer, "ONETAIL"))
 290                 opts.tails = 1;
 291               else if (lex_match_id (lexer, "SIG"))
 292                 opts.sig = false;
 293               else if (lex_match_id (lexer, "NOSIG"))
 294                 opts.sig = true;
 295               else
 296                 {
 297                   lex_error_expecting (lexer, "TWOTAIL", "ONETAIL",
 298                                        "SIG", "NOSIG");
 299                   goto error;
 300                 }
 301
 302               lex_match (lexer, T_COMMA);
 303             }
 304         }
 305       else if (lex_match_id (lexer, "STATISTICS"))
 306         {
 307           lex_match (lexer, T_EQUALS);
 308           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 309             {
 310               if (lex_match_id (lexer, "DESCRIPTIVES"))
 311                 opts.descriptive_stats = true;
 312               else if (lex_match_id (lexer, "XPROD"))
 313                 opts.xprod_stats = true;
 314               else if (lex_token (lexer) == T_ALL)
 315                 {
 316                   opts.descriptive_stats = opts.xprod_stats = true;
 317                   lex_get (lexer);
 318                 }
 319               else
 320                 {
 321                   lex_error_expecting (lexer, "DESCRIPTIVES", "XPROD", "ALL");
 322                   goto error;
 323                 }
 324
 325               lex_match (lexer, T_COMMA);
 326             }
 327         }
 328       else
 329         {
 330           if (lex_match_id (lexer, "VARIABLES"))
 331             lex_match (lexer, T_EQUALS);
 332
 333           const struct variable **vars;
 334           size_t n_vars1;
 335           if (!parse_variables_const (lexer, dict, &vars, &n_vars1, PV_NUMERIC))
 336             goto error;
 337
 338           size_t n_vars_total = n_vars1;
 339           if (lex_match (lexer, T_WITH)
 340               && !parse_variables_const (lexer, dict, &vars, &n_vars_total,
 341                                          PV_NUMERIC | PV_APPEND))
 342             goto error;
 343
 344           if (n_corrs >= allocated_corrs)
 345             corrs = x2nrealloc (corrs, &allocated_corrs, sizeof *corrs);
 346           corrs[n_corrs++] = (struct corr) {
 347             .n_vars1 = n_vars1,
 348             .n_vars_total = n_vars_total,
 349             .vars = vars,
 350           };
 351
 352           n_all_vars += n_vars_total;
 353         }
 354     }
 355   if (n_corrs == 0)
 356     {
 357       lex_ofs_error (lexer, 0, lex_ofs (lexer) - 1,
 358                      _("No variables specified."));
 359       goto error;
 360     }
 361
 362   const struct variable **all_vars = xmalloc (n_all_vars * sizeof *all_vars);
 363   const struct variable **vv = all_vars;
 364   for (size_t i = 0; i < n_corrs; ++i)
 365     {
 366       const struct corr *c = &corrs[i];
 367       for (size_t v = 0; v < c->n_vars_total; ++v)
 368         *vv++ = c->vars[v];
 369     }
 370
 371   struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
 372   struct casereader *group;
 373   while (casegrouper_get_next_group (grouper, &group))
 374     {
 375       for (size_t i = 0; i < n_corrs; ++i)
 376         {
 377           /* FIXME: No need to iterate the data multiple times */
 378           struct casereader *r = casereader_clone (group);
 379
 380           if (opts.missing_type == CORR_LISTWISE)
 381             r = casereader_create_filter_missing (r, all_vars, n_all_vars,
 382                                                   opts.exclude, NULL, NULL);
 383
 384
 385           run_corr (r, &opts, &corrs[i]);
 386           casereader_destroy (r);
 387         }
 388       casereader_destroy (group);
 389     }
 390   bool ok = casegrouper_destroy (grouper);
 391   ok = proc_commit (ds) && ok;
 392
 393   free (all_vars);
 394
 395   /* Done. */
 396   for (size_t i = 0; i < n_corrs; i++)
 397     free (corrs[i].vars);
 398   free (corrs);
 399
 400   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 401
 402 error:
 403   for (size_t i = 0; i < n_corrs; i++)
 404     free (corrs[i].vars);
 405   free (corrs);
 406   return CMD_FAILURE;
 407 }