Games-Howell Post Hoc Test in R

Question

I am doing some data analysis for my masters and I had some data that is normally distributed but does not fit the assumption of homogeneity of variance and has unequal sample sizes. I have been doing some research and have found that the games-howell post hoc test can deal with this type of data but can not find any code/ algorithm for it in R. Does any one know of the code or a different post hoc test that will do the same thing but is supported in R? Thanks

Answer 1

I had the same issue, and I edited the function. It now provides Games-Howell by default, and the output is a bit clearer. I'll put it in the next version of the 'userfriendlyscience' package. Until then, here you go:

posthoc.tgh <- function(y, x, method=c("games-howell", "tukey"), digits=2) {
  ### Based on http://www.psych.yorku.ca/cribbie/6130/games_howell.R
  method <- tolower(method);
  tryCatch(method <- match.arg(method), error=function(err) {
    stop("Argument for 'method' not valid!");
  });

  res <- list(input = list(x=x, y=y, method=method, digits=digits));

  res$intermediate <- list(x = factor(x[complete.cases(x,y)]),
                               y = y[complete.cases(x,y)]);
      res$intermediate$n <- tapply(y, x, length);
      res$intermediate$groups <- length(res$intermediate$n);
      res$intermediate$df <- sum(res$intermediate$n) - res$intermediate$groups;
      res$intermediate$means <- tapply(y, x, mean);
      res$intermediate$variances <- tapply(y, x, var);

  res$intermediate$pairNames <- combn(levels(res$intermediate$x),
                                      2, paste0, collapse=":");

  res$intermediate$descriptives <- cbind(res$intermediate$n,
                                         res$intermediate$means,
                                         res$intermediate$variances);
  rownames(res$intermediate$descriptives) <- levels(res$intermediate$x);
  colnames(res$intermediate$descriptives) <- c('n', 'means', 'variances');

  ### Start on Tukey
  res$intermediate$errorVariance <-
    sum((res$intermediate$n-1) * res$intermediate$variances) /
    res$intermediate$df;
  res$intermediate$t <- combn(res$intermediate$groups, 2, function(ij) {
    abs(diff(res$intermediate$means[ij]))/
      sqrt(res$intermediate$errorVariance*sum(1/res$intermediate$n[ij]));
  } );
  res$intermediate$p.tukey <- ptukey(res$intermediate$t*sqrt(2),
                                     res$intermediate$groups,
                                     res$intermediate$df,
                                     lower.tail=FALSE);
  res$output <- list();
      res$output$tukey <- cbind(res$intermediate$t,
                                res$intermediate$df,
                                res$intermediate$p.tukey)                                     
      rownames(res$output$tukey) <- res$intermediate$pairNames;
      colnames(res$output$tukey) <- c('t', 'df', 'p');

  ### Start on Games-Howell
  res$intermediate$df.corrected <- combn(res$intermediate$groups, 2, function(ij) {               
    sum(res$intermediate$variances[ij] /
          res$intermediate$n[ij])^2 / 
      sum((res$intermediate$variances[ij] /
             res$intermediate$n[ij])^2 / 
            (res$intermediate$n[ij]-1));
  } );
  res$intermediate$t.corrected <- combn(res$intermediate$groups, 2, function(ij) {               
    abs(diff(res$intermediate$means[ij]))/
      sqrt(sum(res$intermediate$variances[ij] /
                 res$intermediate$n[ij]));
  } );    
  res$intermediate$p.gameshowell <- ptukey(res$intermediate$t.corrected*sqrt(2),
                                           res$intermediate$groups,
                                           res$intermediate$df.corrected,
                                           lower.tail=FALSE)  
  res$output$games.howell <- cbind(res$intermediate$t.corrected,
                                   res$intermediate$df.corrected,
                                   res$intermediate$p.gameshowell);
  rownames(res$output$games.howell) <- res$intermediate$pairNames;
  colnames(res$output$games.howell) <- c('t', 'df', 'p');

  ### Set class and return object
  class(res) <- 'posthocTukeyGamesHowell';
  return(res);

}

print.posthocTukeyGamesHowell <- function(x, digits=x$input$digits, ...) {
  print(x$intermediate$descriptives, digits=digits);
  cat('\n');
  if (x$input$method == 'tukey') {
    print(x$output$tukey);
  }
  else if (x$input$method == 'games-howell') {
    print(x$output$games.howell, digits=digits);
  }
}

An example:

> posthoc.tgh(y=diamonds$y, x=diamonds$cut);
              n means variances
Fair       1610   6.2      0.91
Good       4906   5.9      1.11
Very Good 12082   5.8      1.22
Premium   13791   5.9      1.59
Ideal     21551   5.5      1.15

                     t    df       p
Fair:Good         11.8  2985 0.0e+00
Fair:Very Good    16.0  2221 5.9e-11
Fair:Premium       9.1  2316 1.6e-11
Fair:Ideal        26.6  1925 3.1e-11
Good:Very Good     4.5  9497 7.7e-05
Good:Premium       5.1 10243 3.4e-06
Good:Ideal        19.8  7419 1.3e-11
Very Good:Premium 11.9 25871 0.0e+00
Very Good:Ideal   20.1 24473 0.0e+00
Premium:Ideal     32.7 26011 0.0e+00