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

# '#' denotes a comment

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

r<-0.5		#Correlation (r)

N<-1000		#Number of pairs of scores

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

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#													#
#				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) 
{ 
z<-rnorm(1,0,1)
if (z>0) pairs[i,2] <- NA
}
write.table(pairs,file=outfile,row.names=FALSE,col.names=FALSE,na=".")