src/language/stats/glm.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 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/case.h"
  24 #include "data/casegrouper.h"
  25 #include "data/casereader.h"
  26 #include "data/dataset.h"
  27 #include "data/dictionary.h"
  28 #include "data/format.h"
  29 #include "data/value.h"
  30 #include "language/command.h"
  31 #include "language/dictionary/split-file.h"
  32 #include "language/lexer/lexer.h"
  33 #include "language/lexer/value-parser.h"
  34 #include "language/lexer/variable-parser.h"
  35 #include "libpspp/ll.h"
  36 #include "libpspp/message.h"
  37 #include "libpspp/misc.h"
  38 #include "libpspp/taint.h"
  39 #include "linreg/sweep.h"
  40 #include "math/categoricals.h"
  41 #include "math/covariance.h"
  42 #include "math/moments.h"
  43 #include "output/tab.h"
  44
  45 #include "gettext.h"
  46 #define _(msgid) gettext (msgid)
  47
  48 struct glm_spec
  49 {
  50   size_t n_dep_vars;
  51   const struct variable **dep_vars;
  52
  53   size_t n_factor_vars;
  54   const struct variable **factor_vars;
  55
  56   enum mv_class exclude;
  57
  58   /* The weight variable */
  59   const struct variable *wv;
  60
  61   bool intercept;
  62 };
  63
  64 struct glm_workspace
  65 {
  66   double total_ssq;
  67   struct moments *totals;
  68 };
  69
  70 static void output_glm (const struct glm_spec *, const struct glm_workspace *ws);
  71 static void run_glm (const struct glm_spec *cmd, struct casereader *input, const struct dataset *ds);
  72
  73 int
  74 cmd_glm (struct lexer *lexer, struct dataset *ds)
  75 {
  76   const struct dictionary *dict = dataset_dict (ds);
  77   struct glm_spec glm ;
  78   glm.n_dep_vars = 0;
  79   glm.n_factor_vars = 0;
  80   glm.dep_vars = NULL;
  81   glm.factor_vars = NULL;
  82   glm.exclude = MV_ANY;
  83   glm.intercept = true;
  84   glm.wv = dict_get_weight (dict);
  85
  86
  87   if (!parse_variables_const (lexer, dict,
  88                               &glm.dep_vars, &glm.n_dep_vars,
  89                               PV_NO_DUPLICATE | PV_NUMERIC))
  90     goto error;
  91
  92   lex_force_match (lexer, T_BY);
  93
  94   if (!parse_variables_const (lexer, dict,
  95                               &glm.factor_vars, &glm.n_factor_vars,
  96                               PV_NO_DUPLICATE | PV_NUMERIC))
  97     goto error;
  98
  99   if ( glm.n_dep_vars > 1)
 100     {
 101       msg (ME, _("Multivariate analysis is not yet implemented"));
 102       return CMD_FAILURE;
 103     }
 104
 105   struct const_var_set *factors = const_var_set_create_from_array (glm.factor_vars, glm.n_factor_vars);
 106
 107
 108   while (lex_token (lexer) != T_ENDCMD)
 109     {
 110       lex_match (lexer, T_SLASH);
 111
 112       if (lex_match_id (lexer, "MISSING"))
 113         {
 114           lex_match (lexer, T_EQUALS);
 115           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 116             {
 117               if (lex_match_id (lexer, "INCLUDE"))
 118                 {
 119                   glm.exclude = MV_SYSTEM;
 120                 }
 121               else if (lex_match_id (lexer, "EXCLUDE"))
 122                 {
 123                   glm.exclude = MV_ANY;
 124                 }
 125               else
 126                 {
 127                   lex_error (lexer, NULL);
 128                   goto error;
 129                 }
 130             }
 131         }
 132       else if (lex_match_id (lexer, "INTERCEPT"))
 133         {
 134           lex_match (lexer, T_EQUALS);
 135           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 136             {
 137               if (lex_match_id (lexer, "INCLUDE"))
 138                 {
 139                   glm.intercept = true;
 140                 }
 141               else if (lex_match_id (lexer, "EXCLUDE"))
 142                 {
 143                   glm.intercept = false;
 144                 }
 145               else
 146                 {
 147                   lex_error (lexer, NULL);
 148                   goto error;
 149                 }
 150             }
 151         }
 152       else if (lex_match_id (lexer, "DESIGN"))
 153         {
 154           size_t n_des;
 155           const struct variable **des;
 156           lex_match (lexer, T_EQUALS);
 157
 158           parse_const_var_set_vars (lexer, factors, &des, &n_des, 0);
 159         }
 160       else
 161         {
 162           lex_error (lexer, NULL);
 163           goto error;
 164         }
 165     }
 166
 167
 168   {
 169     struct casegrouper *grouper;
 170     struct casereader *group;
 171     bool ok;
 172
 173     grouper = casegrouper_create_splits (proc_open (ds), dict);
 174     while (casegrouper_get_next_group (grouper, &group))
 175       run_glm (&glm, group, ds);
 176     ok = casegrouper_destroy (grouper);
 177     ok = proc_commit (ds) && ok;
 178   }
 179
 180   return CMD_SUCCESS;
 181
 182  error:
 183   return CMD_FAILURE;
 184 }
 185
 186 static  void dump_matrix (const gsl_matrix *m);
 187
 188 static void
 189 run_glm (const struct glm_spec *cmd, struct casereader *input, const struct dataset *ds)
 190 {
 191   int v;
 192   struct taint *taint;
 193   struct dictionary *dict = dataset_dict (ds);
 194   struct casereader *reader;
 195   struct ccase *c;
 196
 197   struct glm_workspace ws;
 198
 199   struct categoricals *cats = categoricals_create (cmd->factor_vars, cmd->n_factor_vars,
 200                                                    cmd->wv, cmd->exclude,
 201                                                    NULL, NULL,
 202                                                    NULL, NULL);
 203
 204   struct covariance *cov = covariance_2pass_create (cmd->n_dep_vars, cmd->dep_vars,
 205                                                cats,
 206                                                cmd->wv, cmd->exclude);
 207
 208
 209   c = casereader_peek (input, 0);
 210   if (c == NULL)
 211     {
 212       casereader_destroy (input);
 213       return;
 214     }
 215   output_split_file_values (ds, c);
 216   case_unref (c);
 217
 218   taint = taint_clone (casereader_get_taint (input));
 219
 220   ws.totals = moments_create (MOMENT_VARIANCE);
 221
 222   bool warn_bad_weight = true;
 223   for (reader = casereader_clone (input);
 224        (c = casereader_read (reader)) != NULL; case_unref (c))
 225     {
 226       double weight = dict_get_case_weight (dict, c, &warn_bad_weight);
 227
 228       for ( v = 0; v < cmd->n_dep_vars; ++v)
 229         moments_pass_one (ws.totals, case_data (c, cmd->dep_vars[v])->f, weight);
 230
 231       covariance_accumulate_pass1 (cov, c);
 232     }
 233   casereader_destroy (reader);
 234
 235   categoricals_done (cats);
 236
 237   for (reader = casereader_clone (input);
 238        (c = casereader_read (reader)) != NULL; case_unref (c))
 239     {
 240       double weight = dict_get_case_weight (dict, c, &warn_bad_weight);
 241
 242       for ( v = 0; v < cmd->n_dep_vars; ++v)
 243         moments_pass_two (ws.totals, case_data (c, cmd->dep_vars[v])->f, weight);
 244
 245       covariance_accumulate_pass2 (cov, c);
 246     }
 247   casereader_destroy (reader);
 248
 249   {
 250     gsl_matrix *cm = covariance_calculate_unnormalized (cov);
 251
 252     dump_matrix (cm);
 253
 254     ws.total_ssq = gsl_matrix_get (cm, 0, 0);
 255
 256     reg_sweep (cm, 0);
 257
 258     dump_matrix (cm);
 259
 260     gsl_matrix_free (cm);
 261   }
 262
 263   if (!taint_has_tainted_successor (taint))
 264     output_glm (cmd, &ws);
 265
 266   taint_destroy (taint);
 267 }
 268
 269 static void
 270 output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
 271 {
 272   const struct fmt_spec *wfmt = cmd->wv ? var_get_print_format (cmd->wv) : &F_8_0;
 273
 274   int f;
 275   int r;
 276   const int heading_columns = 1;
 277   const int heading_rows = 1;
 278   struct tab_table *t ;
 279
 280   const int nc = 6;
 281   int nr = heading_rows + 4 + cmd->n_factor_vars;
 282   if (cmd->intercept)
 283     nr++;
 284
 285   t = tab_create (nc, nr);
 286   tab_title (t, _("Tests of Between-Subjects Effects"));
 287
 288   tab_headers (t, heading_columns, 0, heading_rows, 0);
 289
 290   tab_box (t,
 291            TAL_2, TAL_2,
 292            -1, TAL_1,
 293            0, 0,
 294            nc - 1, nr - 1);
 295
 296   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
 297   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
 298
 299   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
 300
 301   /* TRANSLATORS: The parameter is a roman numeral */
 302   tab_text_format (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Type %s Sum of Squares"), "III");
 303   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df"));
 304   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 305   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("F"));
 306   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
 307
 308   r = heading_rows;
 309   tab_text (t, 0, r++, TAB_LEFT | TAT_TITLE, _("Corrected Model"));
 310
 311   double intercept, n_total;
 312   if (cmd->intercept)
 313     {
 314       double mean;
 315       moments_calculate (ws->totals, &n_total, &mean, NULL, NULL, NULL);
 316       intercept = pow2 (mean * n_total) / n_total;
 317
 318       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Intercept"));
 319       tab_double (t, 1, r, 0, intercept, NULL);
 320       tab_double (t, 2, r, 0, 1.00, wfmt);
 321
 322       tab_double (t, 3, r, 0, intercept / 1.0 , NULL);
 323       r++;
 324     }
 325
 326   for (f = 0; f < cmd->n_factor_vars; ++f)
 327     {
 328       tab_text (t, 0, r++, TAB_LEFT | TAT_TITLE,
 329                 var_to_string (cmd->factor_vars[f]));
 330     }
 331
 332   tab_text (t, 0, r++, TAB_LEFT | TAT_TITLE, _("Error"));
 333
 334   if (cmd->intercept)
 335     {
 336       double ssq = intercept + ws->total_ssq;
 337       double mse = ssq / n_total;
 338       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Total"));
 339       tab_double (t, 1, r, 0, ssq, NULL);
 340       tab_double (t, 2, r, 0, n_total, wfmt);
 341
 342       r++;
 343     }
 344
 345   tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Total"));
 346
 347   tab_double (t, 1, r, 0, ws->total_ssq, NULL);
 348   tab_double (t, 2, r, 0, n_total - 1.0, wfmt);
 349
 350   tab_submit (t);
 351 }
 352
 353 static
 354 void dump_matrix (const gsl_matrix *m)
 355 {
 356   size_t i, j;
 357   for (i = 0; i < m->size1; ++i)
 358     {
 359       for (j = 0; j < m->size2; ++j)
 360         {
 361           double x = gsl_matrix_get (m, i, j);
 362           printf ("%.3f ", x);
 363         }
 364       printf ("\n");
 365     }
 366   printf ("\n");
 367 }