mxRefModels error
Posted on
Cindy.s
Joined: 07/04/2015
Forums
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.
Be careful and an Idea
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.
# Create full and na-omitted data
fdata <- mxData(YourFullDataSet, 'raw')
odata <- mxData(na.omit(YourFullDataSet), 'raw')
# Swap na-omitted data into model
QualAceFit_rg$data <- odata
# Create ref models on the swapped model
ref <- mxRefModels(QualAceFit_rg, run = FALSE)
# Replace the na-omitted data with the full data
ref[[1]]$data <- fdata
ref[[2]]$data <- fdata
# Run the reference models
ref[[1]] <- mxRun(ref[[1]])
ref[[2]] <- mxRun(ref[[2]])
Log in or register to post comments
In reply to Be careful and an Idea by mhunter
Why not just start with means & variances?
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.)
Log in or register to post comments
In reply to Why not just start with means & variances? by neale
Fine idea
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.
Log in or register to post comments
In reply to Fine idea by mhunter
Not negative definite I think.
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).
Log in or register to post comments
In reply to Be careful and an Idea by mhunter
Thak you Michael for your
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!
Log in or register to post comments
In reply to Thak you Michael for your by Cindy.s
Hmm ... no missing data.
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
QualAceFit_rg
round(cov(QualAceFit_rg$data$observed), 3)
print out?
Log in or register to post comments
In reply to Hmm ... no missing data. by mhunter
It prints out the
It prints out the following:
>mxRefModels(QualAceFit_rg, run = F)
Error in chol.default(sampcov) :
the leading minor of order 1 is not positive definite
> round(cov(QualAceFit_rg$data$observed), 3)
Error in cov(QualAceFit_rcrg$data$observed) :
supply both 'x' and 'y' or a matrix-like 'x'
I'm happy to provide the code, it is an extension of the usual 5 group GxE model, to a 6 group case (presence of a trait or not). Since I have discordant MZ pairs also, I try not only to calculate rg, but also rc simultaneously. If the data would be also needed, I can send that privately.
Log in or register to post comments