Hi all,
I'm running a Qualitative Sex Differences ACE model and I'm using the mxRefModels to get the reference models, but for some variables I get a error:
>mxRefModels(QualAceFit_rg, run = T) Error in chol.default(sampcov) : the leading minor of order 1 is not positive definite > mxVersion() OpenMx version: 2.3.1 [GIT v2.3.1] R version: R version 3.2.1 (2015-06-18) Platform: x86_64-w64-mingw32 Default optimiser: NPSOL
If I define the Saturated and the Independence models by hand, then there doesn’t seem to be a problem, it runs and converges correctly. Could someone help, what might be the problem with mxRefModels? It's a great function, really speeds up things, I hope I can use it for these models also.
Printer-friendly version
Like the documentation for mxRefModels says, be careful using it for twin and family modeling. It may or may not give you the saturated model you want (depends on what you want), particularly if it has definition variables.
With that disclaimer aside, it looks to me like the function is having a hard time setting starting values for your model. The error (Error in chol.default(sampcov)) occurs when it's trying to find good starting values by computing the sample covariance matrix on the pairwise complete data (sampcov) and then perform the Cholesky decomposition on that. In your case, this pairwise complete covariance matrix is not positive definite. It can happen due to missing data.
On the development side, we may want to have a fallback/failsafe for this situation so you don't get this error.
On the user side, here's something you can try. First, try swapping out listwise deleted data in the model before running mxRefModels. Run mxRefModels with run=FALSE. Then overwrite the data with the full data set. Here's a code snippet.
Although Cholesky might give decent starting values for the saturated model in the case that the observed covariance matrix is positive definite, it seems to me that there are sensible things to do if it is not. I would just get the means right and use the sample-wide standard deviations to populate the diagonal of the Cholesky. Unless a variable has zero variance, or there's a terrible outlier, optimization should be able to get off the mark. This behavior should make the function more robust. It would be robust against collinearity also (though the estimated saturated model covariance might try to go non-positive definite, but that's a different error message.)
That's a fine idea! As a side note, the mxRefModels function always estimates the Cholesky of the covariance matrix. So it is not possible for it to become non-positive definite.
A Cholesky only guarantees non-negative definiteness. If you put a >= SMALL lower bound on the diagonal elements of the Cholesky, where SMALL is the smallest number recognized as >0 in machine precision, then its product with its transpose (expCov) will be positive definite. Otherwise, expCov can be singular (a zero in any Cholesky diagonal element will do).
Thak you Michael for your quick responce!
I'm not using any definition variables. With other parameters, where I don't get this error, my hand calculated Saturated and Independence model values are identical to the ones mxRefModels provides.
I checked my data, but I don't have any missing data. I tried setting "run = FALSE", but that doesn't change anything, the error message is the same, while using the hand built models, both the Indepenence and the Saturated converge with code 0 on the same data set. Could it be that the correlation values between the genetic components (rg) is what is causing the problem using mxRefModels? If so, why isn't it a problem only with some datasets?
Thanks for your help!
Hmm ... no missing data. That's odd. I'm wondering now if the issue is the model structure. There is an issue with mxRefModels for multigroup models with inherited data. http://openmx.psyc.virginia.edu/issue/2014/12/mxrefmodels-does-not-work-inherited-data But that tends to give a different error message.
Are you able to share (perhaps privately) more about your model, script, and data? What does
print out?
It prints out the following: