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 }
Printer-friendly version
Without looking at your script in detail, I am pretty sure the key part is this:
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.
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.
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.
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!