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