A Trip from Covariance to Correlation and Back
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
When we have a multi-trait table of data, we might be interested in investigating variability structure and the reciprocal relationships (whether causal or non-causal) among variables. For example, let’s consider the ‘mtcars’ dataset, shown below:
data(mtcars) head(mtcars) ## mpg cyl disp hp drat wt qsec vs am gear carb ## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 ## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 ## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 ## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 ## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 ## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
Considering the different cars (along the rows), we might be interested in investigating the structure of variability for fuel consumption (mpg), number of cylinders (cyl), displacement (disp) and horsepower (hp). How much does each variable change across the different cars? How much does each variable change when another variable changes? Is there any kind of joint variation, either positive or negative?
One suitable statistic to answer these question is the covariance, that is:
where and are the means for each variable and is the number of individuals. In words, covariance is the sum of products of the ith deviations from the mean for both variables. In this setting, the variance is just the covariance of a variable with itself.
For a data table, variances and covariances can be calculated by using the cov() function, that returns a matrix (the variance-covariance matrix, indeed):
matr <- mtcars[,1:4] # Covariances Sigma <- cov(matr) Sigma ## mpg cyl disp hp ## mpg 36.324103 -9.172379 -633.0972 -320.7321 ## cyl -9.172379 3.189516 199.6603 101.9315 ## disp -633.097208 199.660282 15360.7998 6721.1587 ## hp -320.732056 101.931452 6721.1587 4700.8669
While covariances indicate direction (positive or negative association), their values depend on the scale of measurement, making it difficult to gauge the strength of the relationship at a glance. For this reason, covariances are usually ‘standardized’ by dividing by the product of the standard deviations to obtain Pearson’s correlation coefficient:
where and are the standard deviations for and . Unlike covariances, correlations go from -1 to 1 and the correlation matrix has 1s along its main diagonal (since the correlation of any variable with itself is 1). In R, the matrix of correlation is obtained by using the cor() function:
# Correlations R <- cor(matr) R ## mpg cyl disp hp ## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## cyl -0.8521620 1.0000000 0.9020329 0.8324475 ## disp -0.8475514 0.9020329 1.0000000 0.7909486 ## hp -0.7761684 0.8324475 0.7909486 1.0000000
This brings us to the main question of this post: Given a variance-covariance matrix, can we compute the correlation matrix? And can we do the reverse—moving from a correlation matrix back to a variance-covariance matrix? With reference to the matrices above, can we go from ‘Sigma’ to ‘R’ and from ‘R’ to ‘Sigma’? I often find myself needing to make these conversions, especially when working with mixed models or time series analysis.
From covariance to correlation
Considering the previous expression, it is easy to see that, if we take the covariance between ‘disp’ and ‘mpg’ above (-633.0972), the correlation can be obtained as:
-633.097208 / (sqrt(36.324103) * sqrt(15360.7998)) ## [1] -0.8475514
and all the necessary information is contained in the variance-covariance matrix.
Of course, doing these calculations for all covariances one by one is tedious. The question is: is there any simple way to go from a covariance matrix to a correlation matrix with one single operation?
The answer is found in matrix algebra; let’s take a look at the variance-covariane matrix: what we have to do is to take each element in this matrix and divide it by the square roots of the diagonal elements in the same column and in the same row (see figure below). Consequently, our task is:
- take the diagonal elements
- take each ith row and divide all its elements by the ith diagonal element (row-scaling)
- take each ith column of the resulting matrix and divide all its elements by the ith diagonal element (column-scaling)
Row-scaling can be performed by using matrix-algebra; in particular, we can take a diagonal matrix of weights and multiply it by the square matrix we want to scale, as shown below.
We see that all elements of the first row were multiplied by the first weighting element, all elements of the second row were multiplied by the second weighting element and so forth.
Column scaling is done similarly, but in reverse order: we take the target matrix and post-multiply it by the diagonal matrix of weights, as shown below.
We see that all elements of the first column were multiplied by the first weighting element, all elements of the second row were multiplied by the second weighting element, and so forth.
Consequently, we can take our variance-covariance matrix and:
- Take the square roots of the diagonal elements (standard deviations) and load them into a diagonal matrix
- Invert this diagonal matrix
- Pre-multiply and post-multiply the variance-covariance matrix by this diagonal matrix of inverted standard deviations
# from covariance to correlation - 1 StDev <- sqrt(diag(Sigma)) StDevMat <- diag(StDev) InvStDev <- solve(StDevMat) InvStDev %*% Sigma %*% InvStDev ## [,1] [,2] [,3] [,4] ## [1,] 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## [2,] -0.8521620 1.0000000 0.9020329 0.8324475 ## [3,] -0.8475514 0.9020329 1.0000000 0.7909486 ## [4,] -0.7761684 0.8324475 0.7909486 1.0000000
The previous solution is correct, but, in R, it is rather inefficient. The best way is to use the cov2cor() function in the ‘stat’ package.
# From covariance to correlation - 2 cov2cor(Sigma) ## mpg cyl disp hp ## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## cyl -0.8521620 1.0000000 0.9020329 0.8324475 ## disp -0.8475514 0.9020329 1.0000000 0.7909486 ## hp -0.7761684 0.8324475 0.7909486 1.0000000
With base R, we can also sweep() twice:
# From covariance to correlation - 3 sweep(sweep(Sigma, 1, StDev, FUN = "/"), 2, StDev, FUN = "/") ## mpg cyl disp hp ## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## cyl -0.8521620 1.0000000 0.9020329 0.8324475 ## disp -0.8475514 0.9020329 1.0000000 0.7909486 ## hp -0.7761684 0.8324475 0.7909486 1.0000000
Alternatively, we can combine scale() and t() (transpose) twice, though it is noticeably less intuitive:
# From covariance to correlation - 4
scale(t(scale(t(Sigma), center = F, scale = StDev)),
center = F, scale = StDev)
## mpg cyl disp hp
## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl -0.8521620 1.0000000 0.9020329 0.8324475
## disp -0.8475514 0.9020329 1.0000000 0.7909486
## hp -0.7761684 0.8324475 0.7909486 1.0000000
## attr(,"scaled:scale")
## mpg cyl disp hp
## 6.026948 1.785922 123.938694 68.562868
Another very clean approach in R is based on the outer product. If we have a vectors of standard deviations ($), the outer product is:
Thus, the correlation is calculated as:
# From covariance to correlation - 5 Sigma/outer(StDev, StDev) ## mpg cyl disp hp ## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## cyl -0.8521620 1.0000000 0.9020329 0.8324475 ## disp -0.8475514 0.9020329 1.0000000 0.7909486 ## hp -0.7761684 0.8324475 0.7909486 1.0000000
From correlation to covariance
It is easy to see that:
Thus, if we have the correlation between ‘disp’ and ‘mpg’ (-0.8521620), we can easily calculate the covariance as:
-0.8475514 * sqrt(36.324103) * sqrt(15360.7998) ## [1] -633.0972
The only problem is that the variances are not contained in a correlation matrix, so we need to get them from elsewhere. Once we have them, we
- Take their square roots and place them into a diagonal matrix
- Pre-multiply and post-multiply the correlation matrix for this diagonal matrix of standard deviations
The solution is as follows:
# from correlation to covariance - 1 StDev <- sqrt(diag(Sigma)) StDevMat <- diag(StDev) StDevMat %*% R %*% StDevMat ## [,1] [,2] [,3] [,4] ## [1,] 36.324103 -9.172379 -633.0972 -320.7321 ## [2,] -9.172379 3.189516 199.6603 101.9315 ## [3,] -633.097208 199.660282 15360.7998 6721.1587 ## [4,] -320.732056 101.931452 6721.1587 4700.8669
More efficiently, we can use the cor2cov() function in the ‘lavaan’ package.
# From correlation to covariance - 2 lavaan::cor2cov(R, StDev) ## mpg cyl disp hp ## mpg 36.324103 -9.172379 -633.0972 -320.7321 ## cyl -9.172379 3.189516 199.6603 101.9315 ## disp -633.097208 199.660282 15360.7998 6721.1587 ## hp -320.732056 101.931452 6721.1587 4700.8669
Or, we can sweep() twice:
# From correlation to covariance - 3 sweep(sweep(R, 1, StDev, FUN = "*"), 2, StDev, FUN = "*") ## mpg cyl disp hp ## mpg 36.324103 -9.172379 -633.0972 -320.7321 ## cyl -9.172379 3.189516 199.6603 101.9315 ## disp -633.097208 199.660282 15360.7998 6721.1587 ## hp -320.732056 101.931452 6721.1587 4700.8669
or scale() and t() twice:
# From correlation to covariance - 5
scale(t(scale(t(R), center = F, scale = 1/StDev)),
center = F, scale = 1/StDev)
## mpg cyl disp hp
## mpg 36.324103 -9.172379 -633.0972 -320.7321
## cyl -9.172379 3.189516 199.6603 101.9315
## disp -633.097208 199.660282 15360.7998 6721.1587
## hp -320.732056 101.931452 6721.1587 4700.8669
## attr(,"scaled:scale")
## mpg cyl disp hp
## 0.165921457 0.559934979 0.008068505 0.014585154
Finally, we can use the outer product as follows:
# From correlation to covariance - 6 R * outer(StDev, StDev) ## mpg cyl disp hp ## mpg 36.324103 -9.172379 -633.0972 -320.7321 ## cyl -9.172379 3.189516 199.6603 101.9315 ## disp -633.097208 199.660282 15360.7998 6721.1587 ## hp -320.732056 101.931452 6721.1587 4700.8669
Do you have a cleaner or faster trick for this in R? Drop a comment below and I’ll be happy to feature it here!
Thanks for reading, and happy coding!
(P.S. If you enjoyed this post, check out my book “Field Research Methods in Agriculture: An Introduction with R” via the link below!)
Prof. Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Send comments to: [email protected]
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.