Regression: Use alternative formula for adjusted R-square
[pspp] / src / language / stats / regression.c
index 9974ea4c2c2155c140e072e48c8e27701b2e8033..e32abb3eb57485313a3a13a6aa94b82ae71f76da 100644 (file)
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2005, 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
+   Copyright (C) 2005, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -196,16 +196,15 @@ reg_get_name (const struct dictionary *dict, const char *prefix)
 static bool
 regression_trns_free (void *t_)
 {
-  bool result = true;
   struct reg_trns *t = t_;
 
   if (t->trns_id == t->n_trns)
     {
-      result = linreg_free (t->c);
+      linreg_unref (t->c);
     }
   free (t);
 
-  return result;
+  return true;
 }
 
 static void
@@ -249,6 +248,7 @@ subcommand_save (const struct regression *cmd)
         {
           if ((*lc)->depvar != NULL)
             {
+             (*lc)->refcnt++;
               if (cmd->resid)
                 {
                   reg_save_var (cmd->ds, "RES", regression_trns_resid_proc, *lc,
@@ -267,8 +267,10 @@ subcommand_save (const struct regression *cmd)
 int
 cmd_regression (struct lexer *lexer, struct dataset *ds)
 {
+  int k;
   struct regression regression;
   const struct dictionary *dict = dataset_dict (ds);
+  bool save;
 
   memset (&regression, 0, sizeof (struct regression));
 
@@ -387,29 +389,47 @@ cmd_regression (struct lexer *lexer, struct dataset *ds)
 
   regression.models = xcalloc (regression.n_dep_vars, sizeof *regression.models);
 
+  save = regression.pred || regression.resid;
+  if (save)
+    {
+      if (proc_make_temporary_transformations_permanent (ds))
+        msg (SW, _("REGRESSION with SAVE ignores TEMPORARY.  "
+                   "Temporary transformations will be made permanent."));
+    }
+
   {
     struct casegrouper *grouper;
     struct casereader *group;
     bool ok;
     
-    grouper = casegrouper_create_splits (proc_open (ds), dict);
+    grouper = casegrouper_create_splits (proc_open_filtering (ds, !save),
+                                         dict);
     while (casegrouper_get_next_group (grouper, &group))
       run_regression (&regression, group);
     ok = casegrouper_destroy (grouper);
     ok = proc_commit (ds) && ok;
   }
 
-  if (regression.pred || regression.resid )
-    subcommand_save (&regression);
+  if (save)
+    {
+      subcommand_save (&regression);
+    }
  
 
+  for (k = 0; k < regression.n_dep_vars; k++)
+    linreg_unref (regression.models[k]);
   free (regression.models);
   free (regression.vars);
   free (regression.dep_vars);
   return CMD_SUCCESS;
   
  error:
+  if (regression.models)
+   {
+  for (k = 0; k < regression.n_dep_vars; k++)
+    linreg_unref (regression.models[k]);
   free (regression.models);
+   }
   free (regression.vars);
   free (regression.dep_vars);
   return CMD_FAILURE;
@@ -495,7 +515,7 @@ identify_indep_vars (const struct regression *cmd,
        There is only one independent variable, and it is the same
        as the dependent variable. Print a warning and continue.
       */
-      msg (SE,
+      msg (SW,
           gettext ("The dependent variable is equal to the independent variable." 
                    "The least squares line is therefore Y=X." 
                    "Standard errors and related statistics may be meaningless."));
@@ -581,23 +601,24 @@ fill_covariance (gsl_matrix *cov, struct covariance *all_cov,
 /*
   STATISTICS subcommand output functions.
 */
-static void reg_stats_r (linreg *, void *);
-static void reg_stats_coeff (linreg *, void *);
-static void reg_stats_anova (linreg *, void *);
-static void reg_stats_bcov (linreg *, void *);
+static void reg_stats_r (linreg *, void *, const struct variable *);
+static void reg_stats_coeff (linreg *, void *, const struct variable *);
+static void reg_stats_anova (linreg *, void *, const struct variable *);
+static void reg_stats_bcov (linreg *, void *, const struct variable *);
 
-static void statistics_keyword_output (void (*)(linreg *, void *),
-                                      bool, linreg *, void *);
+static void statistics_keyword_output (void (*)(linreg *, void *, const struct variable *),
+                                      bool, linreg *, void *, const struct variable *);
 
 
 
 static void
-subcommand_statistics (const struct regression *cmd , linreg * c, void *aux)
+subcommand_statistics (const struct regression *cmd , linreg * c, void *aux,
+                      const struct variable *var)
 {
-  statistics_keyword_output (reg_stats_r, cmd->r, c, aux);
-  statistics_keyword_output (reg_stats_anova, cmd->anova, c, aux);
-  statistics_keyword_output (reg_stats_coeff, cmd->coeff, c, aux);
-  statistics_keyword_output (reg_stats_bcov, cmd->bcov, c, aux);
+  statistics_keyword_output (reg_stats_r, cmd->r, c, aux, var);
+  statistics_keyword_output (reg_stats_anova, cmd->anova, c, aux, var);
+  statistics_keyword_output (reg_stats_coeff, cmd->coeff, c, aux, var);
+  statistics_keyword_output (reg_stats_bcov, cmd->bcov, c, aux, var);
 }
 
 
@@ -612,7 +633,6 @@ run_regression (const struct regression *cmd, struct casereader *input)
   struct covariance *cov;
   const struct variable **vars;
   const struct variable **all_vars;
-  const struct variable *dep_var;
   struct casereader *reader;
   size_t n_all_vars;
 
@@ -639,9 +659,9 @@ run_regression (const struct regression *cmd, struct casereader *input)
   for (k = 0; k < cmd->n_dep_vars; k++)
     {
       double n_data;
-
+      const struct variable *dep_var = cmd->dep_vars[k];
       gsl_matrix *this_cm;
-      dep_var = cmd->dep_vars[k];
+
       n_indep = identify_indep_vars (cmd, vars, dep_var);
       
       this_cm = gsl_matrix_alloc (n_indep + 1, n_indep + 1);
@@ -672,14 +692,14 @@ run_regression (const struct regression *cmd, struct casereader *input)
          
          if (!taint_has_tainted_successor (casereader_get_taint (input)))
            {
-             subcommand_statistics (cmd, models[k], this_cm);
+             subcommand_statistics (cmd, models[k], this_cm, dep_var);
            }
        }
       else
        {
           msg (SE,
               _("No valid data found. This command was skipped."));
-          linreg_free (models[k]);
+          linreg_unref (models[k]);
           models[k] = NULL;
        }
       gsl_matrix_free (this_cm);
@@ -698,7 +718,7 @@ run_regression (const struct regression *cmd, struct casereader *input)
 
 
 static void
-reg_stats_r (linreg *c, void *aux UNUSED)
+reg_stats_r (linreg *c, void *aux UNUSED, const struct variable *var)
 {
   struct tab_table *t;
   int n_rows = 2;
@@ -709,7 +729,8 @@ reg_stats_r (linreg *c, void *aux UNUSED)
 
   assert (c != NULL);
   rsq = linreg_ssreg (c) / linreg_sst (c);
-  adjrsq = 1.0 - (1.0 - rsq) * (linreg_n_obs (c) - 1.0) / (linreg_n_obs (c) - linreg_n_coeffs (c));
+  adjrsq = rsq - 
+    (1.0 - rsq) * linreg_n_coeffs (c) / (linreg_n_obs (c) - linreg_n_coeffs (c) - 1);
   std_error = sqrt (linreg_mse (c));
   t = tab_create (n_cols, n_rows);
   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
@@ -725,7 +746,7 @@ reg_stats_r (linreg *c, void *aux UNUSED)
   tab_double (t, 2, 1, TAB_RIGHT, rsq, NULL);
   tab_double (t, 3, 1, TAB_RIGHT, adjrsq, NULL);
   tab_double (t, 4, 1, TAB_RIGHT, std_error, NULL);
-  tab_title (t, _("Model Summary"));
+  tab_title (t, _("Model Summary (%s)"), var_to_string (var));
   tab_submit (t);
 }
 
@@ -733,7 +754,7 @@ reg_stats_r (linreg *c, void *aux UNUSED)
   Table showing estimated regression coefficients.
 */
 static void
-reg_stats_coeff (linreg * c, void *aux_)
+reg_stats_coeff (linreg * c, void *aux_, const struct variable *var)
 {
   size_t j;
   int n_cols = 7;
@@ -812,11 +833,11 @@ reg_stats_coeff (linreg * c, void *aux_)
       */
       pval =
        2 * gsl_cdf_tdist_Q (fabs (t_stat),
-                            (double) (linreg_n_obs (c) - linreg_n_coeffs (c)));
+                            (double) (linreg_n_obs (c) - linreg_n_coeffs (c) - 1));
       tab_double (t, 6, this_row, 0, pval, NULL);
       ds_destroy (&tstr);
     }
-  tab_title (t, _("Coefficients"));
+  tab_title (t, _("Coefficients (%s)"), var_to_string (var));
   tab_submit (t);
 }
 
@@ -824,7 +845,7 @@ reg_stats_coeff (linreg * c, void *aux_)
   Display the ANOVA table.
 */
 static void
-reg_stats_anova (linreg * c, void *aux UNUSED)
+reg_stats_anova (linreg * c, void *aux UNUSED, const struct variable *var)
 {
   int n_cols = 7;
   int n_rows = 4;
@@ -874,13 +895,13 @@ reg_stats_anova (linreg * c, void *aux UNUSED)
 
   tab_double (t, 6, 1, 0, pval, NULL);
 
-  tab_title (t, _("ANOVA"));
+  tab_title (t, _("ANOVA (%s)"), var_to_string (var));
   tab_submit (t);
 }
 
 
 static void
-reg_stats_bcov (linreg * c, void *aux UNUSED)
+reg_stats_bcov (linreg * c, void *aux UNUSED, const struct variable *var)
 {
   int n_cols;
   int n_rows;
@@ -916,16 +937,16 @@ reg_stats_bcov (linreg * c, void *aux UNUSED)
                       gsl_matrix_get (c->cov, row, col), NULL);
        }
     }
-  tab_title (t, _("Coefficient Correlations"));
+  tab_title (t, _("Coefficient Correlations (%s)"), var_to_string (var));
   tab_submit (t);
 }
 
 static void
-statistics_keyword_output (void (*function) (linreg *, void *),
-                          bool keyword, linreg * c, void *aux)
+statistics_keyword_output (void (*function) (linreg *, void *, const struct variable *var),
+                          bool keyword, linreg * c, void *aux, const struct variable *var)
 {
   if (keyword)
     {
-      (*function) (c, aux);
+      (*function) (c, aux, var);
     }
 }