GLM: Mimic spss output with /INTERCEPT=EXCLUDE

diff --git a/src/language/stats/glm.c b/src/language/stats/glm.c
index d01d41805758705adae7cf6f94f5ffcf20cc5c17..92e35ccd27e1f5f4c1f8082d63d402e90f05ca4d 100644 (file)
--- a/src/language/stats/glm.c
+++ b/src/language/stats/glm.c
@@ -717,7 +717,7 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
     cmd->wv ? var_get_print_format (cmd->wv) : &F_8_0;
 
   double n_total, mean;
-  double df_corr = 0.0;
+  double df_corr = 1.0;
   double mse = 0;
 
   int f;
@@ -727,9 +727,9 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
   struct tab_table *t;
 
   const int nc = 6;
-  int nr = heading_rows + 4 + cmd->n_interactions;
+  int nr = heading_rows + 3 + cmd->n_interactions;
   if (cmd->intercept)
-    nr++;
+    nr += 2;
 
   msg (MW, "GLM is experimental.  Do not rely on these results.");
   t = tab_create (nc, nr);
@@ -755,27 +755,29 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
 
   moments_calculate (ws->totals, &n_total, &mean, NULL, NULL, NULL);
 
-  if (cmd->intercept)
-    df_corr += 1.0;
-
   df_corr += categoricals_df_total (ws->cats);
 
-  mse = gsl_vector_get (ws->ssq, 0) / (n_total - df_corr);
-
   r = heading_rows;
-  tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Model"));
+  if (cmd->intercept)
+    tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Model"));
+  else
+    tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Model"));
 
   r++;
 
+  mse = gsl_vector_get (ws->ssq, 0) / (n_total - df_corr);
+
+  const double intercept_ssq = pow2 (mean * n_total) / n_total;
+
+  double ssq_effects = 0.0;
   if (cmd->intercept)
     {
-      const double intercept = pow2 (mean * n_total) / n_total;
       const double df = 1.0;
-      const double F = intercept / df / mse;
+      const double F = intercept_ssq / df / mse;
       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Intercept"));
-      tab_double (t, 1, r, 0, intercept, NULL);
+      tab_double (t, 1, r, 0, intercept_ssq, NULL);
       tab_double (t, 2, r, 0, 1.00, wfmt);
-      tab_double (t, 3, r, 0, intercept / df, NULL);
+      tab_double (t, 3, r, 0, intercept_ssq / df, NULL);
       tab_double (t, 4, r, 0, F, NULL);
       tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr),
                  NULL);
@@ -785,8 +787,17 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
   for (f = 0; f < cmd->n_interactions; ++f)
     {
       struct string str = DS_EMPTY_INITIALIZER;
-      const double df = categoricals_df (ws->cats, f);
-      const double ssq = gsl_vector_get (ws->ssq, f + 1);
+      double df = categoricals_df (ws->cats, f);
+
+      double ssq = gsl_vector_get (ws->ssq, f + 1);
+      ssq_effects += ssq;
+
+      if (! cmd->intercept) 
+       {
+         df++;
+         ssq += intercept_ssq;
+       }
+
       const double F = ssq / df / mse;
       interaction_to_string (cmd->interactions[f], &str);
       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, ds_cstr (&str));
@@ -803,9 +814,14 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
     }
 
   {
-    /* Corrected Model */
-    const double df = df_corr - 1.0;
-    const double ssq = ws->total_ssq - gsl_vector_get (ws->ssq, 0);
+    /* Model / Corrected Model */
+    double df = df_corr;
+    double ssq = ws->total_ssq - gsl_vector_get (ws->ssq, 0);
+    if ( cmd->intercept )
+      df --;
+    else
+      ssq += intercept_ssq;
+
     const double F = ssq / df / mse;
     tab_double (t, 1, heading_rows, 0, ssq, NULL);
     tab_double (t, 2, heading_rows, 0, df, wfmt);
@@ -826,24 +842,21 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
     tab_double (t, 3, r++, 0, mse, NULL);
   }
 
+  {
+    tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Total"));
+    tab_double (t, 1, r, 0, ws->total_ssq + intercept_ssq, NULL);
+    tab_double (t, 2, r, 0, n_total, wfmt);
+    
+    r++;
+  }
+
   if (cmd->intercept)
     {
-      const double intercept = pow2 (mean * n_total) / n_total;
-      const double ssq = intercept + ws->total_ssq;
-
-      tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Total"));
-      tab_double (t, 1, r, 0, ssq, NULL);
-      tab_double (t, 2, r, 0, n_total, wfmt);
-
-      r++;
+      tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Total"));
+      tab_double (t, 1, r, 0, ws->total_ssq, NULL);
+      tab_double (t, 2, r, 0, n_total - 1.0, wfmt);
     }
 
-  tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Total"));
-
-
-  tab_double (t, 1, r, 0, ws->total_ssq, NULL);
-  tab_double (t, 2, r, 0, n_total - 1.0, wfmt);
-
   tab_submit (t);
 }
 

diff --git a/tests/language/stats/glm.at b/tests/language/stats/glm.at
index 28a86d359f66d969e69c3a74136d420731ad1262..ec9a6799e379cd31f4852f07ac256e1b2d93bb35 100644 (file)
--- a/tests/language/stats/glm.at
+++ b/tests/language/stats/glm.at
@@ -263,3 +263,80 @@ Corrected Total,399.474,37,,,
 ])
 
 AT_CLEANUP
+
+
+
+AT_SETUP([GLM excluded intercept])
+
+dnl  The following example comes from 
+dnl
+dnl Rudolf N. Cardinal
+dnl Graduate-level statistics for psychology and neuroscience
+dnl ANOVA in practice, and complex ANOVA designs
+dnl Version of 2 May 2004
+dnl
+dnl Downloaded from: http://egret.psychol.cam.ac.uk/psychology/graduate/Guide_to_ANOVA.pdf
+
+AT_DATA([intercept-exclude.sps], [dnl
+set format = F20.3.
+
+data list notable list /depvar * A *.
+begin data.
+10     1
+14     1
+8      1
+7      1
+2      1
+10     1
+1      1
+3      1
+2      1
+8.5    1
+14.29  2
+18.49  2
+12.46  2
+11.63  2
+6.66   2
+14.02  2
+5.66   2
+7.06   2
+6.37   2
+13.26  2
+end data.
+
+GLM depvar by A
+   /intercept = exclude
+  .
+
+
+GLM depvar by A
+   /intercept = include
+  .
+
+])
+
+AT_CHECK([pspp -O format=csv intercept-exclude.sps], [0],
+  [dnl
+warning: GLM is experimental.  Do not rely on these results.
+
+Table: Tests of Between-Subjects Effects
+Source,Type III Sum of Squares,df,Mean Square,F,Sig.
+Model,1636.826,2,818.413,43.556,.000
+A,1636.826,2,818.413,43.556,.000
+Error,338.216,18,18.790,,
+Total,1975.042,20,,,
+
+warning: GLM is experimental.  Do not rely on these results.
+
+Table: Tests of Between-Subjects Effects
+Source,Type III Sum of Squares,df,Mean Square,F,Sig.
+Corrected Model,98.568,1,98.568,5.246,.034
+Intercept,1538.258,1,1538.258,81.867,.000
+A,98.568,1,98.568,5.246,.034
+Error,338.216,18,18.790,,
+Total,1975.042,20,,,
+Corrected Total,436.784,19,,,
+])
+
+AT_CLEANUP
+