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