s.ratti gwas power calculations

sratti_power.R

library(pwr)
library(ggplot2)

# calculate the minumum r2 that could be detected

r2 <- pwr.r.test(n=c(100,200), power=0.8, sig.level=0.05/100000)$r^2

# 

dat <- expand.grid(
	n = seq(100, 2000, by=20),
	maf = c(0.05, 0.1, 0.15, 0.3, 0.5),
	nsnp = c(10000, 50000, 100000)
)

for(i in 1:nrow(dat))
{
	dat$r2[i] <- pwr.r.test(n=dat$n[i], power=0.8, sig.level=0.05 / dat$nsnp[i])$r^2
}

dat$eff <- dat$r2 / (2 * dat$maf - 2 * dat$maf^2)



ggplot(dat, aes(y=eff, x=n)) +
geom_line(aes(colour=as.factor(maf))) +
facet_grid(. ~ nsnp) +
scale_colour_brewer(type="div") +
labs(x="Sample size", y="Minimum detectable effect size at 80% power", colour="Minor allele\nfrequency") +
scale_x_continuous(breaks=seq(0,2000,by=100)) +
scale_y_continuous(breaks=seq(0,1,by=0.1), lim=c(0,1)) +
theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))
ggsave(file="~/Desktop/sratti_power.pdf")