Attachment | Size |
---|---|
Contrast.mx [6] | 7.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