5

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.

kjetil b halvorsen
  • 63,378
  • 26
  • 142
  • 467
Patrick
  • 1,381
  • 1
  • 15
  • 21
  • 1
    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](http://stats.stackexchange.com/a/22520/3277) 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. – ttnphns Jun 28 '13 at 17:32
  • 2
    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. – whuber Jun 28 '13 at 17:44
  • 1
    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`. – whuber Jun 28 '13 at 18:10
  • 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. – Patrick Jun 28 '13 at 18:40
  • @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. – ttnphns Jun 28 '13 at 19:29

1 Answers1

6

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.

kjetil b halvorsen
  • 63,378
  • 26
  • 142
  • 467
  • 3
    also `outer(S,S) * R` since `*` performs the Hadamard (elementwise) product ... – Ben Bolker May 12 '19 at 14:13
  • @Ben Bolker: It would be a good exercise to compare speedwise the efficiencies ... – kjetil b halvorsen May 12 '19 at 15:00
  • 3
    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 – whuber May 12 '19 at 15:48