#Simulate Dataset for MR pratical

#Income in 10s of thousands

library(MASS)

nind <- 10000

p <- 0.5
q <- 1 - p

a <- sqrt(0.001/(2*p*q))

rs3091244 <- sample(c(-a,0,a), nind, replace = TRUE, prob = c(p^2, 2*p*q, q^2))  

mu <- c(0, 0, 0, 0)
Sigma <- matrix(c(1, 0.4, -0.3, -0.3, 0.4, 1, -0.1, -0.3, -0.3, -0.3, 1, 0.1, -0.3, -0.3, 0.1, 1),4,4)

x <- mvrnorm(n = nind, mu, Sigma, tol = 1e-6, empirical = FALSE, EISPACK = FALSE)

CRP <- x[,1]
SBP <- x[,2]
INCOME <- x[,3]
HDL <- x[,4]

CRP <- 0.2*CRP + 2
SBP <- 10*SBP + 120
INCOME <- INCOME + 7.5
HDL <- 0.1*HDL + 1.2

CRP <- CRP + rs3091244

rs3091244 <- replace(rs3091244, rs3091244==a, 2)
rs3091244 <- replace(rs3091244, rs3091244==0, 1)
rs3091244 <- replace(rs3091244, rs3091244==-a, 0)

names <- c("rs3091244", "CRP", "SBP", "INCOME", "HDL")

d <- cbind(rs3091244, CRP, SBP, INCOME, HDL)


e <- data.frame(d)

write.table(e, file = "data.txt", row.names = FALSE, quote = FALSE)