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