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