# -----------------------------------------------------------------------
# Program: sim_run.R  
#  Author: Hermine Maes and Benjamin Neale
#    Date: 03 10 2010 
#
# 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)
source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")
#source("GenEpiHelperFunctions.R")


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)


# Print Descriptive Statistics
# -----------------------------------------------------------------------
#summary(mzData)
#colMeans(mzData,na.rm=TRUE)
#cov(mzData,use="complete")
#summary(dzData)
#colMeans(dzData,na.rm=TRUE)
#cov(dzData,use="complete")

# Fit ACE Model with RawData and Matrices Input
# -----------------------------------------------------------------------
univACEModel <- mxModel("univACE",
    mxModel("ACE",
    # Matrices a, c, and e to store a, c, and e path coefficients
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" ), #X
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="c11", name="c" ), #Y
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" ), #Z
    # Matrices A, C, and E compute variance components
        mxAlgebra( expression=a %*% t(a), name="A" ),
        mxAlgebra( expression=c %*% t(c), name="C" ),
        mxAlgebra( expression=e %*% t(e), name="E" ),
    # Algebra to compute total variances and standard deviations (diagonal only)
        mxAlgebra( expression=A+C+E, name="V" ),
        mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
        mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"),
   ## Note that the rest of the mxModel statements do not change for bivariate/multivariate case
    # Matrix & Algebra for expected means vector
        mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 0, label="mean", name="Mean" ),
        mxAlgebra( expression= cbind(Mean,Mean), name="expMean"),
    # Algebra for expected variance/covariance matrix in MZ
        mxAlgebra( expression= rbind  ( cbind(A+C+E , A+C),
                                        cbind(A+C   , A+C+E)), name="expCovMZ" ),
    # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
        mxAlgebra( expression= rbind  ( cbind(A+C+E     , 0.5%x%A+C),
                                        cbind(0.5%x%A+C , A+C+E)),  name="expCovDZ" ) 
    ),
    mxModel("MZ",
        mxData( observed=mzData, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars )
    ),
    mxModel("DZ", 
        mxData( observed=dzData, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars ) 
    ),
    mxAlgebra( expression=MZ.objective + DZ.objective, name="-2sumll" ),
    mxAlgebraObjective("-2sumll")
)

univACEFit <- mxRun(univACEModel)
univACESumm <- summary(univACEFit)
#univACESumm

# Generate ACE Output - we don't generate this as we're just simulating
# -----------------------------------------------------------------------
#parameterSpecifications(univACEFit)
#expectedMeansCovariances(univACEFit)
#tableFitStatistics(univACEFit)

# Generate Table of Parameter Estimates using mxEval
#pathEstimatesACE <- print(round(mxEval(cbind(ACE.a,ACE.c,ACE.e), univACEFit),4))
#varComponentsACE <- print(round(mxEval(cbind(ACE.A/ACE.V,ACE.C/ACE.V,ACE.E/ACE.V), univACEFit),4))
#    rownames(pathEstimatesACE) <- 'pathEstimates'
#    colnames(pathEstimatesACE) <- c('a','c','e')
#    rownames(varComponentsACE) <- 'varComponents'
#    colnames(varComponentsACE) <- c('a^2','c^2','e^2')
#pathEstimatesACE
#varComponentsACE

# Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices
#ACEpathMatrices <- c("ACE.a","ACE.c","ACE.e","ACE.iSD","ACE.iSD %*% ACE.a","ACE.iSD %*% ACE.c","ACE.iSD %*% ACE.e")
#ACEpathLabels <- c("pathEst_a","pathEst_c","pathEst_e","isd","stPathEst_a","stPathEst_c","stPathEst_e")
#formatOutputMatrices(univACEFit,ACEpathMatrices,ACEpathLabels,Vars,4)

#ACEcovMatrices <- c("ACE.A","ACE.C","ACE.E","ACE.V","ACE.A/ACE.V","ACE.C/ACE.V","ACE.E/ACE.V")
#ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E")
#formatOutputMatrices(univACEFit,ACEcovMatrices,ACEcovLabels,Vars,4)

# Fit AE model
# -----------------------------------------------------------------------
univAEModel <- mxModel(univACEFit, name="univAE",
    mxModel(univACEFit$ACE,
        mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) # drop c at 0
    )
)
univAEFit <- mxRun(univAEModel)
univAESumm <- summary(univAEFit)

# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
univACENested <- list(univAEFit)
#tableFitStatistics(univACEFit,univACENested)

univAESumm$Mi-univACESumm$Mi
}