I want to use the mvrnorm function in R, but I get an error that my sigma is negative. I created some data and used log likelihood functions to estimate the best parameters for my model. I now want to test if these parameters are significantly different from zero. I thought to do this by making use of the mvrnorm function and than do a t-test for every parameter. For the mvrnorm I used the estimated parameter as the mean and as covariance matrix the inverse of the information matrix. However I get an error while running the mvrnorm that my sigma matrix is negative. How can I now test if my parameters are significantly different from zero? Should I make the sigma matrix positive definite and if so how can I do this?

  • $\begingroup$ Let's consider a special case of your question, where the dimension is $1.$ You are asking "I'm passing a negative variance to the software, but it gives an error. What would be the best thing to do to make the variance positive?" The answer is to supply a positive variance. There's nothing else we can tell you! $\endgroup$
    – whuber
    Aug 16 '19 at 20:34
  • $\begingroup$ 'm just guessing here but I think the issue is that the numerical optimisation of the log likelihood didn't converge to an optimum but perhaps instead to a saddle point. This can easily happen if you use Newton's method. When you invert the negative of the observed fisher information matrix, you then get something which is not positive definite. Perhaps better starting values will solve the problem $\endgroup$ Aug 17 '19 at 12:40

You can simply perturb your original matrix by adding small diagonal elements, i.e. $A\leftarrow A+\mu I$, where $\mu$ is chosen a small value to guarantee the PD property. You can simply try out some values.

Another alternative can be using make.positive.definite method in package corpcor, which implements the algorithm provided here. In his blog, the author of the algorithm also provides other implementations that you can use.

  • $\begingroup$ another alternative would be to calculate the the generalized inverse $V^-$ of the vcov matrix $V$ (see stat.ethz.ch/R-manual/R-devel/library/MASS/html/ginv.html) and calculate the test statistic by hand $\hat\beta^\prime V^{-}\hat\beta\sim\chi^2_{k}$, where $\hat\beta$ is the $k\times 1$-vector estimates. $\endgroup$ Aug 16 '19 at 11:11
  • $\begingroup$ I tried the make.positive.definite method, however when I than take mvrnorm for n=1000 and let it print it results, it gives 1000 times te same value for one parameter.. Can it be that this problem comes form using make.positve.definite? $\endgroup$ Aug 16 '19 at 14:53
  • $\begingroup$ try changing the tolerance argument, i.e. ‘tol’ $\endgroup$
    – gunes
    Aug 16 '19 at 15:07
  • 1
    $\begingroup$ Although these recommendations can work mathematically, unless they are preceded by careful analysis of the cause of the original problem, I would be concerned that they are just patching over the issue and might be hiding either (1) some other serious difficulty with the analysis, (2) misconceptions, or (3) opportunities for insight. $\endgroup$
    – whuber
    Aug 17 '19 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.