I want to estimate a covariance matrix from data with some missing values. Ideally I'd like an R package but python could be ok.

R has some built in ways of doing this. You can use


Or the same using cor (correlation). The trouble is that if you do this with cov, the matrix will not be guaranteed to be positive definite. If you do cov2cor(cor.mat), you will find correlation coefficients outside of [-1,1]. Using pairwise with cor seems to handle this. Then I could use the diagonal variances to go from cor.mat to cov.mat. Still, this is probably not optimal.

There appears to be a few packages that claim to do this (mvnmle, rsem) but neither appear to work. rsem fails to run for me. mvnmle can only handle up to 50 variables. I need to handle roughly 1500 variables. Would like it to run in a few seconds.

Anyone know of a good package for this?

  • $\begingroup$ What error do you get from rsem? (It might be running out of memory). You could try imputation, but with 1500 variables it's not going to take a few seconds, and it's going to take a big chunk of memory. What's the sample size? $\endgroup$ Aug 5, 2014 at 16:49
  • $\begingroup$ Error in solve.default(sig0) : system is computationally singular: reciprocal condition number = 1.938e-21 $\endgroup$
    – Dave31415
    Aug 5, 2014 at 20:09
  • 1
    $\begingroup$ This is a large p, small n problem if that matters. $\endgroup$
    – Dave31415
    Aug 5, 2014 at 20:10
  • 1
    $\begingroup$ Oh, that makes sense. rsem expects a positive definite matrix, and if n < p, it won't be. $\endgroup$ Aug 5, 2014 at 22:55
  • 1
    $\begingroup$ I asked a very closely related question at stats.stackexchange.com/questions/1781, which concerns missingness due to censoring. If your form of missingness is completely at random, it could be viewed as complete censoring, whence any useful answer to my question (should one ever appear!) ought to apply here too. $\endgroup$
    – whuber
    Jan 21, 2015 at 22:40

1 Answer 1


Here are some of the approaches that come to mind.

You can use nearPD function from the Matrix package to convert the matrix output from cov/cor to a positive definite matrix.

You can use an imputation package such as VIM to fill missing data and then use cov.

You can code up your own routine to compute cov estimate using Expectation Maximization under normality. For details you can refer to Rubin.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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