ACE and ADE model returns the same AIC

Posted on Thu, 05/23/2024 - 15:31

Purple Fish Joined: 05/23/2024

Forums

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

}

Replied on Fri, 05/24/2024 - 09:51

AdminRobK Joined: 01/24/2014

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

Replied on Fri, 05/24/2024 - 10:45

Purple Fish Joined: 05/23/2024

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

Replied on Wed, 05/29/2024 - 09:04

AdminRobK Joined: 01/24/2014

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

Replied on Wed, 05/29/2024 - 21:43

Purple Fish Joined: 05/23/2024

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

Log in or register to post comments

News

Recent Posts

ACE and ADE model returns the same AIC

parameterization

Thanks a lot for your clear

Hermine Maes has scripts here

Thank you for your help!

News

Recent Posts