3
$\begingroup$

I am trying to figure out how to convert a correlation matrix (R) to a covariance matrix (S) for input into a random number generator that only accepts S (rmvnorm("mvtnorm") in R)

library("mvtnorm") 

TRUTH= 0.8 # target correlation value between X1 and X2
R <- as.matrix(data.frame(c(1, TRUTH), c(TRUTH, 1)))
V <- diag(c(sqrt(1), sqrt(1))) # diagonal matrix of sqrt(variances)
S <- V %*% R %*% V
cor(rmvnorm(100, sigma=S) )

# repeat this to get an idea of the variance around Pearson's estimator

Instance where variances are not equal to 1

V <- diag(c(sqrt(3), sqrt(2))) 
S <- V %*% R %*% V
cor(rmvnorm(100, sigma=S) )

This seems to be correct, but I would like expert criticism.

$\endgroup$
  • 1
    $\begingroup$ I can't read your code (I don't know R), but conversion of corr to cov matrix or back is switching an SSCP-type matrix to a new diagonal while preserving the cosine. NewMatrix = Matrix &* (coef*t(coef)) where coef = sqrt(NewDiagonal/Diagonal), * is vector multiplication and &* is usual, elementwise multiplication. $\endgroup$ – ttnphns Jun 28 '13 at 17:32
  • 2
    $\begingroup$ The code looks ok but its comments don't. V had better be the diagonal matrix of square roots of variances for this to work. $\endgroup$ – whuber Jun 28 '13 at 17:44
  • 1
    $\begingroup$ Another way to accomplish what you want--and in some cases it might be numerically a little better--is to generate your data using the correlation matrix and post-multiply them by V. It should be obvious that this works because (1) separate linear transformations in the variables do not change their correlation but (2) rescaling a unit variance variable by a constant scales its variance by the square of that constant. Then if you look at what mvtnorm does behind the scenes to factor R, you can see how it effectively carries out the same post-multiplication by V. $\endgroup$ – whuber Jun 28 '13 at 18:10
  • $\begingroup$ Which can be seen in the 1st example where R = S, because the standard deviations are both 1. So, this has the potential to be different? I am not mathematically inclined, so I will have to "prove" this to myself via coding. Thanks again whuber. $\endgroup$ – Patrick Jun 28 '13 at 18:40
  • $\begingroup$ @Patrick: Note that your formula V %*% R %*% V is equivalent to what I've suggested above. But, I predict, your formula with two matrix multiplications is slower (which must show on big matrices). Elementwise multiplication of matrices (is there such an operation in R?) is faster, AFAIK. $\endgroup$ – ttnphns Jun 28 '13 at 19:29
2
$\begingroup$

Let $R$ be the correlation matrix and $S$ the vector of standard deviations, so that $S\cdot S$ (where $\cdot$ is the componentwise product) is the vector of variances. Then $$ \text{diag}(S) R \text{diag}(S) $$ is the covariance matrix. This is fully explained here.

This can be implemented in R as

cor2cov_1 <- function(R,S){
    diag(S) %*% R %*% diag(S)
}

but is inefficient. An efficient implementation is

cor2cov <- function(R, S) {
 sweep(sweep(R, 1, S, "*"), 2, S, "*")
 }

and you can test yourself they give the same result.

$\endgroup$
  • 2
    $\begingroup$ also outer(S,S) * R since * performs the Hadamard (elementwise) product ... $\endgroup$ – Ben Bolker May 12 at 14:13
  • $\begingroup$ @Ben Bolker: It would be a good exercise to compare speedwise the efficiencies ... $\endgroup$ – kjetil b halvorsen May 12 at 15:00
  • $\begingroup$ On my system (Microsoft's version of R), the outer solution is far faster for large dimensions. It is about eight times faster than the sweep implementation. (If you were to replace the inner sweep by R*S it would be almost twice as fast, but still four times slower than outer.) Replacing outer(S,S) by d <- length(S); (matrix(S,d,1) %*% matrix(S,1,d)) is a little faster. cc @BenBolker $\endgroup$ – whuber May 12 at 15:48

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.