X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Foneway.c;h=70e6a3d0d84c909f9c5f20a477c25a7f0b147c72;hb=ab90b7f68c2151df4caaf376c7200443d1cf9f5e;hp=a63d1a29e532876ebadb1fbf4d2f8100e0f83021;hpb=fe8dc2171009e90d2335f159d05f7e6660e24780;p=pspp

diff --git a/src/language/stats/oneway.c b/src/language/stats/oneway.c
index a63d1a29e5..70e6a3d0d8 100644
--- a/src/language/stats/oneway.c
+++ b/src/language/stats/oneway.c
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 1997-9, 2000, 2007, 2009, 2010, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-9, 2000, 2007, 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
@@ -37,7 +37,9 @@
 #include "libpspp/misc.h"
 #include "libpspp/taint.h"
 #include "linreg/sweep.h"
+#include "tukey/tukey.h"
 #include "math/categoricals.h"
+#include "math/interaction.h"
 #include "math/covariance.h"
 #include "math/levene.h"
 #include "math/moments.h"
@@ -45,7 +47,36 @@
 
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
+#define N_(msgid) msgid
 
+/* Workspace variable for each dependent variable */
+struct per_var_ws
+{
+  struct interaction *iact;
+  struct categoricals *cat;
+  struct covariance *cov;
+  struct levene *nl;
+
+  double n;
+
+  double sst;
+  double sse;
+  double ssa;
+
+  int n_groups;
+
+  double mse;
+};
+
+/* Per category data */
+struct descriptive_data
+{
+  const struct variable *var;
+  struct moments1 *mom;
+
+  double minimum;
+  double maximum;
+};
 
 enum missing_type
   {
@@ -70,8 +101,28 @@ struct contrasts_node
 {
   struct ll ll; 
   struct ll_list coefficient_list;
+};
+
+
+struct oneway_spec;
+
+typedef double df_func (const struct per_var_ws *pvw, const struct moments1 *mom_i, const struct moments1 *mom_j);
+typedef double ts_func (int k, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err);
+typedef double p1tail_func (double ts, double df1, double df2);
 
-  bool bad_count; /* True if the number of coefficients does not equal the number of groups */
+typedef double pinv_func (double std_err, double alpha, double df, int k, const struct moments1 *mom_i, const struct moments1 *mom_j);
+
+
+struct posthoc
+{
+  const char *syntax;
+  const char *label;
+
+  df_func *dff;
+  ts_func *tsf;
+  p1tail_func *p1f;
+
+  pinv_func *pinv;
 };
 
 struct oneway_spec
@@ -92,33 +143,227 @@ struct oneway_spec
   /* The weight variable */
   const struct variable *wv;
 
+  /* The confidence level for multiple comparisons */
+  double alpha;
+
+  int *posthoc;
+  int n_posthoc;
 };
 
-/* Per category data */
-struct descriptive_data
+static double
+df_common (const struct per_var_ws *pvw, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
 {
-  const struct variable *var;
-  struct moments1 *mom;
+  return  pvw->n - pvw->n_groups;
+}
 
-  double minimum;
-  double maximum;
-};
+static double
+df_individual (const struct per_var_ws *pvw UNUSED, const struct moments1 *mom_i, const struct moments1 *mom_j)
+{
+  double n_i, var_i;
+  double n_j, var_j;
+  double nom,denom;
 
-/* Workspace variable for each dependent variable */
-struct per_var_ws
+  moments1_calculate (mom_i, &n_i, NULL, &var_i, 0, 0);  
+  moments1_calculate (mom_j, &n_j, NULL, &var_j, 0, 0);
+  
+  if ( n_i <= 1.0 || n_j <= 1.0)
+    return SYSMIS;
+
+  nom = pow2 (var_i/n_i + var_j/n_j);
+  denom = pow2 (var_i/n_i) / (n_i - 1) + pow2 (var_j/n_j) / (n_j - 1);
+
+  return nom / denom;
+}
+
+static double lsd_pinv (double std_err, double alpha, double df, int k UNUSED, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
 {
-  struct categoricals *cat;
-  struct covariance *cov;
+  return std_err * gsl_cdf_tdist_Pinv (1.0 - alpha / 2.0, df);
+}
 
-  double sst;
-  double sse;
-  double ssa;
+static double bonferroni_pinv (double std_err, double alpha, double df, int k, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
+{
+  const int m = k * (k - 1) / 2;
+  return std_err * gsl_cdf_tdist_Pinv (1.0 - alpha / (2.0 * m), df);
+}
 
-  int n_groups;
+static double sidak_pinv (double std_err, double alpha, double df, int k, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
+{
+  const double m = k * (k - 1) / 2;
+  double lp = 1.0 - exp (log (1.0 - alpha) / m ) ;
+  return std_err * gsl_cdf_tdist_Pinv (1.0 - lp / 2.0, df);
+}
+
+static double tukey_pinv (double std_err, double alpha, double df, int k, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
+{
+  if ( k < 2 || df < 2)
+    return SYSMIS;
+
+  return std_err / sqrt (2.0)  * qtukey (1 - alpha, 1.0, k, df, 1, 0);
+}
+
+static double scheffe_pinv (double std_err, double alpha, double df, int k, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
+{
+  double x = (k - 1) * gsl_cdf_fdist_Pinv (1.0 - alpha, k - 1, df);
+  return std_err * sqrt (x);
+}
+
+static double gh_pinv (double std_err UNUSED, double alpha, double df, int k, const struct moments1 *mom_i, const struct moments1 *mom_j)
+{
+  double n_i, mean_i, var_i;
+  double n_j, mean_j, var_j;
+  double m;
+
+  moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);  
+  moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+  m = sqrt ((var_i/n_i + var_j/n_j) / 2.0);
+
+  if ( k < 2 || df < 2)
+    return SYSMIS;
+
+  return m * qtukey (1 - alpha, 1.0, k, df, 1, 0);
+}
+
+
+static double 
+multiple_comparison_sig (double std_err,
+				       const struct per_var_ws *pvw,
+				       const struct descriptive_data *dd_i, const struct descriptive_data *dd_j,
+				       const struct posthoc *ph)
+{
+  int k = pvw->n_groups;
+  double df = ph->dff (pvw, dd_i->mom, dd_j->mom);
+  double ts = ph->tsf (k, dd_i->mom, dd_j->mom, std_err);
+  if ( df == SYSMIS)
+    return SYSMIS;
+  return  ph->p1f (ts, k - 1, df);
+}
+
+static double 
+mc_half_range (const struct oneway_spec *cmd, const struct per_var_ws *pvw, double std_err, const struct descriptive_data *dd_i, const struct descriptive_data *dd_j, const struct posthoc *ph)
+{
+  int k = pvw->n_groups;
+  double df = ph->dff (pvw, dd_i->mom, dd_j->mom);
+  if ( df == SYSMIS)
+    return SYSMIS;
+
+  return ph->pinv (std_err, cmd->alpha, df, k, dd_i->mom, dd_j->mom);
+}
+
+static double tukey_1tailsig (double ts, double df1, double df2)
+{
+  double twotailedsig;
+
+  if (df2 < 2 || df1 < 1)
+    return SYSMIS;
+
+  twotailedsig = 1.0 - ptukey (ts, 1.0, df1 + 1, df2, 1, 0);
+
+  return twotailedsig / 2.0;
+}
+
+static double lsd_1tailsig (double ts, double df1 UNUSED, double df2)
+{
+  return ts < 0 ? gsl_cdf_tdist_P (ts, df2) : gsl_cdf_tdist_Q (ts, df2);
+}
+
+static double sidak_1tailsig (double ts, double df1, double df2)
+{
+  double ex = (df1 + 1.0) * df1 / 2.0;
+  double lsd_sig = 2 * lsd_1tailsig (ts, df1, df2);
+
+  return 0.5 * (1.0 - pow (1.0 - lsd_sig, ex));
+}
+
+static double bonferroni_1tailsig (double ts, double df1, double df2)
+{
+  const int m = (df1 + 1) * df1 / 2;
+
+  double p = ts < 0 ? gsl_cdf_tdist_P (ts, df2) : gsl_cdf_tdist_Q (ts, df2);
+  p *= m;
+
+  return p > 0.5 ? 0.5 : p;
+}
+
+static double scheffe_1tailsig (double ts, double df1, double df2)
+{
+  return 0.5 * gsl_cdf_fdist_Q (ts, df1, df2);
+}
+
+
+static double tukey_test_stat (int k UNUSED, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err)
+{
+  double ts;
+  double n_i, mean_i, var_i;
+  double n_j, mean_j, var_j;
+
+  moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);  
+  moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+  ts =  (mean_i - mean_j) / std_err;
+  ts = fabs (ts) * sqrt (2.0);
+
+  return ts;
+}
+
+static double lsd_test_stat (int k UNUSED, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err)
+{
+  double n_i, mean_i, var_i;
+  double n_j, mean_j, var_j;
+
+  moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);  
+  moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+  return (mean_i - mean_j) / std_err;
+}
+
+static double scheffe_test_stat (int k, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err)
+{
+  double t;
+  double n_i, mean_i, var_i;
+  double n_j, mean_j, var_j;
+
+  moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);  
+  moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+  t = (mean_i - mean_j) / std_err;
+  t = pow2 (t);
+  t /= k - 1;
+
+  return t;
+}
+
+static double gh_test_stat (int k UNUSED, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err UNUSED)
+{
+  double ts;
+  double thing;
+  double n_i, mean_i, var_i;
+  double n_j, mean_j, var_j;
+
+  moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);  
+  moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+  thing = var_i / n_i + var_j / n_j;
+  thing /= 2.0;
+  thing = sqrt (thing);
+
+  ts = (mean_i - mean_j) / thing;
+
+  return fabs (ts);
+}
+
+
+
+static const struct posthoc ph_tests [] = 
+  {
+    { "LSD",        N_("LSD"),          df_common, lsd_test_stat,     lsd_1tailsig,          lsd_pinv},
+    { "TUKEY",      N_("Tukey HSD"),    df_common, tukey_test_stat,   tukey_1tailsig,        tukey_pinv},
+    { "BONFERRONI", N_("Bonferroni"),   df_common, lsd_test_stat,     bonferroni_1tailsig,   bonferroni_pinv},
+    { "SCHEFFE",    N_("Scheffé"),      df_common, scheffe_test_stat, scheffe_1tailsig,      scheffe_pinv},
+    { "GH",         N_("Games-Howell"), df_individual, gh_test_stat,  tukey_1tailsig,        gh_pinv},
+    { "SIDAK",      N_("Šidák"),        df_common, lsd_test_stat,     sidak_1tailsig,        sidak_pinv}
+  };
 
-  double mse;
-  double levene_w;
-};
 
 struct oneway_workspace
 {
@@ -140,6 +385,37 @@ static void show_homogeneity (const struct oneway_spec *, const struct oneway_wo
 static void output_oneway (const struct oneway_spec *, struct oneway_workspace *ws);
 static void run_oneway (const struct oneway_spec *cmd, struct casereader *input, const struct dataset *ds);
 
+
+static void
+destroy_coeff_list (struct contrasts_node *coeff_list)
+{
+  struct coeff_node *cn = NULL;
+  struct coeff_node *cnx = NULL;
+  struct ll_list *cl = &coeff_list->coefficient_list;
+  
+  ll_for_each_safe (cn, cnx, struct coeff_node, ll, cl)
+    {
+      free (cn);
+    }
+  
+  free (coeff_list);
+}
+
+static void
+oneway_cleanup (struct oneway_spec *cmd)
+{
+  struct contrasts_node *coeff_list  = NULL;
+  struct contrasts_node *coeff_next  = NULL;
+  ll_for_each_safe (coeff_list, coeff_next, struct contrasts_node, ll, &cmd->contrast_list)
+    {
+      destroy_coeff_list (coeff_list);
+    }
+
+  free (cmd->posthoc);
+}
+
+
+
 int
 cmd_oneway (struct lexer *lexer, struct dataset *ds)
 {
@@ -152,6 +428,9 @@ cmd_oneway (struct lexer *lexer, struct dataset *ds)
   oneway.missing_type = MISS_ANALYSIS;
   oneway.exclude = MV_ANY;
   oneway.wv = dict_get_weight (dict);
+  oneway.alpha = 0.05;
+  oneway.posthoc = NULL;
+  oneway.n_posthoc = 0;
 
   ll_init (&oneway.contrast_list);
 
@@ -198,6 +477,45 @@ cmd_oneway (struct lexer *lexer, struct dataset *ds)
 		}
 	    }
 	}
+      else if (lex_match_id (lexer, "POSTHOC"))
+	{
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
+	    {
+	      int p;
+	      bool method = false;
+	      for (p = 0 ; p < sizeof (ph_tests) / sizeof (struct posthoc); ++p)
+		{
+		  if (lex_match_id (lexer, ph_tests[p].syntax))
+		    {
+		      oneway.n_posthoc++;
+		      oneway.posthoc = xrealloc (oneway.posthoc, sizeof (*oneway.posthoc) * oneway.n_posthoc);
+		      oneway.posthoc[oneway.n_posthoc - 1] = p;
+		      method = true;
+		      break;
+		    }
+		}
+	      if ( method == false)
+		{
+		  if (lex_match_id (lexer, "ALPHA"))
+		    {
+		      if ( !lex_force_match (lexer, T_LPAREN))
+			goto error;
+		      lex_force_num (lexer);
+		      oneway.alpha = lex_number (lexer);
+		      lex_get (lexer);
+		      if ( !lex_force_match (lexer, T_RPAREN))
+			goto error;
+		    }
+		  else
+		    {
+		      msg (SE, _("The post hoc analysis method %s is not supported."), lex_tokcstr (lexer));
+		      lex_error (lexer, NULL);
+		      goto error;
+		    }
+		}
+	    }
+	}
       else if (lex_match_id (lexer, "CONTRAST"))
 	{
 	  struct contrasts_node *cl = xzalloc (sizeof *cl);
@@ -219,6 +537,7 @@ cmd_oneway (struct lexer *lexer, struct dataset *ds)
 		}
 	      else
 		{
+		  destroy_coeff_list (cl);
 		  lex_error (lexer, NULL);
 		  goto error;
 		}
@@ -274,10 +593,12 @@ cmd_oneway (struct lexer *lexer, struct dataset *ds)
     ok = proc_commit (ds) && ok;
   }
 
+  oneway_cleanup (&oneway);
   free (oneway.vars);
   return CMD_SUCCESS;
 
  error:
+  oneway_cleanup (&oneway);
   free (oneway.vars);
   return CMD_FAILURE;
 }
@@ -307,7 +628,7 @@ dd_destroy (struct descriptive_data *dd)
 }
 
 static void *
-makeit (void *aux1, void *aux2 UNUSED)
+makeit (const void *aux1, void *aux2 UNUSED)
 {
   const struct variable *var = aux1;
 
@@ -317,27 +638,25 @@ makeit (void *aux1, void *aux2 UNUSED)
 }
 
 static void 
-updateit (void *user_data, 
-	  enum mv_class exclude,
-	  const struct variable *wv, 
-	  const struct variable *catvar UNUSED,
-	  const struct ccase *c,
-	  void *aux1, void *aux2)
+killit (const void *aux1 UNUSED, void *aux2 UNUSED, void *user_data)
 {
   struct descriptive_data *dd = user_data;
 
-  const struct variable *varp = aux1;
+  dd_destroy (dd);
+}
 
-  const union value *valx = case_data (c, varp);
 
-  struct descriptive_data *dd_total = aux2;
+static void 
+updateit (const void *aux1, void *aux2, void *user_data,
+	  const struct ccase *c, double weight)
+{
+  struct descriptive_data *dd = user_data;
 
-  double weight;
+  const struct variable *varp = aux1;
 
-  if ( var_is_value_missing (varp, valx, exclude))
-    return;
+  const union value *valx = case_data (c, varp);
 
-  weight = wv != NULL ? case_data (c, wv)->f : 1.0;
+  struct descriptive_data *dd_total = aux2;
 
   moments1_add (dd->mom, valx->f, weight);
   if (valx->f < dd->minimum)
@@ -384,15 +703,25 @@ run_oneway (const struct oneway_spec *cmd,
 
   for (v = 0; v < cmd->n_vars; ++v)
     {
-      ws.vws[v].cat = categoricals_create (&cmd->indep_var, 1, cmd->wv,
-                                           cmd->exclude, makeit, updateit,
-                                           CONST_CAST (struct variable *,
-                                                       cmd->vars[v]),
-                                           ws.dd_total[v]);
+      struct payload payload;
+      payload.create = makeit;
+      payload.update = updateit;
+      payload.calculate = NULL;
+      payload.destroy = killit;
+
+      ws.vws[v].iact = interaction_create (cmd->indep_var);
+      ws.vws[v].cat = categoricals_create (&ws.vws[v].iact, 1, cmd->wv,
+                                           cmd->exclude, cmd->exclude);
+
+      categoricals_set_payload (ws.vws[v].cat, &payload, 
+				CONST_CAST (struct variable *, cmd->vars[v]),
+				ws.dd_total[v]);
+
 
       ws.vws[v].cov = covariance_2pass_create (1, &cmd->vars[v],
 					       ws.vws[v].cat, 
 					       cmd->wv, cmd->exclude);
+      ws.vws[v].nl = levene_create (var_get_width (cmd->indep_var), NULL);
     }
 
   c = casereader_peek (input, 0);
@@ -413,20 +742,11 @@ run_oneway (const struct oneway_spec *cmd,
                                               cmd->exclude, NULL, NULL);
   input = casereader_create_filter_weight (input, dict, NULL, NULL);
 
-
-  if (cmd->stats & STATS_HOMOGENEITY)
-    for (v = 0; v < cmd->n_vars; ++v)
-      {
-	struct per_var_ws *pvw = &ws.vws[v];
-
-	pvw->levene_w = levene (input, cmd->indep_var, cmd->vars[v], cmd->wv, cmd->exclude);
-      }
-
   reader = casereader_clone (input);
-
   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
     {
       int i;
+      double w = dict_get_case_weight (dict, c, NULL);
 
       for (i = 0; i < cmd->n_vars; ++i)
 	{
@@ -441,13 +761,16 @@ run_oneway (const struct oneway_spec *cmd,
 	    }
 
 	  covariance_accumulate_pass1 (pvw->cov, c);
+	  levene_pass_one (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
 	}
     }
   casereader_destroy (reader);
+
   reader = casereader_clone (input);
   for ( ; (c = casereader_read (reader) ); case_unref (c))
     {
       int i;
+      double w = dict_get_case_weight (dict, c, NULL);
       for (i = 0; i < cmd->n_vars; ++i)
 	{
 	  struct per_var_ws *pvw = &ws.vws[i];
@@ -461,44 +784,83 @@ run_oneway (const struct oneway_spec *cmd,
 	    }
 
 	  covariance_accumulate_pass2 (pvw->cov, c);
+	  levene_pass_two (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
 	}
     }
   casereader_destroy (reader);
 
+  reader = casereader_clone (input);
+  for ( ; (c = casereader_read (reader) ); case_unref (c))
+    {
+      int i;
+      double w = dict_get_case_weight (dict, c, NULL);
+
+      for (i = 0; i < cmd->n_vars; ++i)
+	{
+	  struct per_var_ws *pvw = &ws.vws[i];
+	  const struct variable *v = cmd->vars[i];
+	  const union value *val = case_data (c, v);
+
+	  if ( MISS_ANALYSIS == cmd->missing_type)
+	    {
+	      if ( var_is_value_missing (v, val, cmd->exclude))
+		continue;
+	    }
+
+	  levene_pass_three (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
+	}
+    }
+  casereader_destroy (reader);
+
+
   for (v = 0; v < cmd->n_vars; ++v)
     {
+      const gsl_matrix *ucm;
+      gsl_matrix *cm;
       struct per_var_ws *pvw = &ws.vws[v];
-      gsl_matrix *cm = covariance_calculate_unnormalized (pvw->cov);
       const struct categoricals *cats = covariance_get_categoricals (pvw->cov);
+      const bool ok = categoricals_sane (cats);
 
-      double n;
-      moments1_calculate (ws.dd_total[v]->mom, &n, NULL, NULL, NULL, NULL);
+      if ( ! ok)
+	{
+	  msg (MW, 
+	       _("Dependent variable %s has no non-missing values.  No analysis for this variable will be done."),
+	       var_get_name (cmd->vars[v]));
+	  continue;
+	}
 
-      pvw->sst = gsl_matrix_get (cm, 0, 0);
+      ucm = covariance_calculate_unnormalized (pvw->cov);
 
-      //      gsl_matrix_fprintf (stdout, cm, "%g ");
+      cm = gsl_matrix_alloc (ucm->size1, ucm->size2);
+      gsl_matrix_memcpy (cm, ucm);
+
+      moments1_calculate (ws.dd_total[v]->mom, &pvw->n, NULL, NULL, NULL, NULL);
+
+      pvw->sst = gsl_matrix_get (cm, 0, 0);
 
       reg_sweep (cm, 0);
 
       pvw->sse = gsl_matrix_get (cm, 0, 0);
+      gsl_matrix_free (cm);
 
       pvw->ssa = pvw->sst - pvw->sse;
 
-      pvw->n_groups = categoricals_total (cats);
-
-      pvw->mse = (pvw->sst - pvw->ssa) / (n - pvw->n_groups);
+      pvw->n_groups = categoricals_n_total (cats);
 
-      gsl_matrix_free (cm);
+      pvw->mse = (pvw->sst - pvw->ssa) / (pvw->n - pvw->n_groups);
     }
 
   for (v = 0; v < cmd->n_vars; ++v)
     {
       const struct categoricals *cats = covariance_get_categoricals (ws.vws[v].cov);
 
-      categoricals_done (cats);
-      
-      if (categoricals_total (cats) > ws.actual_number_of_groups)
-	ws.actual_number_of_groups = categoricals_total (cats);
+      if ( ! categoricals_is_complete (cats))
+	{
+	  continue;
+	}
+
+      if (categoricals_n_total (cats) > ws.actual_number_of_groups)
+	ws.actual_number_of_groups = categoricals_n_total (cats);
     }
 
   casereader_destroy (input);
@@ -509,18 +871,22 @@ run_oneway (const struct oneway_spec *cmd,
   taint_destroy (taint);
 
  finish:
+
   for (v = 0; v < cmd->n_vars; ++v)
     {
       covariance_destroy (ws.vws[v].cov);
+      levene_destroy (ws.vws[v].nl);
       dd_destroy (ws.dd_total[v]);
+      interaction_destroy (ws.vws[v].iact);
     }
+
   free (ws.vws);
   free (ws.dd_total);
-
 }
 
 static void show_contrast_coeffs (const struct oneway_spec *cmd, const struct oneway_workspace *ws);
 static void show_contrast_tests (const struct oneway_spec *cmd, const struct oneway_workspace *ws);
+static void show_comparisons (const struct oneway_spec *cmd, const struct oneway_workspace *ws, int depvar);
 
 static void
 output_oneway (const struct oneway_spec *cmd, struct oneway_workspace *ws)
@@ -529,7 +895,8 @@ output_oneway (const struct oneway_spec *cmd, struct oneway_workspace *ws)
 
   /* Check the sanity of the given contrast values */
   struct contrasts_node *coeff_list  = NULL;
-  ll_for_each (coeff_list, struct contrasts_node, ll, &cmd->contrast_list)
+  struct contrasts_node *coeff_next  = NULL;
+  ll_for_each_safe (coeff_list, coeff_next, struct contrasts_node, ll, &cmd->contrast_list)
     {
       struct coeff_node *cn = NULL;
       double sum = 0;
@@ -539,8 +906,11 @@ output_oneway (const struct oneway_spec *cmd, struct oneway_workspace *ws)
       if (ll_count (cl) != ws->actual_number_of_groups)
 	{
 	  msg (SW,
-	       _("Number of contrast coefficients must equal the number of groups"));
-	  coeff_list->bad_count = true;
+	       _("In contrast list %zu, the number of coefficients (%zu) does not equal the number of groups (%d). This contrast list will be ignored."),
+	       i, ll_count (cl), ws->actual_number_of_groups);
+
+	  ll_remove (&coeff_list->ll);
+	  destroy_coeff_list (coeff_list);
 	  continue;
 	}
 
@@ -564,6 +934,18 @@ output_oneway (const struct oneway_spec *cmd, struct oneway_workspace *ws)
       show_contrast_coeffs (cmd, ws);
       show_contrast_tests (cmd, ws);
     }
+
+  if ( cmd->posthoc )
+    {
+      int v;
+      for (v = 0 ; v < cmd->n_vars; ++v)
+	{
+	  const struct categoricals *cats = covariance_get_categoricals (ws->vws[v].cov);
+
+	  if ( categoricals_is_complete (cats))
+	    show_comparisons (cmd, ws, v);
+	}
+    }
 }
 
 
@@ -593,7 +975,7 @@ show_anova_table (const struct oneway_spec *cmd, const struct oneway_workspace *
   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
-  tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
+  tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
 
 
   for (i = 0; i < cmd->n_vars; ++i)
@@ -718,7 +1100,7 @@ show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace
       if ( v > 0)
 	tab_hline (t, TAL_1, 0, n_cols - 1, row);
 
-      for (count = 0; count < categoricals_total (cats); ++count)
+      for (count = 0; count < categoricals_n_total (cats); ++count)
 	{
 	  double T;
 	  double n, mean, variance;
@@ -726,8 +1108,8 @@ show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace
 
 	  struct string vstr;
 
-	  const union value *gval = categoricals_get_value_by_subscript (cats, count);
-	  const struct descriptive_data *dd = categoricals_get_user_data_by_subscript (cats, count);
+	  const struct ccase *gcc = categoricals_get_case_by_category (cats, count);
+	  const struct descriptive_data *dd = categoricals_get_user_data_by_category (cats, count);
 
 	  moments1_calculate (dd->mom, &n, &mean, &variance, NULL, NULL);
 
@@ -736,7 +1118,7 @@ show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace
 
 	  ds_init_empty (&vstr);
 
-	  var_append_value_name (cmd->indep_var, gval, &vstr);
+	  var_append_value_name (cmd->indep_var, case_data (gcc, cmd->indep_var), &vstr);
 
 	  tab_text (t, 1, row + count,
 		    TAB_LEFT | TAT_TITLE,
@@ -771,6 +1153,7 @@ show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace
 	  tab_double (t, 9, row + count, 0,  dd->maximum, fmt);
 	}
 
+      if (categoricals_is_complete (cats))
       {
 	double T;
 	double n, mean, variance;
@@ -802,12 +1185,13 @@ show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace
 	tab_double (t, 7, row + count, 0,
 		    mean + T * std_error, NULL);
 
+
 	/* Min and Max */
 	tab_double (t, 8, row + count, 0,  ws->dd_total[v]->minimum, fmt);
 	tab_double (t, 9, row + count, 0,  ws->dd_total[v]->maximum, fmt);
       }
 
-      row += categoricals_total (cats) + 1;
+      row += categoricals_n_total (cats) + 1;
     }
 
   tab_submit (t);
@@ -838,7 +1222,7 @@ show_homogeneity (const struct oneway_spec *cmd, const struct oneway_workspace *
   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Levene Statistic"));
   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df1"));
   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df2"));
-  tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
+  tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
 
   tab_title (t, _("Test of Homogeneity of Variances"));
 
@@ -846,7 +1230,7 @@ show_homogeneity (const struct oneway_spec *cmd, const struct oneway_workspace *
     {
       double n;
       const struct per_var_ws *pvw = &ws->vws[v];
-      double F = pvw->levene_w;
+      double F = levene_calculate (pvw->nl);
 
       const struct variable *var = cmd->vars[v];
       const char *s = var_to_string (var);
@@ -936,25 +1320,19 @@ show_contrast_coeffs (const struct oneway_spec *cmd, const struct oneway_workspa
 	   ++count, coeffi = ll_next (coeffi))
 	{
 	  const struct categoricals *cats = covariance_get_categoricals (cov);
-	  const union value *val = categoricals_get_value_by_subscript (cats, count);
+	  const struct ccase *gcc = categoricals_get_case_by_category (cats, count);
+	  struct coeff_node *coeffn = ll_data (coeffi, struct coeff_node, ll);
 	  struct string vstr;
 
 	  ds_init_empty (&vstr);
 
-	  var_append_value_name (cmd->indep_var, val, &vstr);
+	  var_append_value_name (cmd->indep_var, case_data (gcc, cmd->indep_var), &vstr);
 
 	  tab_text (t, count + 2, 1, TAB_CENTER | TAT_TITLE, ds_cstr (&vstr));
 
 	  ds_destroy (&vstr);
 
-	  if (cn->bad_count)
-	    tab_text (t, count + 2, c_num + 2, TAB_RIGHT, "?" );
-	  else
-	    {
-	      struct coeff_node *coeffn = ll_data (coeffi, struct coeff_node, ll);
-
-	      tab_text_format (t, count + 2, c_num + 2, TAB_RIGHT, "%g", coeffn->coeff);
-	    }
+	  tab_text_format (t, count + 2, c_num + 2, TAB_RIGHT, "%g", coeffn->coeff);
 	}
       ++c_num;
     }
@@ -1065,15 +1443,12 @@ show_contrast_tests (const struct oneway_spec *cmd, const struct oneway_workspac
                            (v * lines_per_variable) + i + 1 + n_contrasts,
                            TAB_CENTER | TAT_TITLE, "%d", i + 1);
 
-	  if (cn->bad_count)
-	    continue;
-
 	  for (coeffi = ll_head (&cn->coefficient_list);
 	       coeffi != ll_null (&cn->coefficient_list);
 	       ++ci, coeffi = ll_next (coeffi))
 	    {
 	      double n, mean, variance;
-	      const struct descriptive_data *dd = categoricals_get_user_data_by_subscript (cats, ci);
+	      const struct descriptive_data *dd = categoricals_get_user_data_by_category (cats, ci);
 	      struct coeff_node *cn = ll_data (coeffi, struct coeff_node, ll);
 	      const double coef = cn->coeff; 
 	      double winv ;
@@ -1165,3 +1540,127 @@ show_contrast_tests (const struct oneway_spec *cmd, const struct oneway_workspac
 }
 
 
+
+static void
+show_comparisons (const struct oneway_spec *cmd, const struct oneway_workspace *ws, int v)
+{
+  const int n_cols = 8;
+  const int heading_rows = 2;
+  const int heading_cols = 3;
+
+  int p;
+  int r = heading_rows ;
+
+  const struct per_var_ws *pvw = &ws->vws[v];
+  const struct categoricals *cat = pvw->cat;
+  const int n_rows = heading_rows + cmd->n_posthoc * pvw->n_groups * (pvw->n_groups - 1);
+
+  struct tab_table *t = tab_create (n_cols, n_rows);
+
+  tab_headers (t, heading_cols, 0, heading_rows, 0);
+
+  /* Put a frame around the entire box, and vertical lines inside */
+  tab_box (t,
+	   TAL_2, TAL_2,
+	   -1, -1,
+	   0, 0,
+	   n_cols - 1, n_rows - 1);
+
+  tab_box (t,
+	   -1, -1,
+	   -1, TAL_1,
+	   heading_cols, 0,
+	   n_cols - 1, n_rows - 1);
+
+  tab_vline (t, TAL_2, heading_cols, 0, n_rows - 1);
+
+  tab_title (t, _("Multiple Comparisons (%s)"), var_to_string (cmd->vars[v]));
+
+  tab_text_format (t,  1, 1, TAB_LEFT | TAT_TITLE, _("(I) %s"), var_to_string (cmd->indep_var));
+  tab_text_format (t,  2, 1, TAB_LEFT | TAT_TITLE, _("(J) %s"), var_to_string (cmd->indep_var));
+  tab_text (t,  3, 0, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
+  tab_text (t,  3, 1, TAB_CENTER | TAT_TITLE, _("(I - J)"));
+  tab_text (t,  4, 1, TAB_CENTER | TAT_TITLE, _("Std. Error"));
+  tab_text (t,  5, 1, TAB_CENTER | TAT_TITLE, _("Sig."));
+
+  tab_joint_text_format (t, 6, 0, 7, 0, TAB_CENTER | TAT_TITLE,
+                         _("%g%% Confidence Interval"),
+                         (1 - cmd->alpha) * 100.0);
+
+  tab_text (t,  6, 1, TAB_CENTER | TAT_TITLE, _("Lower Bound"));
+  tab_text (t,  7, 1, TAB_CENTER | TAT_TITLE, _("Upper Bound"));
+
+
+  for (p = 0; p < cmd->n_posthoc; ++p)
+    {
+      int i;
+      const struct posthoc *ph = &ph_tests[cmd->posthoc[p]];
+
+      tab_hline (t, TAL_2, 0, n_cols - 1, r);
+
+      tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, gettext (ph->label));
+
+      for (i = 0; i < pvw->n_groups ; ++i)
+	{
+	  double weight_i, mean_i, var_i;
+	  int rx = 0;
+	  struct string vstr;
+	  int j;
+	  struct descriptive_data *dd_i = categoricals_get_user_data_by_category (cat, i);
+	  const struct ccase *gcc = categoricals_get_case_by_category (cat, i);
+	  
+
+	  ds_init_empty (&vstr);
+	  var_append_value_name (cmd->indep_var, case_data (gcc, cmd->indep_var), &vstr);
+
+	  if ( i != 0)
+	    tab_hline (t, TAL_1, 1, n_cols - 1, r);
+	  tab_text (t, 1, r, TAB_LEFT | TAT_TITLE, ds_cstr (&vstr));
+
+	  moments1_calculate (dd_i->mom, &weight_i, &mean_i, &var_i, 0, 0);
+
+	  for (j = 0 ; j < pvw->n_groups; ++j)
+	    {
+	      double std_err;
+	      double weight_j, mean_j, var_j;
+	      double half_range;
+	      const struct ccase *cc;
+	      struct descriptive_data *dd_j = categoricals_get_user_data_by_category (cat, j);
+	      if (j == i)
+		continue;
+
+	      ds_clear (&vstr);
+	      cc = categoricals_get_case_by_category (cat, j);
+	      var_append_value_name (cmd->indep_var, case_data (cc, cmd->indep_var), &vstr);
+	      tab_text (t, 2, r + rx, TAB_LEFT | TAT_TITLE, ds_cstr (&vstr));
+
+	      moments1_calculate (dd_j->mom, &weight_j, &mean_j, &var_j, 0, 0);
+
+	      tab_double  (t, 3, r + rx, 0, mean_i - mean_j, 0);
+
+	      std_err = pvw->mse;
+	      std_err *= weight_i + weight_j;
+	      std_err /= weight_i * weight_j;
+	      std_err = sqrt (std_err);
+
+	      tab_double  (t, 4, r + rx, 0, std_err, 0);
+	  
+	      tab_double (t, 5, r + rx, 0, 2 * multiple_comparison_sig (std_err, pvw, dd_i, dd_j, ph), 0);
+
+	      half_range = mc_half_range (cmd, pvw, std_err, dd_i, dd_j, ph);
+
+	      tab_double (t, 6, r + rx, 0,
+			   (mean_i - mean_j) - half_range, 0 );
+
+	      tab_double (t, 7, r + rx, 0,
+			   (mean_i - mean_j) + half_range, 0 );
+
+	      rx++;
+	    }
+	  ds_destroy (&vstr);
+	  r += pvw->n_groups - 1;
+	}
+    }
+
+  tab_submit (t);
+}