You are here

ACE and ADE model returns the same AIC

5 posts / 0 new
Last post
Purple Fish's picture
Offline
Joined: 05/23/2024 - 12:30
ACE and ADE model returns the same AIC

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
 
}
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
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.

Purple Fish's picture
Offline
Joined: 05/23/2024 - 12:30
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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
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.

Purple Fish's picture
Offline
Joined: 05/23/2024 - 12:30
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!