Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Produce statistical power contour plot for simple nested design w/ participants (pl) nested in labs (l)
# replace this with a path on your own machine
path <- "/Users/Jake/Desktop/nested_power.png"
# define function to compute statistical power
pow <- function(d, E, LC, l, pl){
t <- d/2/sqrt(E/(2^pl)/(2^l) + LC/(2^l))
DF <- (2^l) - 1
return(pt(qt(.975, DF), DF, ncp=t, lower.tail=F) +
pt(qt(.025, DF), DF, ncp=t))
}
# define parameters to compute power for
lRange <- seq(log2(4), log2(64), length.out=50)
plRange <- seq(log2(8), log2(256), length.out=50)
params <- expand.grid(d = c(.2, .4, .6),
LC = c(.01, .05, .1),
l = lRange,
pl = plRange)
# plot!
png(path, units="in", height=6, width=6, res=200, pointsize=13)
layout(matrix(1:9, nrow=3))
par(mar=c(3, 3, .5, .5)+.1)
par(oma=c(0,0,2,2))
by(params, list(params$d, params$LC), function(dat){
values <- matrix(mapply(pow, l=dat$l, pl=dat$pl,
MoreArgs=list(E=.5, d=unique(dat$d),
LC=unique(dat$LC))),
ncol = 50, nrow = 50, byrow=T)
contour(x=plRange, y=lRange, z=values, xaxt="n", yaxt="n", mgp=2:0,
xlab="Participants per Lab", ylab="Number of Labs",
levels=c(seq(.1, .9, .1), .99))
axis(side=1, at=pretty(plRange), labels=2^pretty(plRange))
axis(side=2, at=pretty(lRange), labels=2^pretty(lRange))
})
mtext(side=3, at=c(.2,.53,.87), outer=TRUE, line=.5, cex=1,
text=paste("Lab variance =",c("1%","5%","10%")))
mtext(side=4, at=c(.2,.53,.87), outer=TRUE, line=.5, cex=1,
text=paste("Cohen's d =",c(".6",".4",".2")))
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment