Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
power graph
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.