First time poster here. I'm learning generally about SEM and openMX at the same time and I want to run an ACE bivariate Cholesky decomposition with a sample of twins. I'm using the script below to test assumptions of the model and then run the Cholesky. (This is very similar to existing scripts that are posted, I just added some elements to get confidence intervals for the genetic/shared environmental correlation, and also the phenotypic correlation.)
My question is: the bivariate saturated model and the existing submodels do not appear to constrain the cross-twin cross-trait correlations so that they are symmetric across the phenotypes (e.g., that CTCT of Twin 1 Variable 1 and Twin 2 Variable 2 = CTCT of Twin 1 Variable 2 and Twin 2 Variable 1). Since the CTCT correlation is typically reported in manuscripts for MZ and DZ twins, are authors typically adding additional constraints beyond what's listed below to get this statistic?
# 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
In a saturated model, it is necessary to constrain correlations (usually via an mxConstraint() call, e.g.,
In a model with variance components (the second one fitted above, bivACEModel) it isn't necessary because the constraint is implicit.
A good question from a novice user, that!