Dunnett's test using means, standard deviations, and sample sizes

dunnet.R

###########
# Calculate a Dunnett's test
# using information from a meta-analysis
# using definitions at http://davidmlane.com/hyperstat/B112114.html
#
# Jarrett Byrnes
# 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 agains
# first value in the vector
# x = vector of means, n = vector of sample size
# and standard formula for MSE from ANOVA
dunnet.test <- function(x, sdVec, n){
  mse_test <- MSE(sdVec, n)
  df_test <- sum(n) - length(n) 
  
  ret <- sapply(2:length(x), function(i){
    nmean <- harmonic.mean(c(n[i], n[1]))
    tdun <- (x[i] - x[1])/sqrt(2*mse_test/nmean)
    c(mean = x[i], diff = x[i] - x[1], t = tdun, p=2*pt(abs(tdun), df_test, lower.tail=F))
  })
  
  return(t(ret))
}
 
#example
dunnet.test(c(3,2,8,1,6,1,3), c(2,2,2,0.5,5,4,4), n=c(5,5,5,5,4,3,5))