EWAS power simulations

ewas_power.r

library(devtools)
install_github("explodecomputer/simulateGP")
library(simulateGP)
library(ggplot2)

param <- rbind(
	expand.grid(
		ncase=100,
		ncontrol=100,
		prev=0.001,
		eff=c(0, 0.01, 0.05, 0.1, 0.15, 0.2),
		sim=1:500,
		or=NA,
		pval=NA
	),
	expand.grid(
		ncase=200,
		ncontrol=200,
		prev=0.001,
		eff=c(0, 0.01, 0.05, 0.1, 0.15, 0.2),
		sim=1:500,
		or=NA,
		pval=NA
	)
)


for(i in 1:nrow(param))
{
	message(i, " of ", nrow(param))
	cpg <- rbeta(400000, 200, 5) %>% matrix
	a <- risk_simulation(cpg, sqrt(param$eff[i]), param$prev[i], 1)
	a$C <- cpg
	b <- rbind(
		subset(a, disease == 0) %>% sample_n(param$ncontrol[i]),
		subset(a, disease == 1) %>% sample_n(param$ncase[i])
	)
	b$C <- scale(b$C)
	o <- summary(glm(disease ~ C, b, family="binomial"))$coefficients
	param$or[i] <- exp(o[2,1])
	param$pval[i] <- o[2,4]
}

save(param, file="power_simulations.rdata")

power_summary <- group_by(param, eff, ncase, ncontrol, prev) %>%
summarise(or=mean(or), pow=sum(pval < 0.05/450000)/n())

ggplot(power_summary, aes(x=or, y=pow)) +
geom_point(aes(colour=paste(ncase, "/", ncontrol))) +
geom_line(aes(colour=paste(ncase, "/", ncontrol))) +
labs(x="Odds ratio", y="Power", colour="Case/control")
ggsave("power_simulations.pdf")



## Convert OR to standardised mean difference
## https://www.meta-analysis.com/downloads/Meta-analysis%20Converting%20among%20effect%20sizes.pdf
## Calculate t-test power

p2 <- expand.grid(
	n = c(100, 200),
	or = c(1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5),
	d = NA,
	pow = NA
)
p2$d <- log(p2$or) * sqrt(3) / pi

for(i in 1:nrow(p2))
{
	message(i)
	o <- pwr.t.test(n=p2$n[i], d=p2$d[i], sig.level=0.05/450000, type='one.sample')
	p2$pow[i] <- o$power
}

ggplot(p2, aes(x=or, y=pow)) +
geom_point(aes(colour=paste(n, "/", n))) +
geom_line(aes(colour=paste(n, "/", n))) +
labs(x="Odds ratio", y="Power", colour="Case/control")
ggsave("power_simulations.pdf")