Dear experts,
I am running a simple ACE model for twin analysis with age and gender as confounds (115 DZ pairs and 60 MZ pairs). I use the model from International Statistical Genetics Workshop:
https://ibg.colorado.edu/cdrom2022/day2/00_ACEvc_contin.R
When I run the model, I get the following results:
free parameters:
name matrix row col Estimate Std.Error A lbound ubound
1 interC intercept 1 1 1.5044971 0.64957699
2 betaS bS 1 1 -0.5795265 0.12804871
3 betaA bA 1 1 -0.0236737 0.01893570
4 VA11 VA 1 1 0.9179206 0.34811542
5 VC11 VC 1 1 -0.4279909 0.32400691
6 VE11 VE 1 1 0.4132866 0.05129183
confidence intervals:
lbound estimate ubound note
ACEvc.VarC[1,4] 0.3595914 1.0162799 1.7306621
ACEvc.VarC[1,5] -1.1498335 -0.4738520 0.1515639
ACEvc.VarC[1,6] 0.3540964 0.4575721 0.5850078
CI details:
parameter side value fit diagnostic statusCode method interC betaS betaA VA11 VC11 VE11 1 ACEvc.VarC[1,4] lower 0.3595914 921.3400 success OK neale-miller-1997 1.397043 -0.5711129 -0.02000247 0.3277627 0.1253068 0.4584169 2 ACEvc.VarC[1,4] upper 1.7306621 921.3660 success OK neale-miller-1997 1.643095 -0.5960471 -0.02823600 1.6357127 -1.0803400 0.3897642 3 ACEvc.VarC[1,5] lower -1.1498335 921.3463 success OK neale-miller-1997 1.655002 -0.5904675 -0.02897582 1.6119600 -1.0772610 0.4021853 4 ACEvc.VarC[1,5] upper 0.1515639 921.3503 success OK neale-miller-1997 1.382750 -0.5690776 -0.01958674 0.3453931 0.1397735 0.4370419 5 ACEvc.VarC[1,6] lower 0.3540964 921.3554 success OK neale-miller-1997 1.474451 -0.5694362 -0.02322113 1.1770877 -0.5488576 0.3444075 6 ACEvc.VarC[1,6] upper 0.5850078 921.3530 success OK neale-miller-1997 1.522639 -0.5798345 -0.02417389 0.6512583 -0.2910815 0.5077355
Critically, I also get the following error:
Running CE with 5 parameters
Error: The job for model 'CE' exited abnormally with the error message: fit is not finite (The continuous part of the model implied covariance (loc2) is not positive definite in data 'DZ.data' row 59. Detail:
covariance = matrix(c( # 2x2
-0.014704266851264, -0.427990888364092
, -0.427990888364092, -0.014704266851264), byrow=TRUE, nrow=2, ncol=2)
)
Following suggestion in this post: https://openmx.ssri.psu.edu/node/4593 I set the lbound=1e-4 to covC.
As a result, the error disappears, but the results become very different to those I got before. In addition, I am not sure how to treat the confidence interval of VA, given that it is constrained.
free parameters:
name matrix row col Estimate Std.Error A lbound ubound
1 interC intercept 1 1 1.41722921 0.67171698
2 betaS bS 1 1 -0.57166848 0.13325128
3 betaA bA 1 1 -0.02073634 0.01960500
4 VA11 VA 1 1 0.48480850 0.25577544
5 VC11 VC 1 1 0.00010000 0.24940786 ! 1e-04!
6 VE11 VE 1 1 0.42938467 0.05552459
confidence intervals:
lbound estimate ubound note
ACEvc.VarC[1,4] 0.227153 0.5302549731 NA !!!
ACEvc.VarC[1,5] NA 0.0001093741 0.2708948 !!!
ACEvc.VarC[1,6] NA 0.4696356528 NA !!!
CI details:
parameter side value fit diagnostic statusCode method interC betaS betaA VA11 VC11 VE11 1 ACEvc.VarC[1,4] lower 2.271530e-01 923.3903 success OK neale-miller-1997 1.383110 -0.5685890 -0.01955053 0.2076383 0.2297972 0.4766547 2 ACEvc.VarC[1,4] upper 6.390220e-01 923.3817 active box constraint OK neale-miller-1997 1.406398 -0.5723281 -0.02036513 0.6334840 0.0001000 0.3577496 3 ACEvc.VarC[1,5] lower 9.411526e-05 923.4075 active box constraint OK neale-miller-1997 1.417229 -0.5716684 -0.02073514 0.5589758 0.0001000 0.5034512 4 ACEvc.VarC[1,5] upper 2.708948e-01 923.3912 success OK neale-miller-1997 1.359038 -0.5664194 -0.01889902 0.2312665 0.2535373 0.4511215 5 ACEvc.VarC[1,6] lower 3.606429e-01 923.4018 active box constraint OK neale-miller-1997 1.412693 -0.5717267 -0.02061276 0.6350086 0.0001000 0.3582465 6 ACEvc.VarC[1,6] upper 6.013603e-01 923.3829 active box constraint OK neale-miller-1997 1.444580 -0.5744477 -0.02149105 0.3475080 0.0001000 0.5243773
How should I interpret these results? Which of the models is correct? According to Falconer’s formula, h-square is about 1, which is close to the first result. But in the first results, not only that I get an error, but there is also negative VC. I would very much appreciate your help .
Best wishes,
Vadim
Printer-friendly version
Hi
When the model is rerun to fit CE, I suspect that it is doing so from the solution of the ACE model, and the C component happens to be negative, so the CE model has a likelihood that cannot be computed. To avoid this, modify the ACE model before it was run, instead of the result of running the ACE model.
Thank you very much for your answer, but can you please elaborate a bit more. I will very appreciate it.
The fitACE is indeed passed as the parameter to estimation of CE.
fitCE <- omxSetParameters(fitACE, labels = "VA11", values = 0, free = F, name = "CE")
fitCE <- mxRun( fitCE, intervals=F )
So, what should I modify?
Try this:
Thank you for your help.
When I run the original code (without lbound constraint) after I replaced my two lines with your code, the error disappeared. But the variance result is the same as before (with the negative variance VC):
name matrix row col Estimate Std.Error A
1 interC intercept 1 1 1.53784349 0.65261107
2 betaS bS 1 1 -0.58236697 0.12811632
3 betaA bA 1 1 -0.02468047 0.01903097
4 betaM bM 1 1 0.02468528 0.04961712
5 VA11 VA 1 1 0.92411484 0.35230753
6 VC11 VC 1 1 -0.43334699 0.32837219
7 VE11 VE 1 1 0.41226047 0.05112496