You are here

Bivariate Cholesky Assumptions and CTCT Correlations

2 posts / 0 new
Last post
JuliaJoplin's picture
Joined: 10/12/2014 - 16:44
Bivariate Cholesky Assumptions and CTCT Correlations

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",
        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)
        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" ),
    mxCI(c('MZ.expCorMZ', 'DZ.expCorDZ'))
bivTwinSatFit <- mxRun(bivTwinSatModel)
bivTwinSatSumm <- summary(bivTwinSatFit)
# Generate Saturated Output
# -----------------------------------------------------------------------
# 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)
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)
tableFitStatistics(bivTwinSatFit, list(equateMeansVarsTwinFit, equateMeansVarsTwinZygFit))
bivACEModel <- mxModel("bivACE",
    # 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" ) 
        mxData( observed=mzData, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars )
        mxData( observed=dzData, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars ) 
    mxAlgebra( expression=MZ.objective + DZ.objective, name="minus2sumloglikelihood" ),
    mxCI(c('ACE.sta', '', '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)
neale's picture
Joined: 07/31/2009 - 15:14
Yes and no

In a saturated model, it is necessary to constrain correlations (usually via an mxConstraint() call, e.g.,

mxConstraint(expCorMZ[4,1]==expCorMZ[3,2], name="EqCor")

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!