require(OpenMx)

# Read in and setup data. 
dat			<- read.csv ("Module2_Medland_data.txt")
# Make a zygosity variable coded nv for MZ and .5 for DZ
dat$zyg		<-1/dat$zyg
# reformat data for openMx
dataMX		<- mxData( observed=dat, type="raw" )
# variable to analyse 
variables	<-c("Disinhibition") 
# which variables to analyse
analyse		<-c("Disinhibition1","Disinhibition2") 
# number of variables
nv			<- 1 

# Make a matrix with nv row and nv column to hold the estimate of the means
means		<- mxMatrix( "Full", 1, nv*2, free=T, label=variables, name="expMean" )

# Make matrices with nv row and nv column to hold the a, c, and e Path Coefficients
pathA		<-  mxMatrix( "Lower", nv, nv, free=T, values=.5, name="a" )
pathC		<-  mxMatrix( "Lower", nv, nv, free=T, values=.1, name="c" )
pathE		<-  mxMatrix( "Lower", nv, nv, free=T, values=.5, name="e" )

# Calculate the A, C, and E Variance Components
covA		<- mxAlgebra( a %*% t(a), name="A" )
covC		<- mxAlgebra( c %*% t(c), name="C" )
covE		<- mxAlgebra( e %*% t(e), name="E" )

# Specify the Variance/Covariance matrix
covariance	<- mxAlgebra( rbind( cbind(A+C+E, data.zyg %x% A+C),
                                 cbind(data.zyg %x% A+C, A+C+E)), name="expCov" )

# Setup the Model
funML		<- mxFitFunctionML() # specify the optimiser

expectation	<- mxExpectationNormal( covariance="expCov", means="expMean",
				dimnames=analyse ) # specify model expectations

ACEModel	<- mxModel( pathA, pathC, pathE, covA, covC, covE,  
				means, funML, covariance, dataMX, expectation, name="ACEmodel" ) # setup the model

# Run ACE model
ACEFit    	<- mxRun(ACEModel) # run the model

# Look at the output
(ACESumm   	<- summary(ACEFit))

#*****************Run AE model*****************#
# Respecify the model
pathC		<-  mxMatrix( "Lower", nv, nv, free=F, values=0, name="c" )
AEModel	<- mxModel( pathA, pathC, pathE, covA, covC, covE,  
				means, funML, covariance, dataMX, expectation, name="AEmodel" ) # setup the model
AEFit    	<- mxRun(AEModel) # run the model
(AESumm   	<- summary(AEFit))

#*****************Run CE model*****************#
# Respecify the model
pathA		<-  mxMatrix( "Lower", nv, nv, free=F, values=0, name="a" )
pathC		<-  mxMatrix( "Lower", nv, nv, free=T, values=.1, name="c" )
CEModel	<- mxModel( pathA, pathC, pathE, covA, covC, covE,  
				means, funML, covariance, dataMX, expectation, name="CEmodel" ) # setup the model
CEFit    	<- mxRun(CEModel) # run the model
(CESumm   	<- summary(CEFit))


#*****************Run E model*****************#
# Respecify the model
pathA		<-  mxMatrix( "Lower", nv, nv, free=F, values=0, name="a" )
pathC		<-  mxMatrix( "Lower", nv, nv, free=F, values=0, name="c" )
EModel	<- mxModel( pathA, pathC, pathE, covA, covC, covE,  
				means, funML, covariance, dataMX, expectation, name="Emodel" ) # setup the model
EFit    	<- mxRun(EModel) # run the model
(ESumm   	<- summary(EFit))