Produce statistical power contour plot for simple nested design w/ participants (pl) nested in labs (l)

nested_contour.R

# 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()