Gist for 80/20 post on Simply Stats

gistfile1.r

# Load library
library(RSkittleBrewer)

# Differences in mean between 0 and 5
sig = seq(0,5,length=100)

# This is the sample size
ss = 20

# Calculate power from 200 replicated tests
wiltest = ttest = ztest = stest = matrix(0,nrow=100,ncol=200)
for(i in seq(along=sig)){
  for(j in 1:200){
    x = rnorm(ss,mean=sig[i])
    wiltest[i,j] = wilcox.test(x,alternative="greater")$p.value < 0.05
    ttest[i,j] = t.test(x,alternative="greater")$p.value < 0.05
    ztest[i,j] = mean(x) > pnorm(0.95)/sqrt(ss)
    stest[i,j] = binom.test(c(sum(x < 0),sum(x > 0)), p=0.5,alternative="less")$p.value < 0.05
  }
  cat(i)
}

# Plot power ratio
plot(sig,rowMeans(ttest)/rowMeans(ztest),col=tropical[1],
     type="l",lwd=4,xlab="Mean difference",ylab="Ratio of Power Between Method and Z-test",ylim=c(0,1))
lines(sig,rowMeans(wiltest)/rowMeans(ztest),col=tropical[2],lwd=4)
lines(sig,rowMeans(stest)/rowMeans(ztest),col=tropical[3],lwd=4)
legend(2,0.9,col=tropical[1:3],lwd=rep(4,3),legend=c("T-test","Wilcox","Sign test"))

# Find signal size where power is greater than 80%
sig[which.max(rowMeans(wiltest)/rowMeans(ztest) > 0.8)]
sig[which.max(rowMeans(ttest)/rowMeans(ztest) > 0.8)]
sig[which.max(rowMeans(stest)/rowMeans(ztest) > 0.8)]