Longitudinal Cholesky decomposition with contrast effect
Posted on
orangeU
Joined: 05/23/2021
Attachment | Size |
---|---|
Contrast.mx | 7.86 KB |
Forums
Hi all,
I have negative DZ correlations across multiple ages and I read that there may be a way to account for this by including a contrast effect (representing possible effects of twins within a pair on one another) in the model. Could someone please share an openMx script that does this? I have an old Mx script (attached) that accounts for a contrast effect in a classic ACE model, but I don't know how to translate it to openMx or how to integrate it in a Cholesky analysis.
Below is the script I'm currently using, if anyone has any suggestions for ways to incorporate a contrast effect in this script I will greatly appreciate it.
Many thanks in advance!
I have negative DZ correlations across multiple ages and I read that there may be a way to account for this by including a contrast effect (representing possible effects of twins within a pair on one another) in the model. Could someone please share an openMx script that does this? I have an old Mx script (attached) that accounts for a contrast effect in a classic ACE model, but I don't know how to translate it to openMx or how to integrate it in a Cholesky analysis.
Below is the script I'm currently using, if anyone has any suggestions for ways to incorporate a contrast effect in this script I will greatly appreciate it.
Many thanks in advance!
# Set Starting Values
meVals <- c(0.5,0.5,0.5) # start value for means
meanLabs <- paste(Vars,"mean",sep="") # create labels for means
paVal <- .6 # start value for path coefficient
paLBo <- .0001 # start value for lower bounds
paValD <- vech(diag(paVal,nv,nv)) # start values for diagonal of covariance matrix
paLBoD <- diag(paLBo,nv,nv) # lower bounds for diagonal of covariance matrix
paLBoD[lower.tri(paLBoD)] <- -10 # lower bounds for below diagonal elements
paLBoD[upper.tri(paLBoD)] <- NA # lower bounds for above diagonal elements
# Create Labels for Lower Triangular Matrices
aLabs <- paste("a",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")
cLabs <- paste("c",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")
eLabs <- paste("e",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")
# Matrices declared to store a, c, and e Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=paValD, labels=aLabs, name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=paValD, labels=cLabs, name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=paValD, labels=eLabs, name="e" )
# Matrices generated to hold A, C, and E computed Variance Components
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
# Algebra to compute total variances and standard deviations (diagonal only)
covP <- mxAlgebra( expression=A+C+E, name="V" )
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")
# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE,
values=meVals, labels=meanLabs, name="Mean" )
meanT <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" )
covMZ <- mxAlgebra( expression=
rbind( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" )
covDZ <- mxAlgebra( expression=
rbind( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
# Standeraized parameters and CI (getting CIs for squared path estimates)
sta <- mxAlgebra((iSD %*% a)*(iSD %*% a), name="sta")
stc <- mxAlgebra((iSD %*% c)*(iSD %*% c), name="stc")
ste <- mxAlgebra((iSD %*% e)*(iSD %*% e), name="ste")
StandardizedA <- mxAlgebra(A/V, name="StandardizedA")
StandardizedC <- mxAlgebra(C/V, name="StandardizedC")
StandardizedE <- mxAlgebra(E/V, name="StandardizedE")
CI <- mxCI(c('sta', 'stc', 'ste'))
# Data objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean",
dimnames=selVars )
objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean",
dimnames=selVars )
funML <- mxFitFunctionML()
# Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP,
matI, invSD, meanG, meanT, sta, stc, ste,
StandardizedA, StandardizedC, StandardizedE)
modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, funML, name="DZ" )
fitML <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction") )
CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, fitML, CI)
# Run Cholesky Decomposition ACE model
CholAceFit <- mxRun(CholAceModel, intervals = TRUE)
CholAceSumm <- summary(CholAceFit, verbose=T)
CholAceSumm
Try this?
# Algebra to obtain solve(I-B) for sib interaction
bMat <- mxMatrix(type="Full", nrow=2*nv, ncol=2*nv, free=c(F,T,T,F), labels=c(NA,"b","b",NA), name="bMat")
iMat <- mxMatrix(type="Iden", nrow=2*nv, ncol=2*nv, name="iMat")
sibInt <- mxAlgebra(solve(iMat-bMat), name="sibInt")
# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE,
values=meVals, labels=meanLabs, name="Mean" )
meanT <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" )
covMZ <- mxAlgebra( expression= sibInt %&% (
rbind( cbind(A+C+E , A+C),
cbind(A+C , A+C+E))), name="expCovMZ" )
covDZ <- mxAlgebra( expression= sibInt %&% (
rbind( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E))), name="expCovDZ" )
And then add the new objects to the pars list:
# Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP,
matI, invSD, meanG, meanT, sta, stc, ste,
StandardizedA, StandardizedC, StandardizedE, bMat, iMat, sibInt)
I've not tried this out as you provided a partial script, but I hope it helps!
- Mike
Log in or register to post comments
In reply to Try this? by AdminNeale
replying properly
I accidentally replied to my own post instead of yours, so I just wanted to make sure you see it.
Many many thanks again for your help!
Reut
Log in or register to post comments
Thank you!
I included the additions you suggested and the model does run. However, I received the attached Mx script a few days after I posted my inquiry here (it can only handle 2 age groups, and I wasn't able to figure out how to change it to make it fit 4 ages, which is what I need), and the results for 2 age groups seem to differ between Mx and openMx. I'm probably missing something simple here, but I'm not sure what. (I attached the openMx output, the Mx script and the output from Mx.)
Would love to hear your thoughts on this and many thanks again for your time!!
This is the script with your suggestions:
PS. Sorry for not providing the full script before! Wasn't sure what is needed.
Log in or register to post comments