Bivariate Cholesky Assumptions and CTCT Correlations

# Fit Bivariate Saturated Model
# -----------------------------------------------------------------------
bivTwinSatModel <- mxModel("bivTwinSat",
mxModel("MZ",
mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="CholMZ" ),
mxAlgebra( expression=CholMZ %*% t(CholMZ), name="expCovMZ" ),
mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="I"),
mxAlgebra( solve(sqrt(I*expCovMZ)), name="iSDMZ"),
mxAlgebra( iSDMZ%*%expCovMZ%*%iSDMZ, name="expCorMZ"),
mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=20, name="expMeanMZ" ),
mxData( observed=mzData, type="raw" ),
mxFIMLObjective( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars)
),
mxModel("DZ",
mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="CholDZ" ),
mxAlgebra( expression=CholDZ %*% t(CholDZ), name="expCovDZ" ),
mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="I"),
mxAlgebra( solve(sqrt(I*expCovDZ)), name="iSDDZ"),
mxAlgebra( iSDDZ%*%expCovDZ%*%iSDDZ, name="expCorDZ"),
mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=20, name="expMeanDZ" ),
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars)
),
mxAlgebra( MZ.objective + DZ.objective, name="minus2sumloglikelihood" ),
mxAlgebraObjective("minus2sumloglikelihood"),
mxCI(c('MZ.expCorMZ', 'DZ.expCorDZ'))
)

bivTwinSatFit <- mxRun(bivTwinSatModel)
bivTwinSatSumm <- summary(bivTwinSatFit)
bivTwinSatSumm

# Generate Saturated Output
# -----------------------------------------------------------------------
parameterSpecifications(bivTwinSatFit)
expectedMeansCovariances(bivTwinSatFit)
tableFitStatistics(bivTwinSatFit)

# Fit Bivariate Model with Equal Means & Variances across Twin Order and Zygosity
# -----------------------------------------------------------------------
# Constrain expected means and variances to be equal across twin order
equateMeansVarsTwinModel <- mxModel(bivTwinSatFit, name="equateMeansVarsTwin",
mxAlgebra( expression=t(diag2vec(MZ.expCovMZ)), name="expVarMZ"),
mxAlgebra( expression=t(diag2vec(DZ.expCovDZ)), name="expVarDZ"),
mxConstraint( expression=expVarMZ[1,1:nv]==expVarMZ[1,(nv+1):ntv], name="VarMZt1eqt2"),
mxConstraint( expression=expVarDZ[1,1:nv]==expVarDZ[1,(nv+1):ntv], name="VarDZt1eqt2"),
mxConstraint( expression=MZ.expMeanMZ[1,1:nv]==MZ.expMeanMZ[1,(nv+1):ntv], name="MeanMZt1eqt2"),
mxConstraint( expression=DZ.expMeanDZ[1,1:nv]==DZ.expMeanDZ[1,(nv+1):ntv], name="MeanDZt1eqt2")
)
equateMeansVarsTwinFit <- mxRun(equateMeansVarsTwinModel)
equateMeansVarsTwinSumm <- summary(equateMeansVarsTwinFit)
equateMeansVarsTwinSumm
parameterSpecifications(equateMeansVarsTwinFit)
expectedMeansCovariances(equateMeansVarsTwinFit)
tableFitStatistics(bivTwinSatFit, equateMeansVarsTwinFit)

# Constrain expected means and variances to be equal across twin order and zygosity
equateMeansVarsTwinZygModel <- mxModel(equateMeansVarsTwinFit, name="equateMeansVarsTwinZyg",
mxConstraint( expression=expVarMZ[1,1:nv]==expVarDZ[1,1:nv], name="VarMZeqDZ"),
mxConstraint( expression=MZ.expMeanMZ[1,1:nv]==DZ.expMeanDZ[1,1:nv], name="MeanMZeqDZ")
)
equateMeansVarsTwinZygFit <- mxRun(equateMeansVarsTwinZygModel,intervals=T)
equateMeansVarsTwinZygSumm <- summary(equateMeansVarsTwinZygFit)
equateMeansVarsTwinZygSumm
parameterSpecifications(equateMeansVarsTwinZygFit)
tableFitStatistics(bivTwinSatFit, list(equateMeansVarsTwinFit, equateMeansVarsTwinZygFit))

bivACEModel <- mxModel("bivACE",
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=.8, lbound=0, labels=c("a11","a21","a22"), name="a" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=0, labels=c("c11","c21","c22"),name="c" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=0, labels=c("e11","e21","e22"),name="e" ),
# 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"),
mxAlgebra(iSD %*% a, name="sta"),
mxAlgebra(iSD %*% c, name="stc"),
mxAlgebra(iSD %*% e, name="ste"),
mxAlgebra(A/V, name="StandardizedA"),
mxAlgebra(C/V, name="StandardizedC"),
mxAlgebra(E/V, name="StandardizedE"),
mxAlgebra(solve(sqrt(I*A)) %&% A, name="CorA"),
mxAlgebra(solve(sqrt(I*C)) %&% C, name="CorC"),
mxAlgebra(solve(sqrt(I*E)) %&% E, name="CorE"),
mxAlgebra(solve(sqrt(I*V)) %&% V, name="CorP"),
## 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= 20, 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="minus2sumloglikelihood" ),
mxAlgebraObjective("minus2sumloglikelihood"),
mxCI(c('ACE.sta', 'ACE.stc', 'ACE.ste')),
mxCI(c('ACE.StandardizedA', 'ACE.StandardizedC', 'ACE.StandardizedE')),
mxCI(c('ACE.CorA', 'ACE.CorC', 'ACE.CorE', 'ACE.CorP'))
)
bivACEFit <- mxRun(bivACEModel)
bivACESumm <- summary(bivACEFit)
bivACESumm