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   /*
  62      Sums of squares due to different variables. Element 0 is the SSE
  63      for the entire model. For i > 0, element i is the SS due to
  64      variable i.
  65    */
  66   gsl_vector * ssq;
  67
  68   bool intercept;
  69 };
  70
  71 struct glm_workspace
  72 {
  73   double total_ssq;
  74   struct moments *totals;
  75 };
  76
  77 static void output_glm (struct glm_spec *, const struct glm_workspace *ws);
  78 static void run_glm (struct glm_spec *cmd, struct casereader *input, const struct dataset *ds);
  79
  80 int
  81 cmd_glm (struct lexer *lexer, struct dataset *ds)
  82 {
  83   const struct dictionary *dict = dataset_dict (ds);
  84   struct glm_spec glm ;
  85   glm.n_dep_vars = 0;
  86   glm.n_factor_vars = 0;
  87   glm.dep_vars = NULL;
  88   glm.factor_vars = NULL;
  89   glm.exclude = MV_ANY;
  90   glm.intercept = true;
  91   glm.wv = dict_get_weight (dict);
  92
  93
  94   if (!parse_variables_const (lexer, dict,
  95                               &glm.dep_vars, &glm.n_dep_vars,
  96                               PV_NO_DUPLICATE | PV_NUMERIC))
  97     goto error;
  98
  99   lex_force_match (lexer, T_BY);
 100
 101   if (!parse_variables_const (lexer, dict,
 102                               &glm.factor_vars, &glm.n_factor_vars,
 103                               PV_NO_DUPLICATE | PV_NUMERIC))
 104     goto error;
 105
 106   if ( glm.n_dep_vars > 1)
 107     {
 108       msg (ME, _("Multivariate analysis is not yet implemented"));
 109       return CMD_FAILURE;
 110     }
 111
 112   struct const_var_set *factors = const_var_set_create_from_array (glm.factor_vars, glm.n_factor_vars);
 113
 114
 115   while (lex_token (lexer) != T_ENDCMD)
 116     {
 117       lex_match (lexer, T_SLASH);
 118
 119       if (lex_match_id (lexer, "MISSING"))
 120         {
 121           lex_match (lexer, T_EQUALS);
 122           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 123             {
 124               if (lex_match_id (lexer, "INCLUDE"))
 125                 {
 126                   glm.exclude = MV_SYSTEM;
 127                 }
 128               else if (lex_match_id (lexer, "EXCLUDE"))
 129                 {
 130                   glm.exclude = MV_ANY;
 131                 }
 132               else
 133                 {
 134                   lex_error (lexer, NULL);
 135                   goto error;
 136                 }
 137             }
 138         }
 139       else if (lex_match_id (lexer, "INTERCEPT"))
 140         {
 141           lex_match (lexer, T_EQUALS);
 142           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 143             {
 144               if (lex_match_id (lexer, "INCLUDE"))
 145                 {
 146                   glm.intercept = true;
 147                 }
 148               else if (lex_match_id (lexer, "EXCLUDE"))
 149                 {
 150                   glm.intercept = false;
 151                 }
 152               else
 153                 {
 154                   lex_error (lexer, NULL);
 155                   goto error;
 156                 }
 157             }
 158         }
 159       else if (lex_match_id (lexer, "DESIGN"))
 160         {
 161           size_t n_des;
 162           const struct variable **des;
 163           lex_match (lexer, T_EQUALS);
 164
 165           parse_const_var_set_vars (lexer, factors, &des, &n_des, 0);
 166         }
 167       else
 168         {
 169           lex_error (lexer, NULL);
 170           goto error;
 171         }
 172     }
 173
 174
 175   {
 176     struct casegrouper *grouper;
 177     struct casereader *group;
 178     bool ok;
 179
 180     grouper = casegrouper_create_splits (proc_open (ds), dict);
 181     while (casegrouper_get_next_group (grouper, &group))
 182       run_glm (&glm, group, ds);
 183     ok = casegrouper_destroy (grouper);
 184     ok = proc_commit (ds) && ok;
 185   }
 186
 187   return CMD_SUCCESS;
 188
 189  error:
 190   return CMD_FAILURE;
 191 }
 192
 193 static void get_ssq (struct covariance *, gsl_vector *, struct glm_spec *);
 194
 195 static bool
 196 not_dropped (size_t j, size_t * dropped, size_t n_dropped)
 197 {
 198   size_t i;
 199
 200   for (i = 0; i < n_dropped; i++)
 201     {
 202       if (j == dropped [i])
 203         return false;
 204     }
 205   return true;
 206 }
 207
 208 static void
 209 get_ssq (struct covariance * cov, gsl_vector * ssq, struct glm_spec * cmd)
 210 {
 211   const struct variable **vars;
 212   gsl_matrix * small_cov = NULL;
 213   gsl_matrix * cm = covariance_calculate_unnormalized (cov);
 214   size_t i;
 215   size_t j;
 216   size_t k;
 217   size_t n;
 218   size_t m;
 219   size_t * dropped;
 220   size_t n_dropped;
 221
 222   dropped = xcalloc (covariance_dim (cov), sizeof (*dropped));
 223   vars = xcalloc (covariance_dim (cov), sizeof (*vars));
 224   covariance_get_var_indices (cov, vars);
 225
 226   for (k = 0; k < cmd->n_factor_vars; k++)
 227     {
 228       n_dropped = 0;
 229       for (i = 1; i < covariance_dim (cov); i++)
 230         {
 231           if (vars [i] == cmd->factor_vars [k])
 232             {
 233               dropped [n_dropped++] = i;
 234             }
 235         }
 236       small_cov = gsl_matrix_alloc (cm->size1 - n_dropped, cm->size2 - n_dropped);
 237       gsl_matrix_set (small_cov, 0, 0, gsl_matrix_get (cm, 0, 0));
 238       n = 0;
 239       m = 0;
 240       for (i = 0; i < cm->size1; i++)
 241         {
 242           if (not_dropped (i, dropped, n_dropped))
 243             {
 244               m = 0;
 245               for (j = 0; j < cm->size2; j++)
 246                 {
 247                   if (not_dropped (j, dropped, n_dropped))
 248                     {
 249                       gsl_matrix_set (small_cov, n, m, gsl_matrix_get (cm, i, j));
 250                       m++;
 251                     }
 252                 }
 253               n++;
 254             }
 255         }
 256       reg_sweep (small_cov, 0);
 257       gsl_vector_set (ssq, k + 1,
 258                       gsl_matrix_get (small_cov, 0, 0)
 259                       - gsl_vector_get (ssq, 0));
 260       gsl_matrix_free (small_cov);
 261     }
 262
 263   free (dropped);
 264   free (vars);
 265   gsl_matrix_free (cm);
 266
 267 }
 268
 269 static  void dump_matrix (const gsl_matrix *m);
 270
 271 static void
 272 run_glm (struct glm_spec *cmd, struct casereader *input, const struct dataset *ds)
 273 {
 274   int v;
 275   struct taint *taint;
 276   struct dictionary *dict = dataset_dict (ds);
 277   struct casereader *reader;
 278   struct ccase *c;
 279
 280   struct glm_workspace ws;
 281
 282   struct categoricals *cats = categoricals_create (cmd->factor_vars, cmd->n_factor_vars,
 283                                                    cmd->wv, cmd->exclude,
 284                                                    NULL, NULL,
 285                                                    NULL, NULL);
 286
 287   struct covariance *cov = covariance_2pass_create (cmd->n_dep_vars, cmd->dep_vars,
 288                                                cats,
 289                                                cmd->wv, cmd->exclude);
 290
 291
 292   c = casereader_peek (input, 0);
 293   if (c == NULL)
 294     {
 295       casereader_destroy (input);
 296       return;
 297     }
 298   output_split_file_values (ds, c);
 299   case_unref (c);
 300
 301   taint = taint_clone (casereader_get_taint (input));
 302
 303   ws.totals = moments_create (MOMENT_VARIANCE);
 304
 305   bool warn_bad_weight = true;
 306   for (reader = casereader_clone (input);
 307        (c = casereader_read (reader)) != NULL; case_unref (c))
 308     {
 309       double weight = dict_get_case_weight (dict, c, &warn_bad_weight);
 310
 311       for ( v = 0; v < cmd->n_dep_vars; ++v)
 312         moments_pass_one (ws.totals, case_data (c, cmd->dep_vars[v])->f, weight);
 313
 314       covariance_accumulate_pass1 (cov, c);
 315     }
 316   casereader_destroy (reader);
 317
 318   categoricals_done (cats);
 319
 320   for (reader = casereader_clone (input);
 321        (c = casereader_read (reader)) != NULL; case_unref (c))
 322     {
 323       double weight = dict_get_case_weight (dict, c, &warn_bad_weight);
 324
 325       for ( v = 0; v < cmd->n_dep_vars; ++v)
 326         moments_pass_two (ws.totals, case_data (c, cmd->dep_vars[v])->f, weight);
 327
 328       covariance_accumulate_pass2 (cov, c);
 329     }
 330   casereader_destroy (reader);
 331
 332   {
 333     gsl_matrix *cm = covariance_calculate_unnormalized (cov);
 334
 335     dump_matrix (cm);
 336
 337     ws.total_ssq = gsl_matrix_get (cm, 0, 0);
 338
 339     reg_sweep (cm, 0);
 340
 341     /*
 342       Store the overall SSE.
 343      */
 344     cmd->ssq = gsl_vector_alloc (cm->size1);
 345     gsl_vector_set (cmd->ssq, 0, gsl_matrix_get (cm, 0, 0));
 346     get_ssq (cov, cmd->ssq, cmd);
 347
 348     gsl_vector_free (cmd->ssq);
 349     dump_matrix (cm);
 350
 351     gsl_matrix_free (cm);
 352   }
 353
 354   if (!taint_has_tainted_successor (taint))
 355     output_glm (cmd, &ws);
 356
 357   taint_destroy (taint);
 358 }
 359
 360 static void
 361 output_glm (struct glm_spec *cmd, const struct glm_workspace *ws)
 362 {
 363   const struct fmt_spec *wfmt = cmd->wv ? var_get_print_format (cmd->wv) : &F_8_0;
 364
 365   int f;
 366   int r;
 367   const int heading_columns = 1;
 368   const int heading_rows = 1;
 369   struct tab_table *t ;
 370
 371   const int nc = 6;
 372   int nr = heading_rows + 4 + cmd->n_factor_vars;
 373   if (cmd->intercept)
 374     nr++;
 375
 376   t = tab_create (nc, nr);
 377   tab_title (t, _("Tests of Between-Subjects Effects"));
 378
 379   tab_headers (t, heading_columns, 0, heading_rows, 0);
 380
 381   tab_box (t,
 382            TAL_2, TAL_2,
 383            -1, TAL_1,
 384            0, 0,
 385            nc - 1, nr - 1);
 386
 387   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
 388   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
 389
 390   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
 391
 392   /* TRANSLATORS: The parameter is a roman numeral */
 393   tab_text_format (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Type %s Sum of Squares"), "III");
 394   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df"));
 395   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 396   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("F"));
 397   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
 398
 399   r = heading_rows;
 400   tab_text (t, 0, r++, TAB_LEFT | TAT_TITLE, _("Corrected Model"));
 401
 402   double intercept, n_total;
 403   if (cmd->intercept)
 404     {
 405       double mean;
 406       moments_calculate (ws->totals, &n_total, &mean, NULL, NULL, NULL);
 407       intercept = pow2 (mean * n_total) / n_total;
 408
 409       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Intercept"));
 410       tab_double (t, 1, r, 0, intercept, NULL);
 411       tab_double (t, 2, r, 0, 1.00, wfmt);
 412
 413       tab_double (t, 3, r, 0, intercept / 1.0 , NULL);
 414       r++;
 415     }
 416
 417   for (f = 0; f < cmd->n_factor_vars; ++f)
 418     {
 419       tab_text (t, 0, r++, TAB_LEFT | TAT_TITLE,
 420                 var_to_string (cmd->factor_vars[f]));
 421     }
 422
 423   tab_text (t, 0, r++, TAB_LEFT | TAT_TITLE, _("Error"));
 424
 425   if (cmd->intercept)
 426     {
 427       double ssq = intercept + ws->total_ssq;
 428       double mse = ssq / n_total;
 429       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Total"));
 430       tab_double (t, 1, r, 0, ssq, NULL);
 431       tab_double (t, 2, r, 0, n_total, wfmt);
 432
 433       r++;
 434     }
 435
 436   tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Total"));
 437
 438   tab_double (t, 1, r, 0, ws->total_ssq, NULL);
 439   tab_double (t, 2, r, 0, n_total - 1.0, wfmt);
 440
 441   tab_submit (t);
 442 }
 443
 444 static
 445 void dump_matrix (const gsl_matrix *m)
 446 {
 447   size_t i, j;
 448   for (i = 0; i < m->size1; ++i)
 449     {
 450       for (j = 0; j < m->size2; ++j)
 451         {
 452           double x = gsl_matrix_get (m, i, j);
 453           printf ("%.3f ", x);
 454         }
 455       printf ("\n");
 456     }
 457   printf ("\n");
 458 }