ACE and ADE model returns the same AIC
Purple Fish
Joined: 05/23/2024
Hello everyone:
I'm currently working on a project with twin data, and I encountered some weird results when fitting ACE and ADE models with OpenMx.
I fit both bivariate ACE, AE, and ADE models based on the script from Dr.Hermine Maes' website. However, our model returns unusual results. Specifically, the AIC values for ACE and ADE are the same. Also, the parameter estimates for the ADE model do not match the phenotypic correlations. Has anyone encountered similar problem before? Any kind of help is appreciated!!
Below please find the code I used for ACE and ADE (I write them into for loops so it's easier to output results):
# ACE:
for (i in 1:12){
vars <- Vars_Twin[i,] # list of variables names
nv <- 2 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
mzData <- subset(Twin_data, ZYGOSITY.x == 1, selVars)
dzData <- subset(Twin_data, ZYGOSITY.x == 2, selVars)
svMe <- c(8,5) # start value for means
svPa <- .5 # start value for path coefficient
svPe <- .6 # start value for path coefficient for e
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe,
labels= labVars("mean",vars) , name="meanG" )
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv),
label=labLower("VA",nv), name="VA" )
covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv),
label=labLower("VC",nv), name="VC" )
covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv),
label=labLower("VE",nv), name="VE" )
covP <- mxAlgebra( expression= VA+VC+VE, name="V" )
covMZ <- mxAlgebra( expression= VA+VC, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%VA+ VC, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
pars <- list( meanG, covA, covC, covE, covP )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I" )
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD" )
corA <- mxAlgebra( expression=solve(sqrt(I*VA))%&%VA, name ="rA" ) #cov2cor()
corC <- mxAlgebra( expression=solve(sqrt(I*VC))%&%VC, name ="rC" )
corE <- mxAlgebra( expression=solve(sqrt(I*VE))%&%VE, name ="rE" )
rowUS <- rep('US',nv)
colUS <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv)
estUS <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="US", dimnames=list(rowUS,colUS) )
odd <- seq(1+3*nv,2*3*nv,nv)
even <- seq(2+3*nv,2*3*nv,nv)
ciACE <- mxCI( c("US[1,odd]","US[2,odd]","US[2,even]") )
calc <- list( matI, invSD, corA, corC, corE, estUS, ciACE )
modelACE <- mxModel( "twoACEvc", pars, modelMZ, modelDZ, multi, calc )
fitACE <- mxTryHard( modelACE, intervals=T )
sumACE <- summary( fitACE, verbose = T)
A11_ACE[i] <- mxEval(VA11, fitACE)
A21_ACE[i] <- mxEval(VA21, fitACE)
A22_ACE[i] <- mxEval(VA22, fitACE)
C11_ACE[i] <- mxEval(VC11, fitACE)
C21_ACE[i] <- mxEval(VC21, fitACE)
C22_ACE[i] <- mxEval(VC22, fitACE)
E11_ACE[i] <- mxEval(VE11, fitACE)
E21_ACE[i] <- mxEval(VE21, fitACE)
E22_ACE[i] <- mxEval(VE22, fitACE)
a1_sqr_ACE[i] <- A11_ACE[i]^2 / (A11_ACE[i]^2 + C11_ACE[i]^2 + E11_ACE[i]^2)
a2_sqr_ACE[i] <- (A21_ACE[i]^2 + A22_ACE[i]^2) / (A21_ACE[i]^2 + A22_ACE[i]^2 + C21_ACE[i]^2 + C22_ACE[i]^2 + E21_ACE[i]^2 + E22_ACE[i]^2)
c1_sqr_ACE[i] <- C11_ACE[i]^2 / (A11_ACE[i]^2 + C11_ACE[i]^2 + E11_ACE[i]^2)
c2_sqr_ACE[i] <- (C21_ACE[i]^2 + C22_ACE[i]^2) / (A21_ACE[i]^2 + A22_ACE[i]^2 + C21_ACE[i]^2 + C22_ACE[i]^2 + E21_ACE[i]^2 + E22_ACE[i]^2)
e1_sqr_ACE[i] <- E11_ACE[i]^2 / (A11_ACE[i]^2 + C11_ACE[i]^2 + E11_ACE[i]^2)
e2_sqr_ACE[i] <- (E21_ACE[i]^2 + E22_ACE[i]^2) / (A21_ACE[i]^2 + A22_ACE[i]^2 + C21_ACE[i]^2 + C22_ACE[i]^2 + E21_ACE[i]^2 + E22_ACE[i]^2)
r_a_ACE[i] <- A21_ACE[i] / sqrt(A21_ACE[i]^2 + A22_ACE[i]^2)
r_c_ACE[i] <- C21_ACE[i] / sqrt(C21_ACE[i]^2 + C22_ACE[i]^2)
r_e_ACE[i] <- E21_ACE[i] / sqrt(E21_ACE[i]^2 + E22_ACE[i]^2)
sumACE_par[[i]] <- sumACE$parameters
sumACE_AIC[i] <- sumACE$informationCriteria[1,3]
sumACE_CI[[i]] <- sumACE$CI
}
# ADE:
for (i in 1:12){
vars <- Vars_Twin[i,] # list of variables names
nv <- 2 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
mzData <- subset(Twin_data, ZYGOSITY.x == 1, selVars)
dzData <- subset(Twin_data, ZYGOSITY.x == 2, selVars)
svMe <- c(8,5) # start value for means
svPa <- .5 # start value for path coefficient
svPe <- .6 # start value for path coefficient for e
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" )
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VA",nv), name="VA" )
covD <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VD",nv), name="VD" )
covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VE",nv), name="VE" )
covP <- mxAlgebra( expression= VA+VD+VE, name="V" )
covMZ <- mxAlgebra( expression= VA+VD, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%VA+ 0.25%x%VD, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind(cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind(cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
pars <- list( meanG, covA, covD, covE, covP )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I" )
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD" )
corA <- mxAlgebra( expression=solve(sqrt(I*VA))%&%VA, name ="rA" ) #cov2cor()
corD <- mxAlgebra( expression=solve(sqrt(I*VD))%&%VD, name ="rD" )
corE <- mxAlgebra( expression=solve(sqrt(I*VE))%&%VE, name ="rE" )
rowUS <- rep('US',nv)
colUS <- rep(c('VA','VD','VE','SA','SD','SE'),each=nv)
estUS <- mxAlgebra( expression=cbind(VA,VD,VE,VA/V,VD/V,VE/V), name="US", dimnames=list(rowUS,colUS) )
odd <- seq(1+3*nv,2*3*nv,nv)
even <- seq(2+3*nv,2*3*nv,nv)
ciADE <- mxCI( c("US[1,odd]","US[2,odd]","US[2,even]") )
calc <- list( matI, invSD, corA, corD, corE, estUS, ciADE )
modelADE <- mxModel( "twoADEvc", pars, modelMZ, modelDZ, multi, calc )
fitADE <- mxTryHard( modelADE, intervals=T )
sumADE <- summary( fitADE, verbose = T)
A11_ADE[i] <- mxEval(VA11, fitADE)
A21_ADE[i] <- mxEval(VA21, fitADE)
A22_ADE[i] <- mxEval(VA22, fitADE)
D11_ADE[i] <- mxEval(VD11, fitADE)
D21_ADE[i] <- mxEval(VD21, fitADE)
D22_ADE[i] <- mxEval(VD22, fitADE)
E11_ADE[i] <- mxEval(VE11, fitADE)
E21_ADE[i] <- mxEval(VE21, fitADE)
E22_ADE[i] <- mxEval(VE22, fitADE)
a1_sqr_ADE[i] <- A11_ADE[i]^2 / (A11_ADE[i]^2 + D11_ADE[i]^2 + E11_ADE[i]^2)
a2_sqr_ADE[i] <- (A21_ADE[i]^2 + A22_ADE[i]^2) / (A21_ADE[i]^2 + A22_ADE[i]^2 + D21_ADE[i]^2 + D22_ADE[i]^2 + E21_ADE[i]^2 + E22_ADE[i]^2)
d1_sqr_ADE[i] <- D11_ADE[i]^2 / (A11_ADE[i]^2 + D11_ADE[i]^2 + E11_ADE[i]^2)
d2_sqr_ADE[i] <- (D21_ADE[i]^2 + D22_ADE[i]^2) / (A21_ADE[i]^2 + A22_ADE[i]^2 + D21_ADE[i]^2 + D22_ADE[i]^2 + E21_ADE[i]^2 + E22_ADE[i]^2)
e1_sqr_ADE[i] <- E11_ADE[i]^2 / (A11_ADE[i]^2 + D11_ADE[i]^2 + E11_ADE[i]^2)
e2_sqr_ADE[i] <- (E21_ADE[i]^2 + E22_ADE[i]^2) / (A21_ADE[i]^2 + A22_ADE[i]^2 + D21_ADE[i]^2 + D22_ADE[i]^2 + E21_ADE[i]^2 + E22_ADE[i]^2)
r_a_ADE[i] <- A21_ADE[i] / sqrt(A21_ADE[i]^2 + A22_ADE[i]^2)
r_d_ADE[i] <- D21_ADE[i] / sqrt(D21_ADE[i]^2 + D22_ADE[i]^2)
r_e_ADE[i] <- E21_ADE[i] / sqrt(E21_ADE[i]^2 + E22_ADE[i]^2)
sumADE_par[[i]] <- sumADE$parameters
sumADE_AIC[i] <- sumADE$informationCriteria[1,3]
sumADE_CI[[i]] <- sumADE$CI
}
parameterization
Without looking at your script in detail, I am pretty sure the key part is this:
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VA",nv), name="VA" )
covD <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VD",nv), name="VD" )
covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VE",nv), name="VE" )
The script uses what is sometimes called the "direct-symmetric" parameterization, under which the *A*, *C*, and *E* covariance matrices are directly specified as symmetric matrices. Under this parameterization--and in contrast to the "Cholesky" parameterization--there is nothing to ensure that they be positive-definite (i.e., they are not guaranteed to be true covariance matrices). So, you may find that 'VA', 'VC', and 'VE' have negative eigenvalues, or even that some variance components are negative! Because such point estimates are allowed by this parameterization, the ACE and the ADE models can, and often do, have the same AIC.
See here for the statistical rationale for the direct-symmetric parameterization.
If you're getting nonsensical results from the ADE model, then that's a sign that it's a poor fit to the data.
Log in or register to post comments
In reply to parameterization by AdminRobK
Thanks a lot for your clear
Thanks a lot for your clear explanation! I’m still pretty new to OpenMx, so I'm trying to figure out if I should experiment with different "type" options, like "Full", or tweak the code for a "Cholesky" parameterization to make sure the matrices are positive-definite. Could you give me some advice on this?
Also, about the ADE model, my professor tried it with some older software, and the results seemed to make more sense with the observed correlations. It would be great to get this sorted out in OpenMx, since the other software is quite outdated. Any suggestions you have would be really helpful.
Log in or register to post comments
In reply to Thanks a lot for your clear by Purple Fish
Hermine Maes has scripts here
Hermine Maes has scripts here that use the Cholesky parameterization; they are marked as "estimating path coefficients".
I can't speak to the analysis in the other software (classic Mx?) without seeing any of its syntax and results.
Log in or register to post comments
In reply to Hermine Maes has scripts here by AdminRobK
Thank you for your help!
I see. I think I'll be able to figure the rest of the analysis out at this point, and maybe compare OpenMx results with the other software to see if they look similar. Thank you so much for being so helpful. Appreciate that a lot!
Log in or register to post comments