power graph

power.R

M = 1000
N = round(sort(rpois(M, 10)^1.5, decreasing = TRUE)) + 20
pows = c(.5,.8,.95)

setup = expand.grid(N=N, pow = pows)

es = mapply(FUN = function(n, power, alpha) power.t.test(n = n, sig.level = alpha, power = power)$delta,
       n = setup$N, power = setup$pow, MoreArgs = list(alpha = .05))

dim(es) = c(setup$es,length(N), length(pows))

par(yaxs = "i", xaxs = "i", las = 1)

plot(0,0,ty='n',ylim = c(0,100), xlim = c(min(es),max(es)*1.1),
     ylab = "% designs with at most specified power", xlab = "Hypothetical effect size")

for(i in 1:length(pows)){
  lines(es[,i], 100*(1:M)/M, lwd = 2, col = i)
  text(max(es[,i]), 100, paste0(pows[i] * 100,"% power"), 3, 
       col = i, xpd = TRUE, adj = .5)
}
axis(1)
axis(2)
abline(h = 50, col = "gray", lty = 2)

quantile(es[,1],p=.5)
mean(es[,3]<.6)