####generate simulated data: nSNP markers, nSNP0 causal ones, n cases/ctls,
#### OR is the causal OR,
#### a haplotype is dichotomized from a MVN with AR-1(rho) corr,

simSNPs<-function(nSNP=40, nSNP0=4, OR=1, n=200, rho=0.8, 
                  MAFlow=0.1, MAFup=0.4){

q=nSNP

R<-matrix(0, nrow=q, ncol=q)
for(i in 1:q)
  for(j in 1:q)
    R[i, j]<-rho^(abs(i-j))
svd.R<-svd(R)
R1<-svd.R$u %*% diag(sqrt(svd.R$d))
#R2<- diag(sqrt(svd.R$d)) %*% t(svd.R$v)
#Note above: R1 %*% t(R1)= R
#R1<-svd.R$u %*% diag(1/sqrt(svd.R$d))
##Note above: R1 %*% t(R1) = Inv(R)
#R2<- diag(1/sqrt(svd.R$d)) %*% t(svd.R$v)

##background disease prev = 0.2
RR<-log(OR); b0<-log(0.2/0.8)

MAFs<-runif(q, MAFlow, MAFup)
cutoff<-qnorm(MAFs)

bs<-rep(0, q)
loc0s<-sample(1:q, nSNP0)
bs[loc0s]<-RR

X<-matrix(0, nrow=n+n, ncol=q)
Y<-rep(0, n+n); Y[(n+1):(2*n)]<-1
i<-1
#sampling controls:
while ( i <= n){
  X0<-rnorm(q, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1 %*% X0   #: X1 ~ MVN(0, R)
  X2<-ifelse(X1<cutoff, 1, 0)
  X0<-rnorm(q, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1 %*% X0   #: X1 ~ MVN(0, R)
  X3<-ifelse(X1<cutoff, 1, 0)
  X4<-X2+ X3
  pr<-1/(1 + exp(-(b0 + sum(bs * X4) )))
  Y1<-sample(c(0, 1), 1, prob=c(1-pr, pr))
  if (Y1==0){
    X[i, ]<-X4
    i<-i+1
    }
  }
#sampling cases:
while ( i <= 2*n){
  X0<-rnorm(q, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1 %*% X0   #: X1 ~ MVN(0, R)
  X2<-ifelse(X1<cutoff, 1, 0)
  X0<-rnorm(q, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1 %*% X0   #: X1 ~ MVN(0, R)
  X3<-ifelse(X1<cutoff, 1, 0)
  X4<-X2+ X3
  pr<-1/(1 + exp(-(b0 + sum(bs * X4) )))
  Y1<-sample(c(0, 1), 1, prob=c(1-pr, pr))
  if (Y1==1){
    X[i, ]<-X4
    i<-i+1
    }
  }

list(Y=Y, X=X)
}