ACE: negative variance component from H.Maes’ published script and slides
lemire
Joined: 12/13/2019
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
boundary conditions
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 argument `verbose=TRUE` to `summary()`, 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".
Log in or register to post comments