src/language/stats/glm.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2007, 2009 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/hash.h>
  41 #include <libpspp/message.h>
  42 #include <math/covariance.h>
  43 #include <math/coefficient.h>
  44 #include <math/linreg.h>
  45 #include <math/moments.h>
  46 #include <output/table.h>
  47
  48 #include "xalloc.h"
  49 #include "gettext.h"
  50
  51 /* (headers) */
  52
  53 /* (specification)
  54    "GLM" (glm_):
  55    *dependent=custom;
  56    design=custom;
  57    by=varlist;
  58    with=varlist.
  59 */
  60 /* (declarations) */
  61 /* (functions) */
  62 static struct cmd_glm cmd;
  63
  64
  65 /*
  66   Moments for each of the variables used.
  67  */
  68 struct moments_var
  69 {
  70   struct moments1 *m;
  71   double *weight;
  72   double *mean;
  73   double *variance;
  74   const struct variable *v;
  75 };
  76
  77
  78 /*
  79   Dependent variable used.
  80  */
  81 static const struct variable **v_dependent;
  82
  83 /*
  84   Number of dependent variables.
  85  */
  86 static size_t n_dependent;
  87
  88 size_t n_inter; /* Number of interactions. */
  89 size_t n_members; /* Number of memebr variables in an interaction. */
  90
  91 struct interaction_variable **interactions;
  92
  93 int cmd_glm (struct lexer *lexer, struct dataset *ds);
  94
  95 static bool run_glm (struct casereader *,
  96                      struct cmd_glm *,
  97                      const struct dataset *);
  98 /*
  99   If the DESIGN subcommand was not specified, specify all possible
 100   two-way interactions.
 101  */
 102 static void
 103 check_interactions (struct dataset *ds, struct cmd_glm *cmd)
 104 {
 105   size_t i;
 106   size_t j;
 107   size_t k = 0;
 108   struct variable **interaction_vars;
 109
 110   /*
 111      User did not specify the design matrix, so we
 112      specify it here.
 113   */
 114   n_inter = (cmd->n_by + cmd->n_with) * (cmd->n_by + cmd->n_with - 1) / 2;
 115   interactions = xnmalloc (n_inter, sizeof (*interactions));
 116   interaction_vars = xnmalloc (2, sizeof (*interaction_vars));
 117   for (i = 0; i < cmd->n_by; i++)
 118     {
 119       for (j = i + 1; j < cmd->n_by; j++)
 120         {
 121           interaction_vars[0] = cmd->v_by[i];
 122           interaction_vars[1] = cmd->v_by[j];
 123           interactions[k] = interaction_variable_create (interaction_vars, 2);
 124           k++;
 125         }
 126       for (j = 0; j < cmd->n_with; j++)
 127         {
 128           interaction_vars[0] = cmd->v_by[i];
 129           interaction_vars[1] = cmd->v_with[j];
 130           interactions[k] = interaction_variable_create (interaction_vars, 2);
 131           k++;
 132         }
 133     }
 134   for (i = 0; i < cmd->n_with; i++)
 135     {
 136       for (j = i + 1; j < cmd->n_with; j++)
 137         {
 138           interaction_vars[0] = cmd->v_with[i];
 139           interaction_vars[1] = cmd->v_with[j];
 140           interactions[k] = interaction_variable_create (interaction_vars, 2);
 141           k++;
 142         }
 143     }
 144 }
 145 int
 146 cmd_glm (struct lexer *lexer, struct dataset *ds)
 147 {
 148   struct casegrouper *grouper;
 149   struct casereader *group;
 150
 151   bool ok;
 152
 153   if (!parse_glm (lexer, ds, &cmd, NULL))
 154     return CMD_FAILURE;
 155
 156   if (!lex_match_id (lexer, "DESIGN"))
 157     {
 158       check_interactions (ds, &cmd);
 159     }
 160    /* Data pass. */
 161   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 162   while (casegrouper_get_next_group (grouper, &group))
 163     {
 164       run_glm (group, &cmd, ds);
 165     }
 166   ok = casegrouper_destroy (grouper);
 167   ok = proc_commit (ds) && ok;
 168
 169   free (v_dependent);
 170   return ok ? CMD_SUCCESS : CMD_FAILURE;
 171 }
 172
 173 static int
 174 parse_interactions (struct lexer *lexer, const struct variable **interaction_vars, int n_members,
 175                     int max_members, struct dataset *ds)
 176 {
 177   if (lex_match (lexer, '*'))
 178     {
 179       if (n_members > max_members)
 180         {
 181           max_members *= 2;
 182           xnrealloc (interaction_vars, max_members, sizeof (*interaction_vars));
 183         }
 184       interaction_vars[n_members] = parse_variable (lexer, dataset_dict (ds));
 185       parse_interactions (lexer, interaction_vars, n_members++, max_members, ds);
 186     }
 187   return n_members;
 188 }
 189 /* Parser for the design subcommand. */
 190 static int
 191 glm_custom_design (struct lexer *lexer, struct dataset *ds,
 192                    struct cmd_glm *cmd UNUSED, void *aux UNUSED)
 193 {
 194   size_t n_allocated = 2;
 195   size_t n_members;
 196   struct variable **interaction_vars;
 197   struct variable *this_var;
 198
 199   interactions = xnmalloc (n_allocated, sizeof (*interactions));
 200   n_inter = 1;
 201   while (lex_token (lexer) != T_STOP && lex_token (lexer) != '.')
 202     {
 203       this_var = parse_variable (lexer, dataset_dict (ds));
 204       if (lex_match (lexer, '('))
 205         {
 206           lex_force_match (lexer, ')');
 207         }
 208       else if (lex_match (lexer, '*'))
 209         {
 210           interaction_vars = xnmalloc (2 * n_inter, sizeof (*interaction_vars));
 211           n_members = parse_interactions (lexer, interaction_vars, 1, 2 * n_inter, ds);
 212           if (n_allocated < n_inter)
 213             {
 214               n_allocated *= 2;
 215               xnrealloc (interactions, n_allocated, sizeof (*interactions));
 216             }
 217           interactions [n_inter - 1] =
 218             interaction_variable_create (interaction_vars, n_members);
 219           n_inter++;
 220           free (interaction_vars);
 221         }
 222     }
 223   return 1;
 224 }
 225 /* Parser for the dependent sub command */
 226 static int
 227 glm_custom_dependent (struct lexer *lexer, struct dataset *ds,
 228                       struct cmd_glm *cmd UNUSED, void *aux UNUSED)
 229 {
 230   const struct dictionary *dict = dataset_dict (ds);
 231
 232   if ((lex_token (lexer) != T_ID
 233        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 234       && lex_token (lexer) != T_ALL)
 235     return 2;
 236
 237   if (!parse_variables_const
 238       (lexer, dict, &v_dependent, &n_dependent, PV_NONE))
 239     {
 240       free (v_dependent);
 241       return 0;
 242     }
 243   assert (n_dependent);
 244   if (n_dependent > 1)
 245     msg (SE, _("Multivariate GLM not yet supported"));
 246   n_dependent = 1;              /* Drop this line after adding support for multivariate GLM. */
 247
 248   return 1;
 249 }
 250
 251 /*
 252   COV is the covariance matrix for variables included in the
 253   model. That means the dependent variable is in there, too.
 254  */
 255 static void
 256 coeff_init (pspp_linreg_cache * c, const struct design_matrix *cov)
 257 {
 258   c->coeff = xnmalloc (cov->m->size2, sizeof (*c->coeff));
 259   c->n_coeffs = cov->m->size2 - 1;
 260   pspp_coeff_init (c->coeff, cov);
 261 }
 262
 263
 264 static pspp_linreg_cache *
 265 fit_model (const struct covariance *cov,
 266            const struct variable *dep_var,
 267            const struct variable ** indep_vars,
 268            size_t n_data, size_t n_indep)
 269 {
 270   pspp_linreg_cache *result = NULL;
 271
 272   return result;
 273 }
 274
 275 static bool
 276 run_glm (struct casereader *input,
 277          struct cmd_glm *cmd,
 278          const struct dataset *ds)
 279 {
 280   casenumber row;
 281   const struct variable **numerics = NULL;
 282   const struct variable **categoricals = NULL;
 283   int n_indep = 0;
 284   pspp_linreg_cache *model = NULL;
 285   pspp_linreg_opts lopts;
 286   struct ccase *c;
 287   size_t i;
 288   size_t n_data;                /* Number of valid cases. */
 289   size_t n_categoricals = 0;
 290   size_t n_numerics;
 291   struct casereader *reader;
 292   struct covariance *cov;
 293
 294   c = casereader_peek (input, 0);
 295   if (c == NULL)
 296     {
 297       casereader_destroy (input);
 298       return true;
 299     }
 300   output_split_file_values (ds, c);
 301   case_unref (c);
 302
 303   if (!v_dependent)
 304     {
 305       dict_get_vars (dataset_dict (ds), &v_dependent, &n_dependent,
 306                      1u << DC_SYSTEM);
 307     }
 308
 309   lopts.get_depvar_mean_std = 1;
 310
 311   lopts.get_indep_mean_std = xnmalloc (n_dependent, sizeof (int));
 312
 313
 314   n_numerics = cmd->n_with + n_dependent;
 315   for (i = 0; i < cmd->n_by; i++)
 316     {
 317       if (var_is_alpha (cmd->v_by[i]))
 318         {
 319           n_categoricals++;
 320         }
 321       else
 322         {
 323           n_numerics++;
 324         }
 325     }
 326   numerics = xnmalloc (n_categoricals, sizeof *numerics);
 327   categoricals = xnmalloc (n_categoricals, sizeof (*categoricals));
 328   size_t k = 0;
 329   for (i = 0; i < cmd->n_by; i++)
 330     {
 331       if (var_is_alpha (cmd->v_by[i]))
 332         {
 333           categoricals[k] = cmd->v_by[i];
 334         }
 335       else
 336         {
 337           numerics[i] = cmd->v_by[i];
 338           k++;
 339         }
 340     }
 341   for (i = 0; i < n_dependent; i++)
 342     {
 343       k++;
 344       numerics[k] = v_dependent[i];
 345     }
 346   for (i = 0; i < cmd->n_with; i++)
 347     {
 348       k++;
 349       numerics[k] = v_dependent[i];
 350     }
 351
 352   covariance_2pass_create (n_numerics, numerics, n_categoricals, categoricals, NULL, MV_NEVER);
 353
 354   reader = casereader_clone (input);
 355   reader = casereader_create_filter_missing (reader, numerics, n_numerics,
 356                                              MV_ANY, NULL, NULL);
 357   reader = casereader_create_filter_missing (reader, categoricals, n_categoricals,
 358                                              MV_ANY, NULL, NULL);
 359   struct casereader *r = casereader_clone (reader);
 360
 361   reader = casereader_create_counter (reader, &row, -1);
 362
 363   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
 364     {
 365       covariance_accumulate_pass1 (cov, c);
 366     }
 367   for (; (c = casereader_read (r)) != NULL; case_unref (c))
 368     {
 369       covariance_accumulate_pass1 (cov, c);
 370     }
 371   covariance_destroy (cov);
 372   casereader_destroy (reader);
 373   casereader_destroy (r);
 374
 375   free (numerics);
 376   free (categoricals);
 377   free (lopts.get_indep_mean_std);
 378   casereader_destroy (input);
 379
 380   return true;
 381 }
 382
 383 /*
 384   Local Variables:
 385   mode: c
 386   End:
 387 */