You are here

Longitudinal Cholesky decomposition with contrast effect

4 posts / 0 new
Last post
orangeU's picture
Offline
Joined: 05/23/2021 - 04:01
Longitudinal Cholesky decomposition with contrast effect
AttachmentSize
Binary Data Contrast.mx7.86 KB

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!

 # 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
AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Try this?

Hi Reut!

# 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

orangeU's picture
Offline
Joined: 05/23/2021 - 04:01
replying properly

Hi Mike,
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

orangeU's picture
Offline
Joined: 05/23/2021 - 04:01
Thank you!

Thank you so so much Mike!
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:

Vars <- c("EXT3", "EXT5")
nvar <- length(Vars)               # number of variables
tnvar <- nvar*2                 # number of variables*max family size
nlower <-nvar*(nvar+1)/2           # number of free elements in a lower matrix nvar*nvar
 
selVars <- paste(Vars, c(rep(1,nvar), rep(2,nvar)), sep="_")
 
## using both DZ-S and DZ-O.
mzData <- as.matrix(subset(twin_data, zygDich==1, selVars))
dzData <- as.matrix(subset(twin_data, zygDich==0, selVars))
 
describe(mzData, skew=F)
describe(dzData, skew=F)
dim(mzData)
dim(dzData)
 
# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
cor(mzData,use="complete")
cor(dzData,use="complete")
 
nv <- length(Vars)                 # number of variables
#tnvar <- nvar*2               # number of variables*max family size
#nlower <-nvar*(nvar+1)/2          # number of free elements in a lower matrix nvar*nvar
 
 
# Set Starting Values
meVals    <- c(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 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" )
 
# 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, bMat, iMat, sibInt)
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

PS. Sorry for not providing the full script before! Wasn't sure what is needed.