pintos-os.org Git - pspp/blob - src/language/stats/regression.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005, 2009, 2010, 2011, 2012, 2013 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 <stdbool.h>
  20
  21 #include <gsl/gsl_cdf.h>
  22 #include <gsl/gsl_matrix.h>
  23
  24 #include <data/dataset.h>
  25
  26 #include "language/command.h"
  27 #include "language/lexer/lexer.h"
  28 #include "language/lexer/value-parser.h"
  29 #include "language/lexer/variable-parser.h"
  30
  31
  32 #include "data/casegrouper.h"
  33 #include "data/casereader.h"
  34 #include "data/dictionary.h"
  35
  36 #include "math/covariance.h"
  37 #include "math/linreg.h"
  38 #include "math/moments.h"
  39
  40 #include "libpspp/message.h"
  41 #include "libpspp/taint.h"
  42
  43 #include "output/tab.h"
  44
  45 #include "gettext.h"
  46 #define _(msgid) gettext (msgid)
  47 #define N_(msgid) msgid
  48
  49
  50 #include <gl/intprops.h>
  51
  52 #define REG_LARGE_DATA 1000
  53
  54 struct regression
  55 {
  56   struct dataset *ds;
  57
  58   const struct variable **vars;
  59   size_t n_vars;
  60
  61   const struct variable **dep_vars;
  62   size_t n_dep_vars;
  63
  64   bool r;
  65   bool coeff;
  66   bool anova;
  67   bool bcov;
  68
  69
  70   bool resid;
  71   bool pred;
  72
  73   linreg **models;
  74 };
  75
  76
  77 static void run_regression (const struct regression *cmd, struct casereader *input);
  78
  79
  80
  81 /*
  82   Transformations for saving predicted values
  83   and residuals, etc.
  84 */
  85 struct reg_trns
  86 {
  87   int n_trns;                   /* Number of transformations. */
  88   int trns_id;                  /* Which trns is this one? */
  89   linreg *c;            /* Linear model for this trns. */
  90 };
  91
  92 /*
  93   Gets the predicted values.
  94 */
  95 static int
  96 regression_trns_pred_proc (void *t_, struct ccase **c,
  97                            casenumber case_idx UNUSED)
  98 {
  99   size_t i;
 100   size_t n_vals;
 101   struct reg_trns *trns = t_;
 102   linreg *model;
 103   union value *output = NULL;
 104   const union value *tmp;
 105   double *vals;
 106   const struct variable **vars = NULL;
 107
 108   assert (trns != NULL);
 109   model = trns->c;
 110   assert (model != NULL);
 111   assert (model->depvar != NULL);
 112   assert (model->pred != NULL);
 113
 114   vars = linreg_get_vars (model);
 115   n_vals = linreg_n_coeffs (model);
 116   vals = xnmalloc (n_vals, sizeof (*vals));
 117   *c = case_unshare (*c);
 118
 119   output = case_data_rw (*c, model->pred);
 120
 121   for (i = 0; i < n_vals; i++)
 122     {
 123       tmp = case_data (*c, vars[i]);
 124       vals[i] = tmp->f;
 125     }
 126   output->f = linreg_predict (model, vals, n_vals);
 127   free (vals);
 128   return TRNS_CONTINUE;
 129 }
 130
 131 /*
 132   Gets the residuals.
 133 */
 134 static int
 135 regression_trns_resid_proc (void *t_, struct ccase **c,
 136                             casenumber case_idx UNUSED)
 137 {
 138   size_t i;
 139   size_t n_vals;
 140   struct reg_trns *trns = t_;
 141   linreg *model;
 142   union value *output = NULL;
 143   const union value *tmp;
 144   double *vals = NULL;
 145   double obs;
 146   const struct variable **vars = NULL;
 147
 148   assert (trns != NULL);
 149   model = trns->c;
 150   assert (model != NULL);
 151   assert (model->depvar != NULL);
 152   assert (model->resid != NULL);
 153
 154   vars = linreg_get_vars (model);
 155   n_vals = linreg_n_coeffs (model);
 156
 157   vals = xnmalloc (n_vals, sizeof (*vals));
 158   *c = case_unshare (*c);
 159   output = case_data_rw (*c, model->resid);
 160   assert (output != NULL);
 161
 162   for (i = 0; i < n_vals; i++)
 163     {
 164       tmp = case_data (*c, vars[i]);
 165       vals[i] = tmp->f;
 166     }
 167   tmp = case_data (*c, model->depvar);
 168   obs = tmp->f;
 169   output->f = linreg_residual (model, obs, vals, n_vals);
 170   free (vals);
 171
 172   return TRNS_CONTINUE;
 173 }
 174
 175
 176 static char *
 177 reg_get_name (const struct dictionary *dict, const char *prefix)
 178 {
 179   char *name;
 180   int i;
 181
 182   /* XXX handle too-long prefixes */
 183   name = xmalloc (strlen (prefix) + INT_BUFSIZE_BOUND (i) + 1);
 184   for (i = 1; ; i++)
 185     {
 186       sprintf (name, "%s%d", prefix, i);
 187       if (dict_lookup_var (dict, name) == NULL)
 188         return name;
 189     }
 190 }
 191
 192 /*
 193   Free the transformation. Free its linear model if this
 194   transformation is the last one.
 195 */
 196 static bool
 197 regression_trns_free (void *t_)
 198 {
 199   struct reg_trns *t = t_;
 200
 201   if (t->trns_id == t->n_trns)
 202     {
 203       linreg_unref (t->c);
 204     }
 205   free (t);
 206
 207   return true;
 208 }
 209
 210 static void
 211 reg_save_var (struct dataset *ds, const char *prefix, trns_proc_func * f,
 212               linreg * c, struct variable **v, int n_trns)
 213 {
 214   struct dictionary *dict = dataset_dict (ds);
 215   static int trns_index = 1;
 216   char *name;
 217   struct variable *new_var;
 218   struct reg_trns *t = NULL;
 219
 220   t = xmalloc (sizeof (*t));
 221   t->trns_id = trns_index;
 222   t->n_trns = n_trns;
 223   t->c = c;
 224
 225   name = reg_get_name (dict, prefix);
 226   new_var = dict_create_var_assert (dict, name, 0);
 227   free (name);
 228
 229   *v = new_var;
 230   add_transformation (ds, f, regression_trns_free, t);
 231   trns_index++;
 232 }
 233
 234 static void
 235 subcommand_save (const struct regression *cmd)
 236 {
 237   linreg **lc;
 238   int n_trns = 0;
 239
 240   if ( cmd->resid ) n_trns++;
 241   if ( cmd->pred ) n_trns++;
 242
 243   n_trns *= cmd->n_dep_vars;
 244
 245   for (lc = cmd->models; lc < cmd->models + cmd->n_dep_vars; lc++)
 246     {
 247       if (*lc != NULL)
 248         {
 249           if ((*lc)->depvar != NULL)
 250             {
 251               (*lc)->refcnt++;
 252               if (cmd->resid)
 253                 {
 254                   reg_save_var (cmd->ds, "RES", regression_trns_resid_proc, *lc,
 255                                 &(*lc)->resid, n_trns);
 256                 }
 257               if (cmd->pred)
 258                 {
 259                   reg_save_var (cmd->ds, "PRED", regression_trns_pred_proc, *lc,
 260                                 &(*lc)->pred, n_trns);
 261                 }
 262             }
 263         }
 264     }
 265 }
 266
 267 int
 268 cmd_regression (struct lexer *lexer, struct dataset *ds)
 269 {
 270   int k;
 271   struct regression regression;
 272   const struct dictionary *dict = dataset_dict (ds);
 273   bool save;
 274
 275   memset (&regression, 0, sizeof (struct regression));
 276
 277   regression.anova = true;
 278   regression.coeff = true;
 279   regression.r = true;
 280
 281   regression.pred = false;
 282   regression.resid = false;
 283
 284   regression.ds = ds;
 285
 286   /* Accept an optional, completely pointless "/VARIABLES=" */
 287   lex_match (lexer, T_SLASH);
 288   if (lex_match_id  (lexer, "VARIABLES"))
 289     {
 290       if (! lex_force_match (lexer, T_EQUALS) )
 291         goto error;
 292     }
 293
 294   if (!parse_variables_const (lexer, dict,
 295                               &regression.vars, &regression.n_vars,
 296                               PV_NO_DUPLICATE | PV_NUMERIC))
 297     goto error;
 298
 299
 300   while (lex_token (lexer) != T_ENDCMD)
 301     {
 302       lex_match (lexer, T_SLASH);
 303
 304       if (lex_match_id  (lexer, "DEPENDENT"))
 305         {
 306           if (! lex_force_match (lexer, T_EQUALS) )
 307             goto error;
 308
 309           if (!parse_variables_const (lexer, dict,
 310                                       &regression.dep_vars, &regression.n_dep_vars,
 311                                       PV_NO_DUPLICATE | PV_NUMERIC))
 312             goto error;
 313         }
 314       else if (lex_match_id (lexer, "METHOD"))
 315         {
 316           lex_match (lexer, T_EQUALS);
 317
 318           if (!lex_force_match_id (lexer, "ENTER"))
 319             {
 320               goto error;
 321             }
 322         }
 323       else if (lex_match_id (lexer, "STATISTICS"))
 324         {
 325           lex_match (lexer, T_EQUALS);
 326
 327           while (lex_token (lexer) != T_ENDCMD
 328                  && lex_token (lexer) != T_SLASH)
 329             {
 330               if (lex_match (lexer, T_ALL))
 331                 {
 332                 }
 333               else if (lex_match_id (lexer, "DEFAULTS"))
 334                 {
 335                 }
 336               else if (lex_match_id (lexer, "R"))
 337                 {
 338                 }
 339               else if (lex_match_id (lexer, "COEFF"))
 340                 {
 341                 }
 342               else if (lex_match_id (lexer, "ANOVA"))
 343                 {
 344                 }
 345               else if (lex_match_id (lexer, "BCOV"))
 346                 {
 347                 }
 348               else
 349                 {
 350                   lex_error (lexer, NULL);
 351                   goto error;
 352                 }
 353             }
 354         }
 355       else if (lex_match_id (lexer, "SAVE"))
 356         {
 357           lex_match (lexer, T_EQUALS);
 358
 359           while (lex_token (lexer) != T_ENDCMD
 360                  && lex_token (lexer) != T_SLASH)
 361             {
 362               if (lex_match_id (lexer, "PRED"))
 363                 {
 364                   regression.pred = true;
 365                 }
 366               else if (lex_match_id (lexer, "RESID"))
 367                 {
 368                   regression.resid = true;
 369                 }
 370               else
 371                 {
 372                   lex_error (lexer, NULL);
 373                   goto error;
 374                 }
 375             }
 376         }
 377       else
 378         {
 379           lex_error (lexer, NULL);
 380           goto error;
 381         }
 382     }
 383
 384   if (!regression.vars)
 385     {
 386       dict_get_vars (dict, &regression.vars, &regression.n_vars, 0);
 387     }
 388
 389
 390   regression.models = xcalloc (regression.n_dep_vars, sizeof *regression.models);
 391
 392   save = regression.pred || regression.resid;
 393   if (save)
 394     {
 395       if (proc_make_temporary_transformations_permanent (ds))
 396         msg (SW, _("REGRESSION with SAVE ignores TEMPORARY.  "
 397                    "Temporary transformations will be made permanent."));
 398     }
 399
 400   {
 401     struct casegrouper *grouper;
 402     struct casereader *group;
 403     bool ok;
 404
 405     grouper = casegrouper_create_splits (proc_open_filtering (ds, !save),
 406                                          dict);
 407     while (casegrouper_get_next_group (grouper, &group))
 408       run_regression (&regression, group);
 409     ok = casegrouper_destroy (grouper);
 410     ok = proc_commit (ds) && ok;
 411   }
 412
 413   if (save)
 414     {
 415       subcommand_save (&regression);
 416     }
 417
 418
 419   for (k = 0; k < regression.n_dep_vars; k++)
 420     linreg_unref (regression.models[k]);
 421   free (regression.models);
 422   free (regression.vars);
 423   free (regression.dep_vars);
 424   return CMD_SUCCESS;
 425
 426  error:
 427   if (regression.models)
 428    {
 429   for (k = 0; k < regression.n_dep_vars; k++)
 430     linreg_unref (regression.models[k]);
 431   free (regression.models);
 432    }
 433   free (regression.vars);
 434   free (regression.dep_vars);
 435   return CMD_FAILURE;
 436 }
 437
 438
 439 static size_t
 440 get_n_all_vars (const struct regression *cmd)
 441 {
 442   size_t result = cmd->n_vars;
 443   size_t i;
 444   size_t j;
 445
 446   result += cmd->n_dep_vars;
 447   for (i = 0; i < cmd->n_dep_vars; i++)
 448     {
 449       for (j = 0; j < cmd->n_vars; j++)
 450         {
 451           if (cmd->vars[j] == cmd->dep_vars[i])
 452             {
 453               result--;
 454             }
 455         }
 456     }
 457   return result;
 458 }
 459
 460 static void
 461 fill_all_vars (const struct variable **vars, const struct regression *cmd)
 462 {
 463   size_t i;
 464   size_t j;
 465   bool absent;
 466
 467   for (i = 0; i < cmd->n_vars; i++)
 468     {
 469       vars[i] = cmd->vars[i];
 470     }
 471   for (i = 0; i < cmd->n_dep_vars; i++)
 472     {
 473       absent = true;
 474       for (j = 0; j < cmd->n_vars; j++)
 475         {
 476           if (cmd->dep_vars[i] == cmd->vars[j])
 477             {
 478               absent = false;
 479               break;
 480             }
 481         }
 482       if (absent)
 483         {
 484           vars[i + cmd->n_vars] = cmd->dep_vars[i];
 485         }
 486     }
 487 }
 488
 489 /*
 490   Is variable k the dependent variable?
 491 */
 492 static bool
 493 is_depvar (const struct regression *cmd, size_t k, const struct variable *v)
 494 {
 495   return v == cmd->vars[k];
 496 }
 497
 498
 499 /* Identify the explanatory variables in v_variables.  Returns
 500    the number of independent variables. */
 501 static int
 502 identify_indep_vars (const struct regression *cmd,
 503                      const struct variable **indep_vars,
 504                      const struct variable *depvar)
 505 {
 506   int n_indep_vars = 0;
 507   int i;
 508
 509   for (i = 0; i < cmd->n_vars; i++)
 510     if (!is_depvar (cmd, i, depvar))
 511       indep_vars[n_indep_vars++] = cmd->vars[i];
 512   if ((n_indep_vars < 1) && is_depvar (cmd, 0, depvar))
 513     {
 514       /*
 515         There is only one independent variable, and it is the same
 516         as the dependent variable. Print a warning and continue.
 517       */
 518       msg (SW,
 519            gettext ("The dependent variable is equal to the independent variable."
 520                     "The least squares line is therefore Y=X."
 521                     "Standard errors and related statistics may be meaningless."));
 522       n_indep_vars = 1;
 523       indep_vars[0] = cmd->vars[0];
 524     }
 525   return n_indep_vars;
 526 }
 527
 528
 529 static double
 530 fill_covariance (gsl_matrix *cov, struct covariance *all_cov,
 531                  const struct variable **vars,
 532                  size_t n_vars, const struct variable *dep_var,
 533                  const struct variable **all_vars, size_t n_all_vars,
 534                  double *means)
 535 {
 536   size_t i;
 537   size_t j;
 538   size_t dep_subscript;
 539   size_t *rows;
 540   const gsl_matrix *ssizes;
 541   const gsl_matrix *mean_matrix;
 542   const gsl_matrix *ssize_matrix;
 543   double result = 0.0;
 544
 545   gsl_matrix *cm = covariance_calculate_unnormalized (all_cov);
 546
 547   if ( cm == NULL)
 548     return 0;
 549
 550   rows = xnmalloc (cov->size1 - 1, sizeof (*rows));
 551
 552   for (i = 0; i < n_all_vars; i++)
 553     {
 554       for (j = 0; j < n_vars; j++)
 555         {
 556           if (vars[j] == all_vars[i])
 557             {
 558               rows[j] = i;
 559             }
 560         }
 561       if (all_vars[i] == dep_var)
 562         {
 563           dep_subscript = i;
 564         }
 565     }
 566   mean_matrix = covariance_moments (all_cov, MOMENT_MEAN);
 567   ssize_matrix = covariance_moments (all_cov, MOMENT_NONE);
 568   for (i = 0; i < cov->size1 - 1; i++)
 569     {
 570       means[i] = gsl_matrix_get (mean_matrix, rows[i], 0)
 571         / gsl_matrix_get (ssize_matrix, rows[i], 0);
 572       for (j = 0; j < cov->size2 - 1; j++)
 573         {
 574           gsl_matrix_set (cov, i, j, gsl_matrix_get (cm, rows[i], rows[j]));
 575           gsl_matrix_set (cov, j, i, gsl_matrix_get (cm, rows[j], rows[i]));
 576         }
 577     }
 578   means[cov->size1 - 1] = gsl_matrix_get (mean_matrix, dep_subscript, 0)
 579     / gsl_matrix_get (ssize_matrix, dep_subscript, 0);
 580   ssizes = covariance_moments (all_cov, MOMENT_NONE);
 581   result = gsl_matrix_get (ssizes, dep_subscript, rows[0]);
 582   for (i = 0; i < cov->size1 - 1; i++)
 583     {
 584       gsl_matrix_set (cov, i, cov->size1 - 1,
 585                       gsl_matrix_get (cm, rows[i], dep_subscript));
 586       gsl_matrix_set (cov, cov->size1 - 1, i,
 587                       gsl_matrix_get (cm, rows[i], dep_subscript));
 588       if (result > gsl_matrix_get (ssizes, rows[i], dep_subscript))
 589         {
 590           result = gsl_matrix_get (ssizes, rows[i], dep_subscript);
 591         }
 592     }
 593   gsl_matrix_set (cov, cov->size1 - 1, cov->size1 - 1,
 594                   gsl_matrix_get (cm, dep_subscript, dep_subscript));
 595   free (rows);
 596   gsl_matrix_free (cm);
 597   return result;
 598 }
 599
 600
 601 /*
 602   STATISTICS subcommand output functions.
 603 */
 604 static void reg_stats_r (linreg *, void *, const struct variable *);
 605 static void reg_stats_coeff (linreg *, void *, const struct variable *);
 606 static void reg_stats_anova (linreg *, void *, const struct variable *);
 607 static void reg_stats_bcov (linreg *, void *, const struct variable *);
 608
 609 static void statistics_keyword_output (void (*)(linreg *, void *, const struct variable *),
 610                                        bool, linreg *, void *, const struct variable *);
 611
 612
 613
 614 static void
 615 subcommand_statistics (const struct regression *cmd , linreg * c, void *aux,
 616                        const struct variable *var)
 617 {
 618   statistics_keyword_output (reg_stats_r, cmd->r, c, aux, var);
 619   statistics_keyword_output (reg_stats_anova, cmd->anova, c, aux, var);
 620   statistics_keyword_output (reg_stats_coeff, cmd->coeff, c, aux, var);
 621   statistics_keyword_output (reg_stats_bcov, cmd->bcov, c, aux, var);
 622 }
 623
 624
 625 static void
 626 run_regression (const struct regression *cmd, struct casereader *input)
 627 {
 628   size_t i;
 629   int n_indep = 0;
 630   int k;
 631   double *means;
 632   struct ccase *c;
 633   struct covariance *cov;
 634   const struct variable **vars;
 635   const struct variable **all_vars;
 636   struct casereader *reader;
 637   size_t n_all_vars;
 638
 639   linreg **models = cmd->models;
 640
 641   n_all_vars = get_n_all_vars (cmd);
 642   all_vars = xnmalloc (n_all_vars, sizeof (*all_vars));
 643   fill_all_vars (all_vars, cmd);
 644   vars = xnmalloc (cmd->n_vars, sizeof (*vars));
 645   means  = xnmalloc (n_all_vars, sizeof (*means));
 646   cov = covariance_1pass_create (n_all_vars, all_vars,
 647                                  dict_get_weight (dataset_dict (cmd->ds)), MV_ANY);
 648
 649   reader = casereader_clone (input);
 650   reader = casereader_create_filter_missing (reader, all_vars, n_all_vars,
 651                                              MV_ANY, NULL, NULL);
 652
 653
 654   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
 655     {
 656       covariance_accumulate (cov, c);
 657     }
 658
 659   for (k = 0; k < cmd->n_dep_vars; k++)
 660     {
 661       double n_data;
 662       const struct variable *dep_var = cmd->dep_vars[k];
 663       gsl_matrix *this_cm;
 664
 665       n_indep = identify_indep_vars (cmd, vars, dep_var);
 666
 667       this_cm = gsl_matrix_alloc (n_indep + 1, n_indep + 1);
 668       n_data = fill_covariance (this_cm, cov, vars, n_indep,
 669                                 dep_var, all_vars, n_all_vars, means);
 670       models[k] = linreg_alloc (dep_var, (const struct variable **) vars,
 671                                 n_data, n_indep);
 672       models[k]->depvar = dep_var;
 673       for (i = 0; i < n_indep; i++)
 674         {
 675           linreg_set_indep_variable_mean (models[k], i, means[i]);
 676         }
 677       linreg_set_depvar_mean (models[k], means[i]);
 678       /*
 679         For large data sets, use QR decomposition.
 680       */
 681       if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
 682         {
 683           models[k]->method = LINREG_QR;
 684         }
 685
 686       if (n_data > 0)
 687         {
 688           /*
 689             Find the least-squares estimates and other statistics.
 690           */
 691           linreg_fit (this_cm, models[k]);
 692
 693           if (!taint_has_tainted_successor (casereader_get_taint (input)))
 694             {
 695               subcommand_statistics (cmd, models[k], this_cm, dep_var);
 696             }
 697         }
 698       else
 699         {
 700           msg (SE,
 701                _("No valid data found. This command was skipped."));
 702           linreg_unref (models[k]);
 703           models[k] = NULL;
 704         }
 705       gsl_matrix_free (this_cm);
 706     }
 707
 708   casereader_destroy (reader);
 709   free (vars);
 710   free (all_vars);
 711   free (means);
 712   casereader_destroy (input);
 713   covariance_destroy (cov);
 714 }
 715
 716
 717 \f
 718
 719
 720 static void
 721 reg_stats_r (linreg *c, void *aux UNUSED, const struct variable *var)
 722 {
 723   struct tab_table *t;
 724   int n_rows = 2;
 725   int n_cols = 5;
 726   double rsq;
 727   double adjrsq;
 728   double std_error;
 729
 730   assert (c != NULL);
 731   rsq = linreg_ssreg (c) / linreg_sst (c);
 732   adjrsq = 1.0 - (1.0 - rsq) * (linreg_n_obs (c) - 1.0) / (linreg_n_obs (c) - linreg_n_coeffs (c));
 733   std_error = sqrt (linreg_mse (c));
 734   t = tab_create (n_cols, n_rows);
 735   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 736   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 737   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 738   tab_vline (t, TAL_0, 1, 0, 0);
 739
 740   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
 741   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
 742   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
 743   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
 744   tab_double (t, 1, 1, TAB_RIGHT, sqrt (rsq), NULL);
 745   tab_double (t, 2, 1, TAB_RIGHT, rsq, NULL);
 746   tab_double (t, 3, 1, TAB_RIGHT, adjrsq, NULL);
 747   tab_double (t, 4, 1, TAB_RIGHT, std_error, NULL);
 748   tab_title (t, _("Model Summary (%s)"), var_to_string (var));
 749   tab_submit (t);
 750 }
 751
 752 /*
 753   Table showing estimated regression coefficients.
 754 */
 755 static void
 756 reg_stats_coeff (linreg * c, void *aux_, const struct variable *var)
 757 {
 758   size_t j;
 759   int n_cols = 7;
 760   int n_rows;
 761   int this_row;
 762   double t_stat;
 763   double pval;
 764   double std_err;
 765   double beta;
 766   const char *label;
 767
 768   const struct variable *v;
 769   struct tab_table *t;
 770   gsl_matrix *cov = aux_;
 771
 772   assert (c != NULL);
 773   n_rows = linreg_n_coeffs (c) + 3;
 774
 775   t = tab_create (n_cols, n_rows);
 776   tab_headers (t, 2, 0, 1, 0);
 777   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 778   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 779   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 780   tab_vline (t, TAL_0, 1, 0, 0);
 781
 782   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
 783   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 784   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
 785   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
 786   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 787   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
 788   tab_double (t, 2, 1, 0, linreg_intercept (c), NULL);
 789   std_err = sqrt (gsl_matrix_get (linreg_cov (c), 0, 0));
 790   tab_double (t, 3, 1, 0, std_err, NULL);
 791   tab_double (t, 4, 1, 0, 0.0, NULL);
 792   t_stat = linreg_intercept (c) / std_err;
 793   tab_double (t, 5, 1, 0, t_stat, NULL);
 794   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), (double) (linreg_n_obs (c) - linreg_n_coeffs (c)));
 795   tab_double (t, 6, 1, 0, pval, NULL);
 796   for (j = 0; j < linreg_n_coeffs (c); j++)
 797     {
 798       struct string tstr;
 799       ds_init_empty (&tstr);
 800       this_row = j + 2;
 801
 802       v = linreg_indep_var (c, j);
 803       label = var_to_string (v);
 804       /* Do not overwrite the variable's name. */
 805       ds_put_cstr (&tstr, label);
 806       tab_text (t, 1, this_row, TAB_CENTER, ds_cstr (&tstr));
 807       /*
 808         Regression coefficients.
 809       */
 810       tab_double (t, 2, this_row, 0, linreg_coeff (c, j), NULL);
 811       /*
 812         Standard error of the coefficients.
 813       */
 814       std_err = sqrt (gsl_matrix_get (linreg_cov (c), j + 1, j + 1));
 815       tab_double (t, 3, this_row, 0, std_err, NULL);
 816       /*
 817         Standardized coefficient, i.e., regression coefficient
 818         if all variables had unit variance.
 819       */
 820       beta = sqrt (gsl_matrix_get (cov, j, j));
 821       beta *= linreg_coeff (c, j) /
 822         sqrt (gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1));
 823       tab_double (t, 4, this_row, 0, beta, NULL);
 824
 825       /*
 826         Test statistic for H0: coefficient is 0.
 827       */
 828       t_stat = linreg_coeff (c, j) / std_err;
 829       tab_double (t, 5, this_row, 0, t_stat, NULL);
 830       /*
 831         P values for the test statistic above.
 832       */
 833       pval =
 834         2 * gsl_cdf_tdist_Q (fabs (t_stat),
 835                              (double) (linreg_n_obs (c) - linreg_n_coeffs (c) - 1));
 836       tab_double (t, 6, this_row, 0, pval, NULL);
 837       ds_destroy (&tstr);
 838     }
 839   tab_title (t, _("Coefficients (%s)"), var_to_string (var));
 840   tab_submit (t);
 841 }
 842
 843 /*
 844   Display the ANOVA table.
 845 */
 846 static void
 847 reg_stats_anova (linreg * c, void *aux UNUSED, const struct variable *var)
 848 {
 849   int n_cols = 7;
 850   int n_rows = 4;
 851   const double msm = linreg_ssreg (c) / linreg_dfmodel (c);
 852   const double mse = linreg_mse (c);
 853   const double F = msm / mse;
 854   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
 855
 856   struct tab_table *t;
 857
 858   assert (c != NULL);
 859   t = tab_create (n_cols, n_rows);
 860   tab_headers (t, 2, 0, 1, 0);
 861
 862   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 863
 864   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 865   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 866   tab_vline (t, TAL_0, 1, 0, 0);
 867
 868   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 869   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 870   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 871   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 872   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 873
 874   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
 875   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
 876   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
 877
 878   /* Sums of Squares */
 879   tab_double (t, 2, 1, 0, linreg_ssreg (c), NULL);
 880   tab_double (t, 2, 3, 0, linreg_sst (c), NULL);
 881   tab_double (t, 2, 2, 0, linreg_sse (c), NULL);
 882
 883
 884   /* Degrees of freedom */
 885   tab_text_format (t, 3, 1, TAB_RIGHT, "%g", c->dfm);
 886   tab_text_format (t, 3, 2, TAB_RIGHT, "%g", c->dfe);
 887   tab_text_format (t, 3, 3, TAB_RIGHT, "%g", c->dft);
 888
 889   /* Mean Squares */
 890   tab_double (t, 4, 1, TAB_RIGHT, msm, NULL);
 891   tab_double (t, 4, 2, TAB_RIGHT, mse, NULL);
 892
 893   tab_double (t, 5, 1, 0, F, NULL);
 894
 895   tab_double (t, 6, 1, 0, pval, NULL);
 896
 897   tab_title (t, _("ANOVA (%s)"), var_to_string (var));
 898   tab_submit (t);
 899 }
 900
 901
 902 static void
 903 reg_stats_bcov (linreg * c, void *aux UNUSED, const struct variable *var)
 904 {
 905   int n_cols;
 906   int n_rows;
 907   int i;
 908   int k;
 909   int row;
 910   int col;
 911   const char *label;
 912   struct tab_table *t;
 913
 914   assert (c != NULL);
 915   n_cols = c->n_indeps + 1 + 2;
 916   n_rows = 2 * (c->n_indeps + 1);
 917   t = tab_create (n_cols, n_rows);
 918   tab_headers (t, 2, 0, 1, 0);
 919   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 920   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 921   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 922   tab_vline (t, TAL_0, 1, 0, 0);
 923   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 924   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 925   for (i = 0; i < linreg_n_coeffs (c); i++)
 926     {
 927       const struct variable *v = linreg_indep_var (c, i);
 928       label = var_to_string (v);
 929       tab_text (t, 2, i, TAB_CENTER, label);
 930       tab_text (t, i + 2, 0, TAB_CENTER, label);
 931       for (k = 1; k < linreg_n_coeffs (c); k++)
 932         {
 933           col = (i <= k) ? k : i;
 934           row = (i <= k) ? i : k;
 935           tab_double (t, k + 2, i, TAB_CENTER,
 936                       gsl_matrix_get (c->cov, row, col), NULL);
 937         }
 938     }
 939   tab_title (t, _("Coefficient Correlations (%s)"), var_to_string (var));
 940   tab_submit (t);
 941 }
 942
 943 static void
 944 statistics_keyword_output (void (*function) (linreg *, void *, const struct variable *var),
 945                            bool keyword, linreg * c, void *aux, const struct variable *var)
 946 {
 947   if (keyword)
 948     {
 949       (*function) (c, aux, var);
 950     }
 951 }