the influence of exposure's winner's curse on the outcome with different amounts of confounding

winners_curse_confounding.r

# no confounder

run <- function(nsnp, nid, ueff, nsim, prop)
{
	g <- matrix(rbinom(nsnp * nid, 2, 0.5), nid, nsnp)
	vg <- apply(g, 2, var)
	conf <- rnorm(nid)

	res1 <- array(0, nsim)
	res2 <- array(0, nsim)
	res3 <- array(0, nsim)
	res4 <- array(0, nsim)
	for(i in 1:nsim)
	{
		message(i)
		eff1 <- rnorm(nsnp)
		eff1[sample(nsnp * prop, 1:nsnp, replace=FALSE)] <- 0
		eff2 <- rnorm(nsnp)
		eff2[sample(nsnp * prop, 1:nsnp, replace=FALSE)] <- 0

		trait1 <- scale(g %*% eff1) + conf * sqrt(ueff)
		trait2 <- scale(g %*% eff2) + conf * sqrt(ueff)

		a <- cov(trait1, g) / vg
		b <- cov(trait2, g) / vg
		ind <- which.max(a)
		res1[i] <- a[ind]
		res2[i] <- b[ind]
		res3[i] <- mean(b, na.rm=T)
		res4[i] <- sd(b, na.rm=T)
	}

	return(data.frame(
		i=1:nsim,
		exp=res1,
		out=res2,
		mout=res3,
		sdout=res4
	))

}


a <- run(2000, 500, 0, 100, 1)
b <- run(2000, 500, 0.5, 100, 1)
c <- run(2000, 500, 1, 100, 1)
d <- run(2000, 500, 4, 100, 1)

t.test(a$out, rnorm(100, a$mout[1], a$sdout[1]))
t.test(b$out, rnorm(100, b$mout[1], b$sdout[1]))
t.test(c$out, rnorm(100, c$mout[1], c$sdout[1]))
t.test(d$out, rnorm(100, d$mout[1], d$sdout[1]))