dunnet.R

###########
# Calculate a Dunnett's test
# using information from a meta-analysis
# using definitions at http://davidmlane.com/hyperstat/B112114.html
#
# Jarrett Byrnes & Jon Lefcheck
# 12/8/2013
###########

#helper functions
harmonic.mean <- function(x) length(x)/sum(1/x)
grandMean <- function(x, n) weighted.mean(x, n)
MST <- function(x, n){
  gm <- grandMean(x,n)
  sst <- sum((x-gm)^2)
  mst <- sse/(length(x)-1)
  mst
}
MSE <- function(sdVec, nVec){
  sse <- sum(sdVec^2/(nVec-1))
  sse/(sum(nVec)-length(nVec))
}
# Function, where we compare values against
# first value in the vector (aka, the "control")
# x = vector of means, n = vector of sample size
# and standard formula for MSE from ANOVA
dunnet.test <- function(x, sdVec, n, one.sided = T, equal.var = T){
  mse_test <- MSE(sdVec, n)
  df_test <- (sum(n)-1) - (length(n)+1)
  ret <- sapply(2:length(x), function(i){
    nmean <- harmonic.mean(c(n[i], n[1]))
    if(one.sided == T & equal.var == T) { #One-sided test if means have equal variances
      #Calculate test statistic
      tdun <- (x[i] - x[1])/sqrt(2*mse_test/nmean)
      #Define correlation matrix (corr) with rho=0.5 (Dunnett 1955, 1964)
      corr <- matrix(0.5, length(x) - 1, length(x) - 1); diag(corr) <- 1
      c(mean = x[i], diff = x[i] - x[1], t = tdun, 
        p = pmvt(lower = rep(-Inf, length(x) - 1), upper = rep(tdun, length(x) - 1), df=df_test, corr = corr) )
    } else 
      if(one.sided == F & equal.var == T) {
        tdun <- (x[i] - x[1])/sqrt(2*mse_test/nmean)
        corr <- matrix(0.5, length(x) - 1, length(x) - 1); diag(corr) <- 1
        c(mean = x[i], diff = x[i] - x[1], t = tdun, 
          p = pmvt(lower = rep(tdun, length(x) - 1), upper = rep(Inf, length(x) - 1), df=df_test, corr = corr) )
    } else   
      if(one.sided == T & equal.var == F) { #If means do not have equal variances
        tdun <- (x[i] - x[1])/sqrt((sdVec[i]/n[i])+sdVec[1]/n[1])
        corr <- matrix(0.5, length(x) - 1, length(x) - 1); diag(corr) <- 1
        c(mean = x[i], diff = x[i] - x[1], t = tdun, 
          p = pmvt(lower = rep(-Inf, length(x) - 1), upper = rep(tdun, length(x) - 1), df=df_test, corr = corr) )
    } else {
        tdun <- (x[i] - x[1])/sqrt((sdVec[i]/n[i])+sdVec[1]/n[1])
        corr <- matrix(0.5, length(x) - 1, length(x) - 1); diag(corr) <- 1
        c(mean = x[i], diff = x[i] - x[1], t = tdun, 
          p = pmvt(lower = rep(tdun, length(x) - 1), upper = rep(tdun, length(x) - 1), df=df_test, corr = corr) )
    }
  } )
  return(t(ret))
}
        
#example
dunnet.test(x=c(3,2,8,1,6,1,3), sdVec=c(2,2,2,0.5,5,4,4), n=c(5,5,5,5,4,3,5), one.sided=T, equal.var=T)
dunnet.test(x=c(3,2,8,1,6,1,3), sdVec=c(2,2,2,0.5,5,4,4), n=c(5,5,5,5,4,3,5), one.sided=F, equal.var=T)
dunnet.test(x=c(3,2,8,1,6,1,3), sdVec=c(2,2,2,0.5,5,4,4), n=c(5,5,5,5,4,3,5), one.sided=T, equal.var=F)
dunnet.test(x=c(3,2,8,1,6,1,3), sdVec=c(2,2,2,0.5,5,4,4), n=c(5,5,5,5,4,3,5), one.sided=F, equal.var=F)