rm(list=ls(all=TRUE))
#
#
# simulate multivariate normal data
################## input
# np and nv
#
#
x<-matrix(c(.5^.5,0,0.8,0.36^.5),2,2)
y<-matrix(c(.2^.5,0,0.5059644,0.3794733),2,2)
z<-matrix(c(.3^.5,0,0,.6^.5),2,2)
A<-(t(x))%*%x
D<-(t(y))%*%y
E<-(t(z))%*%z
sigmz<-cov2cor(matrix(c(rbind((cbind(A+D+E,A+D)),(cbind(A+D,A+D+E)))),4,4))
sigdz<-cov2cor(matrix(c(rbind((cbind(A+D+E,0.5*A+D)),(cbind(0.5*A+D,A+D+E)))),4,4))


nrep=250
nmz=2000*nrep
ndz=4000*nrep
nv=4
# define sigma the cov. matrix
sigma=sigmz
# means
mu=c(0,0,0,0)
##################
npnv<-nmz*nv
np=nmz
rdatmz<-matrix(rnorm(npnv),np,nv)
rdatmz<-rdatmz%*%chol(sigma)+(t(mu)%x%matrix(1,np,1))
##################
npnv<-ndz*nv
np=ndz
sigma=sigdz
rdatdz<-matrix(rnorm(npnv),np,nv)
rdatdz<-rdatdz%*%chol(sigma)+(t(mu)%x%matrix(1,np,1))
#####################################
write(t(rdatmz),file="rdat_ACE_mz",ncolumn=nv)
write(t(rdatdz),file="rdat_ACE_dz",ncolumn=nv)
######################
#concordant low + individual high (selection 1)
thresl=-1
thresh=2.35
select25mz<-matrix(c(rbinom(nmz,1,.25)),nmz,1)
select25dz<-matrix(c(rbinom(ndz,1,.25)),ndz,1)
rdatmzs1<-rdatmz
rdatdzs1<-rdatdz
for (i in 1:nmz)
{
if (
(rdatmz[i,1] < thresl & rdatmz[i,3]<thresl & select25mz[i,1]==1 )
| (rdatmz[i,1] > thresh) 
| (rdatmz[i,3] > thresh)
)
(rdatmzs1[i,2]<-rdatmz[i,2]) & (rdatmzs1[i,4]<-rdatmz[i,4])
else 
(rdatmzs1[i,2]<-99) & (rdatmzs1[i,4]<-99) 
}

for (i in 1:ndz)
{
if (
(rdatdz[i,1] < thresl & rdatdz[i,3]<thresl & select25dz[i,1]==1 )
| (rdatdz[i,1] > thresh) 
| (rdatdz[i,3] > thresh)
)
(rdatdzs1[i,2]<-rdatdz[i,2]) & (rdatdzs1[i,4]<-rdatdz[i,4])
else 
(rdatdzs1[i,2]<-99) & (rdatdzs1[i,4]<-99) 
}

write(t(rdatmzs1),file="rdat_ACE_mz",ncolumn=nv,append=T)
write(t(rdatdzs1),file="rdat_ACE_dz",ncolumn=nv,append=T)

######################
#extreme concordant (selection 2)

thresl=-1.57
thresh=1.57

rdatmzs2<-rdatmz
rdatdzs2<-rdatdz     

for (i in 1:nmz)
{
if (
(rdatmz[i,1] < thresl & rdatmz[i,3]<thresl)
| (rdatmz[i,1] > thresh & rdatmz[i,3]>thresh) 
)
(rdatmzs2[i,2]<-rdatmz[i,2]) & (rdatmzs2[i,4]<-rdatmz[i,4])
else 
(rdatmzs2[i,2]<-99) & (rdatmzs2[i,4]<-99) 
}

for (i in 1:ndz)
{
if (
(rdatdz[i,1] < thresl & rdatdz[i,3]<thresl)
| (rdatdz[i,1] > thresh & rdatdz[i,3]>thresh) 
)
(rdatdzs2[i,2]<-rdatdz[i,2]) & (rdatdzs2[i,4]<-rdatdz[i,4])
else 
(rdatdzs2[i,2]<-99) & (rdatdzs2[i,4]<-99) 
}

write(t(rdatmzs2),file="rdat_ACE_mz",ncolumn=nv,append=T)
write(t(rdatdzs2),file="rdat_ACE_dz",ncolumn=nv,append=T)

######################
#extreme discordant and concordant (selection 3)
thresl=-1.57
thresh=1.57

rdatmzs3<-rdatmz
rdatdzs3<-rdatdz

for (i in 1:nmz)
{
if (
(rdatmz[i,1] < thresl & rdatmz[i,3]<thresl)
| (rdatmz[i,1] > thresh & rdatmz[i,3]>thresh)
| (rdatmz[i,1] > thresh & rdatmz[i,3]<thresl)
| (rdatmz[i,1] < thresl & rdatmz[i,3]>thresh)
)
(rdatmzs3[i,2]<-rdatmz[i,2]) & (rdatmzs3[i,4]<-rdatmz[i,4])
else 
(rdatmzs3[i,2]<-99) & (rdatmzs3[i,4]<-99) 
}

for (i in 1:ndz)
{
if (
(rdatdz[i,1] < thresl & rdatdz[i,3]<thresl)
| (rdatdz[i,1] > thresh & rdatdz[i,3]>thresh)
| (rdatdz[i,1] > thresh & rdatdz[i,3]<thresl)
| (rdatdz[i,1] < thresl & rdatdz[i,3]>thresh)
)
(rdatdzs3[i,2]<-rdatdz[i,2]) & (rdatdzs3[i,4]<-rdatdz[i,4])
else 
(rdatdzs3[i,2]<-99) & (rdatdzs3[i,4]<-99) 
}

write(t(rdatmzs3),file="rdat_ACE_mz",ncolumn=nv,append=T)
write(t(rdatdzs3),file="rdat_ACE_dz",ncolumn=nv,append=T)

######################
#extreme discordant (selection 4)

thresl=-.67
thresh=.67

rdatmzs4<-rdatmz
rdatdzs4<-rdatdz

for (i in 1:nmz)
{
if ( 
(rdatmz[i,1] > thresh & rdatmz[i,3]<thresl)
| (rdatmz[i,1] < thresl & rdatmz[i,3]>thresh)
)
(rdatmzs4[i,2]<-rdatmz[i,2]) & (rdatmzs4[i,4]<-rdatmz[i,4])
else 
(rdatmzs4[i,2]<-99) & (rdatmzs4[i,4]<-99) 
}

for (i in 1:ndz)
{
if ( 
(rdatdz[i,1] > thresh & rdatdz[i,3]<thresl)
| (rdatdz[i,1] < thresl & rdatdz[i,3]>thresh)
)
(rdatdzs4[i,2]<-rdatdz[i,2]) & (rdatdzs4[i,4]<-rdatdz[i,4])
else 
(rdatdzs4[i,2]<-99) & (rdatdzs4[i,4]<-99) 
}

write(t(rdatmzs4),file="rdat_ACE_mz",ncolumn=nv,append=T)
write(t(rdatdzs4),file="rdat_ACE_dz",ncolumn=nv,append=T)

######################
# individual selection (selection 5)

thresl=-2.375
thresh=2.375

rdatmzs5<-rdatmz
rdatdzs5<-rdatdz

for (i in 1:nmz)
{
if ( 
(rdatmz[i,1] > thresh)
| (rdatmz[i,1] < thresl)
| (rdatmz [i,3] > thresh)
| (rdatmz [i,3] < thresl)
)
(rdatmzs5[i,2]<-rdatmz[i,2]) & (rdatmzs5[i,4]<-rdatmz[i,4])
else 
(rdatmzs5[i,2]<-99) & (rdatmzs5[i,4]<-99) 
}

for (i in 1:ndz)
{
if ( 
(rdatdz[i,1] > thresh)
| (rdatdz[i,1] < thresl)
| (rdatdz [i,3] > thresh)
| (rdatdz [i,3] < thresl)
)
(rdatdzs5[i,2]<-rdatdz[i,2]) & (rdatdzs5[i,4]<-rdatdz[i,4])
else 
(rdatdzs5[i,2]<-99) & (rdatdzs5[i,4]<-99) 
}

write(t(rdatmzs5),file="rdat_ACE_mz",ncolumn=nv,append=T)
write(t(rdatdzs5),file="rdat_ACE_dz",ncolumn=nv,append=T)

######################
# MCAR (selection 6)

rdatmzs6<-rdatmz
rdatdzs6<-rdatdz

select25mz<-matrix(c(rbinom(nmz,1,.033)),nmz,1)
select25dz<-matrix(c(rbinom(ndz,1,.033)),ndz,1)

for (i in 1:nmz)
{
if (select25mz[i,]==1)
(rdatmzs6[i,2]<-rdatmz[i,2]) & (rdatmzs6[i,4]<-rdatmz[i,4])
else 
(rdatmzs6[i,2]<-99) & (rdatmzs6[i,4]<-99) 
}

for (i in 1:ndz)
{
if (select25dz[i,]==1)
(rdatdzs6[i,2]<-rdatdz[i,2]) & (rdatdzs6[i,4]<-rdatdz[i,4])
else 
(rdatdzs6[i,2]<-99) & (rdatdzs6[i,4]<-99) 
}

write(t(rdatmzs6),file="rdat_ACE_mz",ncolumn=nv,append=T)
write(t(rdatdzs6),file="rdat_ACE_dz",ncolumn=nv,append=T)

######################
# optimal NMZ/NDZ (selection 7)

rdatmzs7<-rdatmz
rdatdzs7<-rdatdz

selectmz<-matrix(c(rbinom(nmz,1,.039)),nmz,1)
selectdz<-matrix(c(rbinom(ndz,1,.0305)),ndz,1)

for (i in 1:nmz)
{
if (selectmz[i,]==1)
(rdatmzs7[i,2]<-rdatmz[i,2]) & (rdatmzs7[i,4]<-rdatmz[i,4])
else 
(rdatmzs7[i,2]<-99) & (rdatmzs7[i,4]<-99) 
}

for (i in 1:ndz)
{
if (selectdz[i,]==1)
(rdatdzs7[i,2]<-rdatdz[i,2]) & (rdatdzs7[i,4]<-rdatdz[i,4])
else 
(rdatdzs7[i,2]<-99) & (rdatdzs7[i,4]<-99) 
}

write(t(rdatmzs7),file="rdat_ACE_mz",ncolumn=nv,append=T)
write(t(rdatdzs7),file="rdat_ACE_dz",ncolumn=nv,append=T)