#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#			Simulation program for twin pairs in R				#
#				Michael Neale stolen from Benjamin Neale			#
#				March 2006								#
#     Generates conditionally missing values for second variable in list	#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

# '#' denotes a comment

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#													#
#				Parameters								#
#													#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

r<-0.5		      #Correlation (r)

N<-1000		      #Number of pairs of scores

outfile<-"selonx.rec"	#output file for the simulation

for (j in 1:10)
{
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#													#
#				Data generation							#
#													#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

vals<-matrix(rnorm((N*3)),,3);	#3 random normals per pair of observed scores

rootr<-sqrt(r);				#Square root of r
root1mr<-sqrt(1-r);             #Square root of (1-r)

#Data generation follows: 

pairs<-matrix(cbind(rootr*vals[,1]+root1mr*vals[,2],rootr*vals[,1]+root1mr*vals[,3]),,2)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# pairs now contains pairs of scores drawn from population with               #
# unit variances and correlation r							 	#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Now to make some of the scores missing, using a loop
for (i in 1:N) 
{ 
if (pairs[i,1]>0) pairs[i,2] <- NA
}
write.table(pairs,file=outfile,row.names=FALSE,na=".",col.names=FALSE)
# now we run Mx and capture the chisq to a file
system("c:\\mx\\lahey\\mxe.exe selonx.mx")
write(system("c:\\cygwin\\bin\\grep.exe Chi c:\\Mx\\r\\selonx.mxo",intern=TRUE),file="chis.txt",append=TRUE)}