src/language/stats/glm.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2007 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 <gsl/gsl_vector.h>
  22 #include <math.h>
  23 #include <stdlib.h>
  24
  25 #include <data/case.h>
  26 #include <data/category.h>
  27 #include <data/casegrouper.h>
  28 #include <data/casereader.h>
  29 #include <data/dictionary.h>
  30 #include <data/missing-values.h>
  31 #include <data/procedure.h>
  32 #include <data/transformations.h>
  33 #include <data/value-labels.h>
  34 #include <data/variable.h>
  35 #include <language/command.h>
  36 #include <language/dictionary/split-file.h>
  37 #include <language/data-io/file-handle.h>
  38 #include <language/lexer/lexer.h>
  39 #include <libpspp/compiler.h>
  40 #include <libpspp/message.h>
  41 #include <math/covariance-matrix.h>
  42 #include <math/coefficient.h>
  43 #include <math/linreg.h>
  44 #include <math/moments.h>
  45 #include <output/table.h>
  46
  47 #include "xalloc.h"
  48 #include "gettext.h"
  49
  50 #define GLM_LARGE_DATA 10000
  51
  52 /* (headers) */
  53
  54 /* (specification)
  55    "GLM" (glm_):
  56    *dependent=custom;
  57    by=varlist;
  58    with=varlist.
  59 */
  60 /* (declarations) */
  61 /* (functions) */
  62 static struct cmd_glm cmd;
  63
  64 /*
  65   Moments for each of the variables used.
  66  */
  67 struct moments_var
  68 {
  69   struct moments1 *m;
  70   double *weight;
  71   double *mean;
  72   double *variance;
  73   const struct variable *v;
  74 };
  75
  76
  77 /*
  78   Dependent variable used.
  79  */
  80 static const struct variable **v_dependent;
  81
  82 /*
  83   Number of dependent variables.
  84  */
  85 static size_t n_dependent;
  86
  87 #if 0
  88 /*
  89   Return value for the procedure.
  90  */
  91 static int pspp_glm_rc = CMD_SUCCESS;
  92 #else
  93 int cmd_glm (struct lexer *lexer, struct dataset *ds);
  94 #endif
  95
  96 static bool run_glm (struct casereader *,
  97                      struct cmd_glm *,
  98                      const struct dataset *, pspp_linreg_cache *);
  99
 100 int
 101 cmd_glm (struct lexer *lexer, struct dataset *ds)
 102 {
 103   struct casegrouper *grouper;
 104   struct casereader *group;
 105   pspp_linreg_cache *model = NULL;
 106
 107   bool ok;
 108
 109   model = xmalloc (sizeof *model);
 110
 111   if (!parse_glm (lexer, ds, &cmd, NULL))
 112     return CMD_FAILURE;
 113
 114   /* Data pass. */
 115   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 116   while (casegrouper_get_next_group (grouper, &group))
 117     {
 118       run_glm (group, &cmd, ds, model);
 119     }
 120   ok = casegrouper_destroy (grouper);
 121   ok = proc_commit (ds) && ok;
 122
 123   free (model);
 124   free (v_dependent);
 125   return ok ? CMD_SUCCESS : CMD_FAILURE;
 126 }
 127
 128 /* Parser for the dependent sub command */
 129 static int
 130 glm_custom_dependent (struct lexer *lexer, struct dataset *ds,
 131                       struct cmd_glm *cmd UNUSED, void *aux UNUSED)
 132 {
 133   const struct dictionary *dict = dataset_dict (ds);
 134
 135   if ((lex_token (lexer) != T_ID
 136        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 137       && lex_token (lexer) != T_ALL)
 138     return 2;
 139
 140   if (!parse_variables_const
 141       (lexer, dict, &v_dependent, &n_dependent, PV_NONE))
 142     {
 143       free (v_dependent);
 144       return 0;
 145     }
 146   assert (n_dependent);
 147   if (n_dependent > 1)
 148     msg (SE, _("Multivariate GLM not yet supported"));
 149   n_dependent = 1;              /* Drop this line after adding support for multivariate GLM. */
 150
 151   return 1;
 152 }
 153
 154 /*
 155   COV is the covariance matrix for variables included in the
 156   model. That means the dependent variable is in there, too.
 157  */
 158 static void
 159 coeff_init (pspp_linreg_cache * c, struct design_matrix *cov)
 160 {
 161   c->coeff = xnmalloc (cov->m->size2, sizeof (*c->coeff));
 162   c->n_coeffs = cov->m->size2 - 1;
 163   c->coeff[0] = xmalloc (sizeof (*(c->coeff[0])));      /* The first coefficient is the intercept. */
 164   c->coeff[0]->v_info = NULL;   /* Intercept has no associated variable. */
 165   pspp_coeff_init (c->coeff + 1, cov);
 166 }
 167
 168 /*
 169   Put the moments in the linreg cache.
 170  */
 171 static void
 172 compute_moments (pspp_linreg_cache * c, struct moments_var *mom,
 173                  struct design_matrix *dm, size_t n)
 174 {
 175   size_t i;
 176   size_t j;
 177   double weight;
 178   double mean;
 179   double variance;
 180   double skewness;
 181   double kurtosis;
 182   /*
 183      Scan the variable names in the columns of the design matrix.
 184      When we find the variable we need, insert its mean in the cache.
 185    */
 186   for (i = 0; i < dm->m->size2; i++)
 187     {
 188       for (j = 0; j < n; j++)
 189         {
 190           if (design_matrix_col_to_var (dm, i) == (mom + j)->v)
 191             {
 192               moments1_calculate ((mom + j)->m, &weight, &mean, &variance,
 193                                   &skewness, &kurtosis);
 194               gsl_vector_set (c->indep_means, i, mean);
 195               gsl_vector_set (c->indep_std, i, sqrt (variance));
 196             }
 197         }
 198     }
 199 }
 200
 201 /* Encode categorical variables.
 202    Returns number of valid cases. */
 203 static int
 204 data_pass_one (struct casereader *input,
 205                const struct variable **vars, size_t n_vars,
 206                struct moments_var **mom)
 207 {
 208   int n_data;
 209   struct ccase c;
 210   size_t i;
 211
 212   for (i = 0; i < n_vars; i++)
 213     {
 214       mom[i] = xmalloc (sizeof (*mom[i]));
 215       mom[i]->v = vars[i];
 216       mom[i]->mean = xmalloc (sizeof (*mom[i]->mean));
 217       mom[i]->variance = xmalloc (sizeof (*mom[i]->mean));
 218       mom[i]->weight = xmalloc (sizeof (*mom[i]->weight));
 219       mom[i]->m = moments1_create (MOMENT_VARIANCE);
 220       if (var_is_alpha (vars[i]))
 221         cat_stored_values_create (vars[i]);
 222     }
 223
 224   n_data = 0;
 225   for (; casereader_read (input, &c); case_destroy (&c))
 226     {
 227       /*
 228          The second condition ensures the program will run even if
 229          there is only one variable to act as both explanatory and
 230          response.
 231        */
 232       for (i = 0; i < n_vars; i++)
 233         {
 234           const union value *val = case_data (&c, vars[i]);
 235           if (var_is_alpha (vars[i]))
 236             cat_value_update (vars[i], val);
 237           else
 238             moments1_add (mom[i]->m, val->f, 1.0);
 239         }
 240       n_data++;
 241     }
 242   casereader_destroy (input);
 243   for (i = 0; i < n_vars; i++)
 244     {
 245       if (var_is_numeric (mom[i]->v))
 246         {
 247           moments1_calculate (mom[i]->m, mom[i]->weight, mom[i]->mean,
 248                               mom[i]->variance, NULL, NULL);
 249         }
 250     }
 251
 252   return n_data;
 253 }
 254
 255 static bool
 256 run_glm (struct casereader *input,
 257          struct cmd_glm *cmd,
 258          const struct dataset *ds, pspp_linreg_cache * model)
 259 {
 260   size_t i;
 261   size_t j;
 262   int n_indep = 0;
 263   struct ccase c;
 264   const struct variable **indep_vars;
 265   const struct variable **all_vars;
 266   struct design_matrix *X;
 267   struct moments_var **mom;
 268   struct casereader *reader;
 269   casenumber row;
 270   size_t n_all_vars;
 271   size_t n_data;                /* Number of valid cases. */
 272
 273   pspp_linreg_opts lopts;
 274
 275   assert (model != NULL);
 276
 277   if (!casereader_peek (input, 0, &c))
 278     {
 279       casereader_destroy (input);
 280       return true;
 281     }
 282   output_split_file_values (ds, &c);
 283   case_destroy (&c);
 284
 285   if (!v_dependent)
 286     {
 287       dict_get_vars (dataset_dict (ds), &v_dependent, &n_dependent,
 288                      1u << DC_SYSTEM);
 289     }
 290
 291
 292
 293   lopts.get_depvar_mean_std = 1;
 294
 295   lopts.get_indep_mean_std = xnmalloc (n_dependent, sizeof (int));
 296   indep_vars = xnmalloc (cmd->n_by, sizeof *indep_vars);
 297   n_all_vars = cmd->n_by + n_dependent;
 298   all_vars = xnmalloc (n_all_vars, sizeof *all_vars);
 299
 300   for (i = 0; i < n_dependent; i++)
 301     {
 302       all_vars[i] = v_dependent[i];
 303     }
 304   for (i = 0; i < cmd->n_by; i++)
 305     {
 306       indep_vars[i] = cmd->v_by[i];
 307       all_vars[i + n_dependent] = cmd->v_by[i];
 308     }
 309   n_indep = cmd->n_by;
 310   mom = xnmalloc (n_all_vars, sizeof (*mom));
 311
 312
 313   reader = casereader_clone (input);
 314   reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
 315                                              MV_ANY, NULL, NULL);
 316   reader = casereader_create_filter_missing (reader, v_dependent, 1,
 317                                              MV_ANY, NULL, NULL);
 318   n_data = data_pass_one (casereader_clone (reader),
 319                           (const struct variable **) all_vars, n_all_vars,
 320                           mom);
 321
 322   if ((n_data > 0) && (n_indep > 0))
 323     {
 324       X =
 325         covariance_matrix_create (n_all_vars,
 326                                   (const struct variable **) all_vars);
 327       reader = casereader_create_counter (reader, &row, -1);
 328       for (; casereader_read (reader, &c); case_destroy (&c))
 329         {
 330           /*
 331              Accumulate the covariance matrix.
 332            */
 333           for (i = 0; i < n_all_vars; ++i)
 334             {
 335               const struct variable *v = all_vars[i];
 336               const union value *val_v = case_data (&c, v);
 337               for (j = i; j < n_all_vars; j++)
 338                 {
 339                   const struct variable *w = all_vars[j];
 340                   const union value *val_w = case_data (&c, w);
 341                   covariance_pass_two (X, *mom[i]->mean, *mom[j]->mean,
 342                                        (double) n_data,
 343                                        v, w, val_v, val_w);
 344                 }
 345             }
 346         }
 347       model = pspp_linreg_cache_alloc (v_dependent, indep_vars, n_data, n_indep);
 348       /*
 349         For large data sets, use QR decomposition.
 350       */
 351       if (n_data > sqrt (n_indep) && n_data > GLM_LARGE_DATA)
 352         {
 353           model->method = PSPP_LINREG_QR;
 354         }
 355       coeff_init (model, X);
 356       pspp_linreg_with_cov (X, model);
 357       casereader_destroy (reader);
 358       for (i = 0; i < n_all_vars; i++)
 359         {
 360           moments1_destroy (mom[i]->m);
 361           free (mom[i]->mean);
 362           free (mom[i]->variance);
 363           free (mom[i]->weight);
 364           free (mom[i]);
 365         }
 366       free (mom);
 367       covariance_matrix_destroy (X);
 368     }
 369   else
 370     {
 371       msg (SE, gettext ("No valid data found. This command was skipped."));
 372     }
 373   free (indep_vars);
 374   free (lopts.get_indep_mean_std);
 375   pspp_linreg_cache_free (model);
 376   casereader_destroy (input);
 377
 378   return true;
 379 }
 380
 381 /*
 382   Local Variables:
 383   mode: c
 384   End:
 385 */