Simulation of stratifying instrument

stratify.R

# Simulate the smoking thing
# SNP influences smoking status
# Smoking status influences income
# Test for association of SNP on income in smokers and non smokers


nsim <- 100
res <- matrix(0, nsim, 4)
n <- 100000
for(i in 1:nsim)
{
	cat(i, "\n")
	g <- rbinom(n, 2, 0.5)

	b <- scale(g) * sqrt(0.1) + rbeta(n, 1, 1)
	b <- round(b)
	b1 <- scale(b) * sqrt(0.1) + rnorm(n)
	d <- data.frame(g, b, b1)
	res[i, 1] <- coefficients(summary(lm(b ~ g)))[2,4]
	res[i, 2] <- coefficients(summary(lm(b1 ~ g)))[2,4]
	res[i, 3] <- coefficients(summary(lm(b1 ~ g, subset(d, b==0))))[2,4]
	res[i, 4] <- coefficients(summary(lm(b1 ~ g, subset(d, b==1))))[2,4]

}

# Is there an inflation of test statistics when stratifying?
mean(res[,4])
median(res[,4])

mean(res[,3])
median(res[,4])

hist(res[,3])
hist(res[,4])

# No inflation, all good.