# -----------------------------------------------------------------------
# Program: sim_run.R  
#  Author: Hermine Maes and Benjamin Neale
#    Date: 03 05 2014 
#
# Univariate Twin Analysis model to test for the presence of C
# Creates a function that takes asq, csq, esq, nmz, and ndz
# Output is test of C
#
# -----------------------------------------------------------------------

require(OpenMx)
require(MASS)

sim_run<-function(asq,csq,esq,nmz,ndz)
{
# Defining covariance matrices for simulation
#------------------------------------------------------------------------
mzcov     <-matrix(c(asq+csq+esq,asq+csq,asq+csq,asq+csq+esq),2,2)
dzcov     <-matrix(c(asq+csq+esq,0.5*asq+csq,0.5*asq+csq,asq+csq+esq),2,2)

# Simulate Data using the mvrnorm command
# -----------------------------------------------------------------------
Vars      <- 'bmi'
nv        <- 1
selVars   <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") 
ntv       <- nv*2
mzData    <- data.frame(mvrnorm(nmz,mu=c(0,0),Sigma=mzcov))
colnames(mzData)<-(selVars)
dzData    <- data.frame(mvrnorm(ndz,mu=c(0,0),Sigma=dzcov))
colnames(dzData)<-(selVars)

# Fit ACE Model with RawData and Matrices Input
# -----------------------------------------------------------------------
# Matrices a, c, and e to store a, c, and e path coefficients
pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" )
pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="c11", name="c" )
pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" )
# Matrices A, C, and E compute variance components
covA      <- mxAlgebra( expression=a %*% t(a), name="A" )
covC      <- mxAlgebra( expression=c %*% t(c), name="C" )
covE      <- mxAlgebra( expression=e %*% t(e), name="E" )
# Algebra to compute total variances and standard deviations (diagonal only)
covP      <- mxAlgebra( expression=A+C+E, name="V" )
matI      <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
iSD       <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")
# Matrix & Algebra for expected means vector
meanG     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 0, label="mean", name="Mean" )
expMean   <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean")
# Algebra for expected variance/covariance matrix in MZ & DZ twins
expCovMZ  <- mxAlgebra( expression= rbind  ( cbind(V, A+C), cbind(A+C, V)), name="expCovMZ" )
expCovDZ  <- mxAlgebra( expression= rbind  ( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)),  name="expCovDZ" ) 

dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )
objMZ     <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars )
objDZ     <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars ) 

pars      <- list(pathA, pathC, pathE, covA, covC, covE, covP, matI, iSD, meanG)
modelMZ   <- mxModel("MZ", pars, expMean, expCovMZ, dataMZ, objMZ)
modelDZ   <- mxModel("DZ", pars, expMean, expCovDZ, dataDZ, objDZ)
minus2ll  <- mxAlgebra( expression=MZ.objective + DZ.objective, name="minus2loglikelihood" )
obj       <- mxAlgebraObjective("minus2loglikelihood")
twinACE_Model   <- mxModel("twinACE", pars, modelMZ, modelDZ, minus2ll, obj)
twinACE_Fit     <- mxRun(twinACE_Model)

# Fit AE model
# -----------------------------------------------------------------------
twinAE_Model <- mxModel(twinACE_Fit, name="twinAE")
twinAE_Model <- omxSetParameters( twinAE_Model, label="c11", free=FALSE, values=0 ) # drop c at 0

twinAE_Fit   <- mxRun(twinAE_Model)
twinAE_Fit@output$Minus-twinACE_Fit@output$Minus
}