library(boot)

cor.lik=cor.p=cor.chisq=0 # counters for how many time each method identifies X1 correctly.
for (i in 1:10000){

    x1 = rnorm(1000) # start from normal x's, where correlation is easy to express
    x2 = x1*sqrt(7/8) + rnorm(1000,sd=sqrt(1/8)) # can play with the numbers to control level of dependence

    X1 = (x1 > -0.5) + (x1>0.5) # discretize X's to {0,1,2}
    X2 = (x2 > -0.5) + (x2>0.5)
    
    rm(x1,x2) # use X1,X2 in modeling!
    
    p = inv.logit (0.3+0.5*(X1==2)) # recessive mode
    y = rbinom(size=1,n=length(p),p=p)

######### now model with three approaches and compare results

}