9

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
  • 77,844
Patrick
  • 1,571
  • 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 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

10

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.

TRUTH= 0.8 
R <- as.matrix(data.frame(c(1, TRUTH), c(TRUTH, 1)))
S = c(sqrt(1), sqrt(1))

cor2cov_1(R,S)


outer(S,S) * R


smat = as.matrix(S)
R * smat %*% t(smat)

Here is a microbenchmark showing the efficiency of the functions:

library(microbenchmark)
microbenchmark::microbenchmark(outer(S,S) * R ,cor2cov_1(R,S), cor2cov(R,S), R * smat %*% t(smat), times = 10000)

Unit: microseconds
                 expr     min      lq       mean  median      uq      max neval cld
      outer(S, S) * R   1.968   2.214   2.724639   2.337   2.460 3611.362 10000  a 
      cor2cov_1(R, S)   1.722   1.886   2.778045   1.968   2.091 3743.259 10000  a 
        cor2cov(R, S) 113.037 116.071 125.844711 118.039 120.663 5462.020 10000   b
 R * smat %*% t(smat)   1.066   1.230   1.422712   1.435   1.517   12.177 10000  a
M. Beausoleil
  • 1,225
  • 3
  • 14
  • 28
kjetil b halvorsen
  • 77,844
  • 4
    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
  • 5
    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 – whuber May 12 '19 at 15:48
  • 1
    What is the idea behind R * smat %*% t(smat)? – M. Beausoleil Jun 08 '22 at 16:06
  • 1
    @M.Beausoleil, it's a fast way of doing the "sweep". Consider that you are multiplying each entry in the correlation matrix by the product of the component standard deviations. For example, the covariance matrix entry at [1, 2] will be the SD[1] * SD[2] * cor[1, 2]. By turning the SD vector into a matrix and multiplying it by itself (tcrossprod(S) would be even faster), you are "pregenerating" the needed 4x4 matrix of multiplications of SD entries each by each other. All that is left to do is to scale them, elementwise, by the correlations, thus the * R. Which could be first or second. – Avraham Feb 26 '24 at 09:53