Hi,
could you help provide some insights as to why running the following script, which uses data(twinData), gives negative estimates for Vc in an ACE model?
https://ibg.colorado.edu/cdrom2020/maes/twinModeling/oneACEvc.R
These are the most relevant lines related to the model:
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" ) covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" ) covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VE11", name="VE" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins 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" )
Running the script provides the following solution:
fitACE$output$confidenceIntervals lbound estimate ubound oneACEvc.US[1,1] 0.60785108 0.75506776 0.9254683654 oneACEvc.US[1,2] -0.29756250 -0.14471244 -0.0061102627 oneACEvc.US[1,3] 0.15055312 0.16935009 0.1913908419
These estimates are further reported on page 71 of the accompanying slides found at
https://ibg.colorado.edu/cdrom2020/maes/twinModeling/TwinModeling2020c.pdf
I don’t really understand why these results would make sense. Of course, just looking at the slides I don’t have the full narrative, which is why I am asking here.
If VA/VC/VE are known to be variances, wouldn't including lbound=0 in the mxMatrix functions be more appropriate? Was that just an oversight? The slides further discuss negative C in the next slide, so it doesn't appear it was.
Furthermore, if setting lbound=0 is appropriate, I noticed that the confidence intervals failed to be calculated:
lbound estimate ubound oneACEvc.US[1,1] 0.54832905 6.1730227e-01 NA oneACEvc.US[1,2] NA 1.3552200e-12 NA oneACEvc.US[1,3] NA 1.7304629e-01 NA
can this be fixed?
Thanks for your help,
Mathieu
See this other post concerning the negative variance components. The matter of negative variance components has been discussed a number of times on these forums, so the OpenMx website's search feature may be useful for finding previous such discussion.
The Verhulst et al. paper linked in that other post addresses why some of the confidence limits aren't shown when you add the lbounds: boundary conditions. If there is one or more active bounds near the confidence-limit solution, the limit may be wrong in that it makes the CI too narrow. Therefore, OpenMx does not print it in the default
summary()
output. You can see a table of "CI details" if you pass argumentverbose=TRUE
tosummary()
, and I expect that table will show the confidence limits missing from the default output to be flagged as either "active box constraints" or "alpha level not reached".