# Heterogeneity ACE Model 
# -----------------------------------------------------------------------
univHetBinACEModel <- mxModel("univHetBinACE",
    mxModel("ACE",
    # Matrices a, c, and e to store a, c, and e path coefficients
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="am11", name="am" ),
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="cm11", name="cm" ),
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="em11", name="em" ),
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="af11", name="af" ),
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="cf11", name="cf" ),
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="ef11", name="ef" ),
        mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1, label="rg11", name="rg"),
    # Matrices A, C, and E compute variance components
        mxAlgebra( expression=am %*% t(am), name="Am" ),
        mxAlgebra( expression=cm %*% t(cm), name="Cm" ),
        mxAlgebra( expression=em %*% t(em), name="Em" ),
        mxAlgebra( expression=af %*% t(af), name="Af" ),
        mxAlgebra( expression=cf %*% t(cf), name="Cf" ),
        mxAlgebra( expression=ef %*% t(ef), name="Ef" ),
        mxAlgebra( expression=cbind(Am/Vm,Cm/Vm,Em/Vm), name="ACE2m"),
        mxAlgebra( expression=cbind(Af/Vf,Cf/Vf,Ef/Vf), name="ACE2f"),
    # Algebra to compute total variances and standard deviations (diagonal only)
        mxAlgebra( expression=Am+Cm+Em, name="Vm" ),
        mxAlgebra( expression=Af+Cf+Ef, name="Vf" ),
        mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
        mxAlgebra( expression=solve(sqrt(I*Vm)), name="iSDm"),
        mxAlgebra( expression=solve(sqrt(I*Vf)), name="iSDf"),
    # Constraint on variance of Bininal variables
        mxConstraint( expression=Vm==1, name="Var1"),
        mxConstraint( expression=Vf==1, name="Var2"),
    # Matrix & Algebra for expected means vector and expected thresholds
        mxMatrix( type="Zero", nrow=1, ncol=nv, name="Mean" ),
        mxAlgebra( expression= cbind(Mean,Mean), name="expMean" ),
        mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.8, label="thre", name="Threm" ),
        mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.8, label="thre", name="Thref" ),
        mxAlgebra( expression= cbind(Threm,Threm), dimnames=list('th1',selVars), name="expThreZm" ),
        mxAlgebra( expression= cbind(Thref,Thref), dimnames=list('th1',selVars), name="expThreZf" ),
        mxAlgebra( expression= cbind(Thref,Threm), dimnames=list('th1',selVars), name="expThreZfm" ),
    # Algebra for expected variance/covariance matrix in MZ
        mxAlgebra( expression= rbind  ( cbind(Am+Cm+Em , Am+Cm),
                                        cbind(Am+Cm   , Am+Cm+Em)), name="expCovMZm" ),
        mxAlgebra( expression= rbind  ( cbind(Af+Cf+Ef , Af+Cf),
                                        cbind(Af+Cf   , Af+Cf+Ef)), name="expCovMZf" ),
    # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
        mxAlgebra( expression= rbind  ( cbind(Am+Cm+Em     , 0.5%x%Am+Cm),
                                        cbind(0.5%x%Am+Cm , Am+Cm+Em)),  name="expCovDZm" ),
        mxAlgebra( expression= rbind  ( cbind(Af+Cf+Ef     , 0.5%x%Af+Cf),
                                        cbind(0.5%x%Af+Cf , Af+Cf+Ef)),  name="expCovDZf" ),
        mxAlgebra( expression= rbind  ( cbind(Af+Cf+Ef     , 0.5%*%rg%x%(af%*%t(am))+cf%*%t(cm)),
                                        cbind(0.5%*%rg%x%(am%*%t(af))+cm%*%t(cf) , Am+Cm+Em)),  name="expCovDZfm" ) 
    ),
    mxModel("MZm",
        mxData( observed=mzmDataBin, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovMZm", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreZm" )
    ),
    mxModel("DZm", 
        mxData( observed=dzmDataBin, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovDZm", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreZm" )
    ),
    mxModel("MZf",
        mxData( observed=mzfDataBin, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovMZf", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreZf" )
    ),
    mxModel("DZf", 
        mxData( observed=dzfDataBin, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovDZf", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreZf" )
    ),
    mxModel("DZfm", 
        mxData( observed=dzoDataBin, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovDZfm", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreZfm" )
    ),
    mxAlgebra( expression=MZm.objective + DZm.objective + MZf.objective + DZf.objective + DZfm.objective, 
        name="minus2sumloglikelihood" ),
    mxAlgebraObjective("minus2sumloglikelihood")
)

univHetBinACEFit <- mxRun(univHetBinACEModel)
univHetBinACESumm <- summary(univHetBinACEFit)
univHetBinACESumm

# Generate Heterogeneity ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(univHetBinACEFit)
expectedMeansCovariances(univHetBinACEFit)
tableFitStatistics(univHetBinACEFit)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(ACE.am,ACE.cm,ACE.em,ACE.af,ACE.cf,ACE.ef,ACE.rg), univHetBinACEFit),4)
varComponentsACE <- round(mxEval(cbind(ACE.Am/ACE.Vm,ACE.Cm/ACE.Vm,ACE.Em/ACE.Vm,ACE.Af/ACE.Vf,ACE.Cf/ACE.Vf,ACE.Ef/ACE.Vf), univHetBinACEFit),4)
    rownames(pathEstimatesACE) <- 'pathEstimates'
    colnames(pathEstimatesACE) <- c('am','cm','em','af','cf','ef','rg')
    rownames(varComponentsACE) <- 'varComponents'
    colnames(varComponentsACE) <- c('am^2','cm^2','em^2','af^2','cf^2','ef^2')
pathEstimatesACE
varComponentsACE

###############################################################################################################################################
# Homogeneity ACE Model 
# -----------------------------------------------------------------------
HomACEModel <- mxModel(univHetBinACEFit, name="HomACE",
    mxModel(univHetBinACEFit$ACE,
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="a11", name="am" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="c11", name="cm" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="e11", name="em" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="a11", name="af" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="c11", name="cf" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="e11", name="ef" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=FALSE, values=.5, label="rg11", name="rg" ) 
    )
)

HomACEModel <- mxOption(HomACEModel, "Function precision", 1e-9)
HomACEModelFit <- mxRun(HomACEModel)
HomACEModelSumm <- summary(HomACEModelFit)
HomACEModelSumm


# Generate homogeneity ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(HomACEModelFit)
expectedMeansCovariances(HomACEModelFit)
tableFitStatistics(HomACEModelFit)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(ACE.am,ACE.cm,ACE.em,ACE.af,ACE.cf,ACE.ef,ACE.rg), HomACEModelFit),4)
varComponentsACE <- round(mxEval(cbind(ACE.Am/ACE.Vm,ACE.Cm/ACE.Vm,ACE.Em/ACE.Vm,ACE.Af/ACE.Vf,ACE.Cf/ACE.Vf,ACE.Ef/ACE.Vf), HomACEModelFit),4)
    rownames(pathEstimatesACE) <- 'pathEstimates'
    colnames(pathEstimatesACE) <- c('am','cm','em','af','cf','ef','rg')
    rownames(varComponentsACE) <- 'varComponents'
    colnames(varComponentsACE) <- c('am^2','cm^2','em^2','af^2','cf^2','ef^2')
pathEstimatesACE
varComponentsACE


###############################################################################
#AE model

HomAEModel <- mxModel(univHetBinACEFit, name="HomAE",
    mxModel(univHetBinACEFit$ACE,
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="a11", name="am" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=FALSE, values=0, label="c11", name="cm" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="e11", name="em" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="a11", name="af" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=FALSE, values=0, label="c11", name="cf" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, label="e11", name="ef" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=FALSE, values=.5, label="rg11", name="rg" ) ))

HomAEModel <- mxOption(HomAEModel, "Function precision", 1e-9)
HomAEModelFit <- mxRun(HomAEModel)
HomAEModelSumm <- summary(HomAEModelFit)
HomAEModelSumm


# Generate Heterogeneity AE Output
# -----------------------------------------------------------------------
parameterSpecifications(HomAEModel)
expectedMeansCovariances(HomAEModel)
tableFitStatistics(HomAEModel)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(ACE.am,ACE.cm,ACE.em,ACE.af,ACE.cf,ACE.ef,ACE.rg), HomAEModelFit),4)
varComponentsACE <- round(mxEval(cbind(ACE.Am/ACE.Vm,ACE.Cm/ACE.Vm,ACE.Em/ACE.Vm,ACE.Af/ACE.Vf,ACE.Cf/ACE.Vf,ACE.Ef/ACE.Vf), HomAEModelFit),4)
    rownames(pathEstimatesACE) <- 'pathEstimates'
    colnames(pathEstimatesACE) <- c('am','cm','em','af','cf','ef','rg')
    rownames(varComponentsACE) <- 'varComponents'
    colnames(varComponentsACE) <- c('am^2','cm^2','em^2','af^2','cf^2','ef^2')
pathEstimatesACE
varComponentsACE


mxCompare(HomACEModelFit,HomAEModelFit)
mxCompare(univHetBinACEFit,HomAEModelFit)
mxCompare(univHetBinACEFit,HomACEModelFit)