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/design-matrix.h>
  42 #include <math/coefficient.h>
  43 #include <math/linreg/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 1000
  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   const struct variable *v;
  71 };
  72
  73
  74 /*
  75   Dependent variable used.
  76  */
  77 static const struct variable **v_dependent;
  78
  79 /*
  80   Number of dependent variables.
  81  */
  82 static size_t n_dependent;
  83
  84 #if 0
  85 /*
  86   Return value for the procedure.
  87  */
  88 static int pspp_glm_rc = CMD_SUCCESS;
  89 #else
  90 int cmd_glm (struct lexer *lexer, struct dataset *ds);
  91 #endif
  92
  93 static bool run_glm (struct casereader*,
  94                      struct cmd_glm *,
  95                      const struct dataset *,
  96                      pspp_linreg_cache *);
  97
  98 int
  99 cmd_glm (struct lexer *lexer, struct dataset *ds)
 100 {
 101   struct casegrouper *grouper;
 102   struct casereader *group;
 103   pspp_linreg_cache *model = NULL;
 104
 105   bool ok;
 106
 107   model = xmalloc (sizeof *model);
 108
 109   if (!parse_glm (lexer, ds, &cmd, NULL))
 110     return CMD_FAILURE;
 111
 112   /* Data pass. */
 113   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 114   while (casegrouper_get_next_group (grouper, &group))
 115     {
 116       run_glm (group, &cmd, ds, model);
 117     }
 118   ok = casegrouper_destroy (grouper);
 119   ok = proc_commit (ds) && ok;
 120
 121   free (model);
 122   free (v_dependent);
 123   return ok ? CMD_SUCCESS : CMD_FAILURE;
 124 }
 125
 126 /* Parser for the dependent sub command */
 127 static int
 128 glm_custom_dependent (struct lexer *lexer, struct dataset *ds,
 129                       struct cmd_glm *cmd UNUSED,
 130                       void *aux UNUSED)
 131 {
 132   const struct dictionary *dict = dataset_dict (ds);
 133
 134   if ((lex_token (lexer) != T_ID
 135        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 136       && lex_token (lexer) != T_ALL)
 137     return 2;
 138
 139   if (!parse_variables_const
 140       (lexer, dict, &v_dependent, &n_dependent, PV_NONE))
 141     {
 142       free (v_dependent);
 143       return 0;
 144     }
 145   assert (n_dependent);
 146   if (n_dependent > 1)
 147     msg (SE, _("Multivariate GLM not yet supported"));
 148   n_dependent = 1; /* Drop this line after adding support for multivariate GLM. */
 149
 150   return 1;
 151 }
 152
 153 static void
 154 coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
 155 {
 156   c->coeff = xnmalloc (dm->m->size2 + 1, sizeof (*c->coeff));
 157   c->coeff[0] = xmalloc (sizeof (*(c->coeff[0])));      /* The first coefficient is the intercept. */
 158   c->coeff[0]->v_info = NULL;   /* Intercept has no associated variable. */
 159   pspp_coeff_init (c->coeff + 1, dm);
 160 }
 161
 162 /*
 163   Put the moments in the linreg cache.
 164  */
 165 static void
 166 compute_moments (pspp_linreg_cache * c, struct moments_var *mom,
 167                  struct design_matrix *dm, size_t n)
 168 {
 169   size_t i;
 170   size_t j;
 171   double weight;
 172   double mean;
 173   double variance;
 174   double skewness;
 175   double kurtosis;
 176   /*
 177      Scan the variable names in the columns of the design matrix.
 178      When we find the variable we need, insert its mean in the cache.
 179    */
 180   for (i = 0; i < dm->m->size2; i++)
 181     {
 182       for (j = 0; j < n; j++)
 183         {
 184           if (design_matrix_col_to_var (dm, i) == (mom + j)->v)
 185             {
 186               moments1_calculate ((mom + j)->m, &weight, &mean, &variance,
 187                                   &skewness, &kurtosis);
 188               gsl_vector_set (c->indep_means, i, mean);
 189               gsl_vector_set (c->indep_std, i, sqrt (variance));
 190             }
 191         }
 192     }
 193 }
 194 /* Encode categorical variables.
 195    Returns number of valid cases. */
 196 static int
 197 prepare_categories (struct casereader *input,
 198                     const struct variable **vars, size_t n_vars,
 199                     struct moments_var *mom)
 200 {
 201   int n_data;
 202   struct ccase c;
 203   size_t i;
 204
 205   for (i = 0; i < n_vars; i++)
 206     if (var_is_alpha (vars[i]))
 207       cat_stored_values_create (vars[i]);
 208
 209   n_data = 0;
 210   for (; casereader_read (input, &c); case_destroy (&c))
 211     {
 212       /*
 213         The second condition ensures the program will run even if
 214         there is only one variable to act as both explanatory and
 215         response.
 216        */
 217       for (i = 0; i < n_vars; i++)
 218         {
 219           const union value *val = case_data (&c, vars[i]);
 220           if (var_is_alpha (vars[i]))
 221             cat_value_update (vars[i], val);
 222           else
 223             moments1_add (mom[i].m, val->f, 1.0);
 224         }
 225       n_data++;
 226    }
 227   casereader_destroy (input);
 228
 229   return n_data;
 230 }
 231
 232 static bool
 233 run_glm (struct casereader *input,
 234          struct cmd_glm *cmd,
 235          const struct dataset *ds,
 236          pspp_linreg_cache *model)
 237 {
 238   size_t i;
 239   int n_indep = 0;
 240   struct ccase c;
 241   const struct variable **indep_vars;
 242   struct design_matrix *X;
 243   struct moments_var *mom;
 244   gsl_vector *Y;
 245   struct casereader *reader;
 246   casenumber row;
 247   size_t n_data;                /* Number of valid cases. */
 248
 249   pspp_linreg_opts lopts;
 250
 251   assert (model != NULL);
 252
 253   if (!casereader_peek (input, 0, &c))
 254     {
 255       casereader_destroy (input);
 256       return true;
 257     }
 258   output_split_file_values (ds, &c);
 259   case_destroy (&c);
 260
 261   if (!v_dependent)
 262     {
 263       dict_get_vars (dataset_dict (ds), &v_dependent, &n_dependent,
 264                      1u << DC_SYSTEM);
 265     }
 266
 267   for (i = 0; i < n_dependent; i++)
 268     {
 269       if (!var_is_numeric (v_dependent[i]))
 270         {
 271           msg (SE, _("Dependent variable must be numeric."));
 272           return false;
 273         }
 274     }
 275
 276   mom = xnmalloc (n_dependent, sizeof (*mom));
 277   mom->m = moments1_create (MOMENT_VARIANCE);
 278   mom->v = v_dependent[0];
 279   lopts.get_depvar_mean_std = 1;
 280
 281   lopts.get_indep_mean_std = xnmalloc (n_dependent, sizeof (int));
 282   indep_vars = xnmalloc (cmd->n_by, sizeof *indep_vars);
 283
 284   for (i = 0; i < cmd->n_by; i++)
 285     {
 286       indep_vars[i] = cmd->v_by[i];
 287     }
 288   n_indep = cmd->n_by;
 289
 290   reader = casereader_clone (input);
 291   reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
 292                                              MV_ANY, NULL);
 293   reader = casereader_create_filter_missing (reader, v_dependent, 1,
 294                                              MV_ANY, NULL);
 295   n_data = prepare_categories (casereader_clone (reader),
 296                                indep_vars, n_indep, mom);
 297
 298   if ((n_data > 0) && (n_indep > 0))
 299     {
 300       Y = gsl_vector_alloc (n_data);
 301       X =
 302         design_matrix_create (n_indep,
 303                               (const struct variable **) indep_vars,
 304                               n_data);
 305       for (i = 0; i < X->m->size2; i++)
 306         {
 307           lopts.get_indep_mean_std[i] = 1;
 308         }
 309       model = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
 310       model->indep_means = gsl_vector_alloc (X->m->size2);
 311       model->indep_std = gsl_vector_alloc (X->m->size2);
 312       model->depvar = v_dependent[0];
 313       reader = casereader_create_counter (reader, &row, -1);
 314       for (; casereader_read (reader, &c); case_destroy (&c))
 315         {
 316           for (i = 0; i < n_indep; ++i)
 317             {
 318               const struct variable *v = indep_vars[i];
 319               const union value *val = case_data (&c, v);
 320               if (var_is_alpha (v))
 321                 design_matrix_set_categorical (X, row, v, val);
 322               else
 323                 design_matrix_set_numeric (X, row, v, val);
 324             }
 325           gsl_vector_set (Y, row, case_num (&c, model->depvar));
 326         }
 327       casereader_destroy (reader);
 328       /*
 329         Now that we know the number of coefficients, allocate space
 330         and store pointers to the variables that correspond to the
 331         coefficients.
 332       */
 333       coeff_init (model, X);
 334
 335       /*
 336         Find the least-squares estimates and other statistics.
 337       */
 338       pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, model);
 339       compute_moments (model, mom, X, n_dependent);
 340
 341       gsl_vector_free (Y);
 342       design_matrix_destroy (X);
 343     }
 344   else
 345     {
 346       msg (SE, gettext ("No valid data found. This command was skipped."));
 347     }
 348   free (indep_vars);
 349   free (lopts.get_indep_mean_std);
 350   casereader_destroy (input);
 351
 352   return true;
 353 }
 354
 355 /*
 356   Local Variables:
 357   mode: c
 358   End:
 359 */