pintos-os.org Git - pspp/blob - 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   /* In the current implementation, design_vars will
  57      normally be the same as factor_vars.
  58      This will change once interactions, nested variables
  59      and repeated measures become involved.
  60   */
  61   size_t n_design_vars;
  62   const struct variable **design_vars;
  63
  64   enum mv_class exclude;
  65
  66   /* The weight variable */
  67   const struct variable *wv;
  68
  69   const struct dictionary *dict;
  70
  71   bool intercept;
  72
  73   double alpha;
  74 };
  75
  76 struct glm_workspace
  77 {
  78   double total_ssq;
  79   struct moments *totals;
  80
  81   struct categoricals *cats;
  82
  83   /*
  84      Sums of squares due to different variables. Element 0 is the SSE
  85      for the entire model. For i > 0, element i is the SS due to
  86      variable i.
  87    */
  88   gsl_vector *ssq;
  89 };
  90
  91 static void output_glm (const struct glm_spec *,
  92                         const struct glm_workspace *ws);
  93 static void run_glm (struct glm_spec *cmd, struct casereader *input,
  94                      const struct dataset *ds);
  95
  96
  97 static bool parse_design_spec (struct lexer *lexer, struct glm_spec *glm);
  98
  99
 100 int
 101 cmd_glm (struct lexer *lexer, struct dataset *ds)
 102 {
 103   struct const_var_set *factors = NULL;
 104   struct glm_spec glm;
 105   bool design = false;
 106   glm.dict = dataset_dict (ds);
 107   glm.n_dep_vars = 0;
 108   glm.n_factor_vars = 0;
 109   glm.n_design_vars = 0;
 110   glm.dep_vars = NULL;
 111   glm.factor_vars = NULL;
 112   glm.design_vars = NULL;
 113   glm.exclude = MV_ANY;
 114   glm.intercept = true;
 115   glm.wv = dict_get_weight (glm.dict);
 116   glm.alpha = 0.05;
 117
 118   if (!parse_variables_const (lexer, glm.dict,
 119                               &glm.dep_vars, &glm.n_dep_vars,
 120                               PV_NO_DUPLICATE | PV_NUMERIC))
 121     goto error;
 122
 123   lex_force_match (lexer, T_BY);
 124
 125   if (!parse_variables_const (lexer, glm.dict,
 126                               &glm.factor_vars, &glm.n_factor_vars,
 127                               PV_NO_DUPLICATE | PV_NUMERIC))
 128     goto error;
 129
 130   if (glm.n_dep_vars > 1)
 131     {
 132       msg (ME, _("Multivariate analysis is not yet implemented"));
 133       return CMD_FAILURE;
 134     }
 135
 136   factors =
 137     const_var_set_create_from_array (glm.factor_vars, glm.n_factor_vars);
 138
 139   while (lex_token (lexer) != T_ENDCMD)
 140     {
 141       lex_match (lexer, T_SLASH);
 142
 143       if (lex_match_id (lexer, "MISSING"))
 144         {
 145           lex_match (lexer, T_EQUALS);
 146           while (lex_token (lexer) != T_ENDCMD
 147                  && lex_token (lexer) != T_SLASH)
 148             {
 149               if (lex_match_id (lexer, "INCLUDE"))
 150                 {
 151                   glm.exclude = MV_SYSTEM;
 152                 }
 153               else if (lex_match_id (lexer, "EXCLUDE"))
 154                 {
 155                   glm.exclude = MV_ANY;
 156                 }
 157               else
 158                 {
 159                   lex_error (lexer, NULL);
 160                   goto error;
 161                 }
 162             }
 163         }
 164       else if (lex_match_id (lexer, "INTERCEPT"))
 165         {
 166           lex_match (lexer, T_EQUALS);
 167           while (lex_token (lexer) != T_ENDCMD
 168                  && lex_token (lexer) != T_SLASH)
 169             {
 170               if (lex_match_id (lexer, "INCLUDE"))
 171                 {
 172                   glm.intercept = true;
 173                 }
 174               else if (lex_match_id (lexer, "EXCLUDE"))
 175                 {
 176                   glm.intercept = false;
 177                 }
 178               else
 179                 {
 180                   lex_error (lexer, NULL);
 181                   goto error;
 182                 }
 183             }
 184         }
 185       else if (lex_match_id (lexer, "CRITERIA"))
 186         {
 187           lex_match (lexer, T_EQUALS);
 188           if (lex_match_id (lexer, "ALPHA"))
 189             {
 190               if (lex_force_match (lexer, T_LPAREN))
 191                 {
 192                   if (! lex_force_num (lexer))
 193                     {
 194                       lex_error (lexer, NULL);
 195                       goto error;
 196                     }
 197
 198                   glm.alpha = lex_number (lexer);
 199                   lex_get (lexer);
 200                   if ( ! lex_force_match (lexer, T_RPAREN))
 201                     {
 202                       lex_error (lexer, NULL);
 203                       goto error;
 204                     }
 205                 }
 206             }
 207           else
 208             {
 209               lex_error (lexer, NULL);
 210               goto error;
 211             }
 212         }
 213       else if (lex_match_id (lexer, "METHOD"))
 214         {
 215           lex_match (lexer, T_EQUALS);
 216           if ( !lex_force_match_id (lexer, "SSTYPE"))
 217             {
 218               lex_error (lexer, NULL);
 219               goto error;
 220             }
 221
 222           if ( ! lex_force_match (lexer, T_LPAREN))
 223             {
 224               lex_error (lexer, NULL);
 225               goto error;
 226             }
 227
 228           if ( ! lex_force_int (lexer))
 229             {
 230               lex_error (lexer, NULL);
 231               goto error;
 232             }
 233
 234           if (3 != lex_integer (lexer))
 235             {
 236               msg (ME, _("Only type 3 sum of squares are currently implemented"));
 237               goto error;
 238             }
 239
 240           lex_get (lexer);
 241
 242           if ( ! lex_force_match (lexer, T_RPAREN))
 243             {
 244               lex_error (lexer, NULL);
 245               goto error;
 246             }
 247         }
 248       else if (lex_match_id (lexer, "DESIGN"))
 249         {
 250           lex_match (lexer, T_EQUALS);
 251
 252           if (! parse_design_spec (lexer, &glm))
 253             goto error;
 254
 255           if ( glm.n_design_vars == 0)
 256             {
 257               msg (ME, _("One or more design  variables must be given"));
 258               goto error;
 259             }
 260
 261           design = true;
 262         }
 263       else
 264         {
 265           lex_error (lexer, NULL);
 266           goto error;
 267         }
 268     }
 269
 270   if ( ! design )
 271     {
 272       lex_error (lexer, _("/DESIGN is mandatory in GLM"));
 273       goto error;
 274     }
 275
 276   {
 277     struct casegrouper *grouper;
 278     struct casereader *group;
 279     bool ok;
 280
 281     grouper = casegrouper_create_splits (proc_open (ds), glm.dict);
 282     while (casegrouper_get_next_group (grouper, &group))
 283       run_glm (&glm, group, ds);
 284     ok = casegrouper_destroy (grouper);
 285     ok = proc_commit (ds) && ok;
 286   }
 287
 288   const_var_set_destroy (factors);
 289   free (glm.factor_vars);
 290   free (glm.design_vars);
 291   free (glm.dep_vars);
 292
 293
 294   return CMD_SUCCESS;
 295
 296 error:
 297
 298   const_var_set_destroy (factors);
 299   free (glm.factor_vars);
 300   free (glm.design_vars);
 301   free (glm.dep_vars);
 302
 303   return CMD_FAILURE;
 304 }
 305
 306 static void get_ssq (struct covariance *, gsl_vector *,
 307                      const struct glm_spec *);
 308
 309 static bool
 310 not_dropped (size_t j, size_t * dropped, size_t n_dropped)
 311 {
 312   size_t i;
 313
 314   for (i = 0; i < n_dropped; i++)
 315     {
 316       if (j == dropped[i])
 317         return false;
 318     }
 319   return true;
 320 }
 321
 322 static void
 323 get_ssq (struct covariance *cov, gsl_vector * ssq, const struct glm_spec *cmd)
 324 {
 325   const struct variable **vars;
 326   gsl_matrix *small_cov = NULL;
 327   gsl_matrix *cm = covariance_calculate_unnormalized (cov);
 328   size_t i;
 329   size_t j;
 330   size_t k;
 331   size_t n;
 332   size_t m;
 333   size_t *dropped;
 334   size_t n_dropped;
 335
 336   dropped = xcalloc (covariance_dim (cov), sizeof (*dropped));
 337   vars = xcalloc (covariance_dim (cov), sizeof (*vars));
 338   covariance_get_var_indices (cov, vars);
 339
 340   for (k = 0; k < cmd->n_design_vars; k++)
 341     {
 342       n_dropped = 0;
 343       for (i = 1; i < covariance_dim (cov); i++)
 344         {
 345           if (vars[i] == cmd->design_vars[k])
 346             {
 347               dropped[n_dropped++] = i;
 348             }
 349         }
 350       small_cov =
 351         gsl_matrix_alloc (cm->size1 - n_dropped, cm->size2 - n_dropped);
 352       gsl_matrix_set (small_cov, 0, 0, gsl_matrix_get (cm, 0, 0));
 353       n = 0;
 354       m = 0;
 355       for (i = 0; i < cm->size1; i++)
 356         {
 357           if (not_dropped (i, dropped, n_dropped))
 358             {
 359               m = 0;
 360               for (j = 0; j < cm->size2; j++)
 361                 {
 362                   if (not_dropped (j, dropped, n_dropped))
 363                     {
 364                       gsl_matrix_set (small_cov, n, m,
 365                                       gsl_matrix_get (cm, i, j));
 366                       m++;
 367                     }
 368                 }
 369               n++;
 370             }
 371         }
 372       reg_sweep (small_cov, 0);
 373       gsl_vector_set (ssq, k + 1,
 374                       gsl_matrix_get (small_cov, 0, 0)
 375                       - gsl_vector_get (ssq, 0));
 376       gsl_matrix_free (small_cov);
 377     }
 378
 379   free (dropped);
 380   free (vars);
 381   gsl_matrix_free (cm);
 382 }
 383
 384 //static  void dump_matrix (const gsl_matrix *m);
 385
 386 static void
 387 run_glm (struct glm_spec *cmd, struct casereader *input,
 388          const struct dataset *ds)
 389 {
 390   bool warn_bad_weight = true;
 391   int v;
 392   struct taint *taint;
 393   struct dictionary *dict = dataset_dict (ds);
 394   struct casereader *reader;
 395   struct ccase *c;
 396
 397   struct glm_workspace ws;
 398   struct covariance *cov;
 399   ws.cats = categoricals_create (cmd->design_vars, cmd->n_design_vars,
 400                                  cmd->wv, cmd->exclude,
 401                                  NULL, NULL, NULL, NULL);
 402
 403   cov = covariance_2pass_create (cmd->n_dep_vars, cmd->dep_vars,
 404                                  ws.cats, cmd->wv, cmd->exclude);
 405
 406
 407   c = casereader_peek (input, 0);
 408   if (c == NULL)
 409     {
 410       casereader_destroy (input);
 411       return;
 412     }
 413   output_split_file_values (ds, c);
 414   case_unref (c);
 415
 416   taint = taint_clone (casereader_get_taint (input));
 417
 418   ws.totals = moments_create (MOMENT_VARIANCE);
 419
 420   for (reader = casereader_clone (input);
 421        (c = casereader_read (reader)) != NULL; case_unref (c))
 422     {
 423       double weight = dict_get_case_weight (dict, c, &warn_bad_weight);
 424
 425       for (v = 0; v < cmd->n_dep_vars; ++v)
 426         moments_pass_one (ws.totals, case_data (c, cmd->dep_vars[v])->f,
 427                           weight);
 428
 429       covariance_accumulate_pass1 (cov, c);
 430     }
 431   casereader_destroy (reader);
 432
 433   categoricals_done (ws.cats);
 434
 435   for (reader = input;
 436        (c = casereader_read (reader)) != NULL; case_unref (c))
 437     {
 438       double weight = dict_get_case_weight (dict, c, &warn_bad_weight);
 439
 440       for (v = 0; v < cmd->n_dep_vars; ++v)
 441         moments_pass_two (ws.totals, case_data (c, cmd->dep_vars[v])->f,
 442                           weight);
 443
 444       covariance_accumulate_pass2 (cov, c);
 445     }
 446   casereader_destroy (reader);
 447
 448   {
 449     gsl_matrix *cm = covariance_calculate_unnormalized (cov);
 450
 451     //    dump_matrix (cm);
 452
 453     ws.total_ssq = gsl_matrix_get (cm, 0, 0);
 454
 455     reg_sweep (cm, 0);
 456
 457     /*
 458        Store the overall SSE.
 459      */
 460     ws.ssq = gsl_vector_alloc (cm->size1);
 461     gsl_vector_set (ws.ssq, 0, gsl_matrix_get (cm, 0, 0));
 462     get_ssq (cov, ws.ssq, cmd);
 463     //    dump_matrix (cm);
 464
 465     gsl_matrix_free (cm);
 466   }
 467
 468   if (!taint_has_tainted_successor (taint))
 469     output_glm (cmd, &ws);
 470
 471   gsl_vector_free (ws.ssq);
 472
 473   covariance_destroy (cov);
 474   moments_destroy (ws.totals);
 475
 476   taint_destroy (taint);
 477 }
 478
 479 static void
 480 output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
 481 {
 482   const struct fmt_spec *wfmt =
 483     cmd->wv ? var_get_print_format (cmd->wv) : &F_8_0;
 484
 485   double n_total, mean;
 486   double df_corr = 0.0;
 487   double mse = 0;
 488
 489   int f;
 490   int r;
 491   const int heading_columns = 1;
 492   const int heading_rows = 1;
 493   struct tab_table *t;
 494
 495   const int nc = 6;
 496   int nr = heading_rows + 4 + cmd->n_design_vars;
 497   if (cmd->intercept)
 498     nr++;
 499
 500   t = tab_create (nc, nr);
 501   tab_title (t, _("Tests of Between-Subjects Effects"));
 502
 503   tab_headers (t, heading_columns, 0, heading_rows, 0);
 504
 505   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
 506
 507   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
 508   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
 509
 510   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
 511
 512   /* TRANSLATORS: The parameter is a roman numeral */
 513   tab_text_format (t, 1, 0, TAB_CENTER | TAT_TITLE,
 514                    _("Type %s Sum of Squares"), "III");
 515   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df"));
 516   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 517   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("F"));
 518   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
 519
 520   moments_calculate (ws->totals, &n_total, &mean, NULL, NULL, NULL);
 521
 522   if (cmd->intercept)
 523     df_corr += 1.0;
 524
 525   for (f = 0; f < cmd->n_design_vars; ++f)
 526     df_corr += categoricals_n_count (ws->cats, f) - 1.0;
 527
 528   mse = gsl_vector_get (ws->ssq, 0) / (n_total - df_corr);
 529
 530   r = heading_rows;
 531   tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Model"));
 532
 533   r++;
 534
 535   if (cmd->intercept)
 536     {
 537       const double intercept = pow2 (mean * n_total) / n_total;
 538       const double df = 1.0;
 539       const double F = intercept / df / mse;
 540       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Intercept"));
 541       tab_double (t, 1, r, 0, intercept, NULL);
 542       tab_double (t, 2, r, 0, 1.00, wfmt);
 543       tab_double (t, 3, r, 0, intercept / df, NULL);
 544       tab_double (t, 4, r, 0, F, NULL);
 545       tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr),
 546                   NULL);
 547       r++;
 548     }
 549
 550   for (f = 0; f < cmd->n_design_vars; ++f)
 551     {
 552       const double df = categoricals_n_count (ws->cats, f) - 1.0;
 553       const double ssq = gsl_vector_get (ws->ssq, f + 1);
 554       const double F = ssq / df / mse;
 555       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
 556                 var_to_string (cmd->design_vars[f]));
 557
 558       tab_double (t, 1, r, 0, ssq, NULL);
 559       tab_double (t, 2, r, 0, df, wfmt);
 560       tab_double (t, 3, r, 0, ssq / df, NULL);
 561       tab_double (t, 4, r, 0, F, NULL);
 562
 563       tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr),
 564                   NULL);
 565
 566
 567       r++;
 568     }
 569
 570   {
 571     /* Corrected Model */
 572     const double df = df_corr - 1.0;
 573     const double ssq = ws->total_ssq - gsl_vector_get (ws->ssq, 0);
 574     const double F = ssq / df / mse;
 575     tab_double (t, 1, heading_rows, 0, ssq, NULL);
 576     tab_double (t, 2, heading_rows, 0, df, wfmt);
 577     tab_double (t, 3, heading_rows, 0, ssq / df, NULL);
 578     tab_double (t, 4, heading_rows, 0, F, NULL);
 579
 580     tab_double (t, 5, heading_rows, 0,
 581                 gsl_cdf_fdist_Q (F, df, n_total - df_corr), NULL);
 582   }
 583
 584   {
 585     const double df = n_total - df_corr;
 586     const double ssq = gsl_vector_get (ws->ssq, 0);
 587     const double mse = ssq / df;
 588     tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Error"));
 589     tab_double (t, 1, r, 0, ssq, NULL);
 590     tab_double (t, 2, r, 0, df, wfmt);
 591     tab_double (t, 3, r++, 0, mse, NULL);
 592   }
 593
 594   if (cmd->intercept)
 595     {
 596       const double intercept = pow2 (mean * n_total) / n_total;
 597       const double ssq = intercept + ws->total_ssq;
 598
 599       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Total"));
 600       tab_double (t, 1, r, 0, ssq, NULL);
 601       tab_double (t, 2, r, 0, n_total, wfmt);
 602
 603       r++;
 604     }
 605
 606   tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Total"));
 607
 608
 609   tab_double (t, 1, r, 0, ws->total_ssq, NULL);
 610   tab_double (t, 2, r, 0, n_total - 1.0, wfmt);
 611
 612   tab_submit (t);
 613 }
 614
 615 #if 0
 616 static void
 617 dump_matrix (const gsl_matrix * m)
 618 {
 619   size_t i, j;
 620   for (i = 0; i < m->size1; ++i)
 621     {
 622       for (j = 0; j < m->size2; ++j)
 623         {
 624           double x = gsl_matrix_get (m, i, j);
 625           printf ("%.3f ", x);
 626         }
 627       printf ("\n");
 628     }
 629   printf ("\n");
 630 }
 631 #endif
 632
 633
 634 \f
 635
 636 /* Match a variable.
 637    If the match succeeds, the variable will be placed in VAR.
 638    Returns true if successful */
 639 static bool
 640 lex_match_variable (struct lexer *lexer, const struct glm_spec *glm, const struct variable **var)
 641 {
 642   if (lex_token (lexer) !=  T_ID)
 643     return false;
 644
 645   *var = parse_variable_const  (lexer, glm->dict);
 646
 647   if ( *var == NULL)
 648     return false;
 649   return true;
 650 }
 651
 652 /* An interaction is a variable followed by {*, BY} followed by an interaction */
 653 static bool
 654 parse_design_interaction (struct lexer *lexer, struct glm_spec *glm)
 655 {
 656   const struct variable *v = NULL;
 657   if (! lex_match_variable (lexer, glm, &v))
 658     return false;
 659
 660   if ( lex_match (lexer, T_ASTERISK) || lex_match (lexer, T_BY))
 661     {
 662       lex_error (lexer, "Interactions are not yet implemented"); return false;
 663       return parse_design_interaction (lexer, glm);
 664     }
 665
 666   glm->n_design_vars++;
 667   glm->design_vars = xrealloc (glm->design_vars, sizeof (*glm->design_vars) * glm->n_design_vars);
 668   glm->design_vars[glm->n_design_vars - 1] = v;
 669
 670   return true;
 671 }
 672
 673 /* A design term is a varible OR an interaction */
 674 static bool
 675 parse_design_term (struct lexer *lexer, struct glm_spec *glm)
 676 {
 677   const struct variable *v = NULL;
 678   if (parse_design_interaction (lexer, glm))
 679     return true;
 680
 681   /* FIXME: This should accept nexted variables */
 682   if ( lex_match_variable (lexer, glm, &v))
 683     {
 684       return true;
 685     }
 686
 687   return false;
 688 }
 689
 690
 691
 692 /* Parse a complete DESIGN specification.
 693    A design spec is a design term, optionally followed by a comma,
 694    and another design spec.
 695 */
 696 static bool
 697 parse_design_spec (struct lexer *lexer, struct glm_spec *glm)
 698 {
 699   /* Kludge: Return success if end of design spec */
 700   if  (lex_token (lexer) == T_ENDCMD || lex_token (lexer) == T_SLASH)
 701     return true;
 702
 703   if ( ! parse_design_term (lexer, glm))
 704     return false;
 705
 706   lex_match (lexer, T_COMMA);
 707
 708   return parse_design_spec (lexer, glm);
 709 }
 710