########simulate SNPs from a latent multivariate Gaussian variable with
######## an AR1(rho) correlation structure. 
######## V2: the block of causal SNPs and the block of non-causal ones
########     are independent (while in V1, they all from AR-1). 
######## Wei Pan, weip@biostat.umn.edu, 10/25/10

########Input: 
#               ORs: association in OR between the causal SNPs and outcome;
#                   implicitly,  # of causal SNPs = length(ORs);
#               n0: # of controls;
#               n1: # of cases ;
#               nSNP: # of noise/marker SNPs to be generated;
#               rho: the parameter in AR1(rho) for the latent multivariate
#                    Gaussian variable to be discretized to SNPs;
#                    rho=0 means all SNPs are independent;
#               MAF0slow, MAF0sup: MAFs for the causal SNPs from 
#                                   Unif(MAF0slow, MAF0sup);
#               MAFslow, MAFsup: MAF's of the noise SNPs are drawn from
#                                 Unif(MAFslow, MAFsup);
#               p0: background disease prevalence; i.e. true logistic reg model's intercept=log(p0/(1-p0))
########Output: a list of the binary outcome Y (=0 or 1) and SNPs (=0, 1 or 2);
#               Y is a vector of length n=n1+n0; X is a matrix of n by nSNP;
#               and location indices for causal SNPs.

simAR1RareSNP2<-function(ORs, n0=500, n1=500, nSNP=0, rho=0,
                 MAF0slow=0.001, MAF0sup=0.01,
                 MAFslow=0.001, MAFsup=0.01, p0=0.05){

nSNP0<-length(ORs)
q<-nSNP+nSNP0

R0<-matrix(1, nrow=nSNP0, ncol=nSNP0)
for(i in 1:nSNP0)
  for(j in 1:nSNP0)
    if (i!=j) R0[i, j]<-rho^(abs(i-j))
svd.R0<-svd(R0)
R1<-svd.R0$u %*% diag(sqrt(svd.R0$d))

##background disease prev = p0
b0<-log(p0/(1-p0))
logORs<-log(ORs)

MAF0s<-runif(nSNP0, MAF0slow, MAF0sup)
MAFs<-runif(nSNP, MAFslow, MAFsup)
cutoff0<-qnorm(MAF0s)
cutoff<-qnorm(MAFs)

X<-matrix(0, nrow=n0+n1, ncol=nSNP0)
Y<-rep(0, n0+n1); Y[(n0+1):(n0+n1)]<-1
i<-1
#sampling controls:
while ( i <= n0){
  X0<-rnorm(nSNP0, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1 %*% X0   #: X1 ~ MVN(0, R)
  X2<-ifelse(X1<cutoff0, 1, 0)
  X0<-rnorm(nSNP0, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1 %*% X0   #: X1 ~ MVN(0, R)
  X3<-ifelse(X1<cutoff0, 1, 0)
  X4<-X2+ X3
  pr<-1/(1 + exp(-(b0 + sum(logORs * 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 <= n0+n1){
  X0<-rnorm(nSNP0, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1 %*% X0   #: X1 ~ MVN(0, R)
  X2<-ifelse(X1<cutoff0, 1, 0)
  X0<-rnorm(nSNP0, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1 %*% X0   #: X1 ~ MVN(0, R)
  X3<-ifelse(X1<cutoff0, 1, 0)
  X4<-X2+ X3
  pr<-1/(1 + exp(-(b0 + sum(logORs * X4))))
  Y1<-sample(c(0, 1), 1, prob=c(1-pr, pr))
  if (Y1==1){
    X[i, ]<-X4
    i<-i+1
    }
  }
if (nSNP>0){
###############generate noise SNPs:
R0noise<-matrix(1, nrow=nSNP, ncol=nSNP)
for(i in 1:nSNP)
  for(j in 1:nSNP)
    if (i!=j) R0noise[i, j]<-rho^(abs(i-j))
svd.R0noise<-svd(R0noise)
R1noise<-svd.R0noise$u %*% diag(sqrt(svd.R0noise$d))

Xnoise<-matrix(0, nrow=n0+n1, ncol=nSNP)
for(i in 1:(n0+n1)){
  X0<-rnorm(nSNP, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1noise %*% X0   #: X1 ~ MVN(0, R)
  X2<-ifelse(X1<cutoff, 1, 0)
  X0<-rnorm(nSNP, 0, 1) #: X0 ~ MVN(0, I)
  X1<-R1noise %*% X0   #: X1 ~ MVN(0, R)
  X3<-ifelse(X1<cutoff, 1, 0)
  X4<-X2+ X3
    Xnoise[i, ]<-X4
  }

Xall=cbind(X, Xnoise)
} else Xall<-X

###############index for causal SNPs; mix causal and noise SNPs:
############### SNP j is causal iff SNP0index[j]=1 or -1; =0 o/w
signORs<-rep(1, nSNP0)
signORs[ORs<1]<-0-1
indx<-sample(1:q, q, replace = FALSE)
SNP0indx<-c(signORs,  rep(0, nSNP))[indx]

list(Y=Y, X=Xall[,indx], SNP0indx=SNP0indx)

}