r <- cor(mtcars$cyl, mtcars$gear)
r
## [1] -0.4926866In statistics, dependence or association is any statistical relationship, whether causal or not, between two random variables or bivariate data. It is often measured by the Pearson’s correlation coefficient:
\[\rho _{X,Y} =\textrm{corr} (X,Y) = \frac {\textrm{cov}(X,Y) }{ \sigma_X \sigma_Y } = \frac{ \sum_{i = 1}^n [(X - \mu_X)(Y - \mu_Y)] }{ \sigma_X \sigma_Y }\]
Other measures of correlation can be thought of, such as the Spearman \(\rho\) rank correlation coefficient or Kendall \(\tau\) rank correlation coefficient.
Assumptions for the Pearson Correlation Coefficient
The Pearson correlation coefficients makes a few assumptions, which should be carefully checked.
- Interval-level measurement: Both variables should be measured on a quantitative scale.
- Random sampling: Each subject in the sample should contribute one value on \(X\) and one value on \(Y\). The values for both variables should represent a random sample drawn from the population of interest.
- Linearity: The relationship between \(X\) and \(Y\) should be linear.
- Bivariate normal distribution: This means that (i) the values of \(X\) follow a normal distribution at each level of \(Y\), and (ii) the values of \(Y\) follow a normal distribution at each level of \(X\).
Hypothesis testing
It is possible to test whether \(r = 0\) against the alternative $ r $. The test is based on the test statistic:
\[ T = \frac{r \sqrt{n - 2}}{\sqrt{1 - r^2}}\] which follows a Student’s \(t\) distribution with \(n-2\) degrees of freedom.
Let’s take the two variables cyl and mpg from the mtcars data frame. For these variables, the Pearson’s correlation index is (please note, that the Spearman’s correlation doefficient might be preferred for these variables, although I neglected this for the sake of simplicity):
The T statistic is:
T <- r * sqrt(32 - 2) / sqrt(1 - r^2)
T
## [1] -3.101051The p-value for the null hypothesis is:
2 * pt(T, 30)
## [1] 0.004173297which is clearly highly significant. The null can be rejected.
When performing this test, it is important to note that the individuals from whom the pairs of measurements were taken must be independent. If they are not, the ‘t-test’ is invalid. I discuss this aspect in another post on my blog.
Correlation in R
We have already seen that we can use the standard cor(matrix, method = …) function. To obtain significance levels as well, we can use the rcorr() function from the ‘Hmisc’ package:
# Correlations with significance levels
library(Hmisc)
corr2 <- rcorr(as.matrix(mtcars), type="pearson")
print(corr2$r, digits = 2)
## mpg cyl disp hp drat wt qsec vs am gear carb
## mpg 1.00 -0.85 -0.85 -0.78 0.681 -0.87 0.419 0.66 0.600 0.48 -0.551
## cyl -0.85 1.00 0.90 0.83 -0.700 0.78 -0.591 -0.81 -0.523 -0.49 0.527
## disp -0.85 0.90 1.00 0.79 -0.710 0.89 -0.434 -0.71 -0.591 -0.56 0.395
## hp -0.78 0.83 0.79 1.00 -0.449 0.66 -0.708 -0.72 -0.243 -0.13 0.750
## drat 0.68 -0.70 -0.71 -0.45 1.000 -0.71 0.091 0.44 0.713 0.70 -0.091
## wt -0.87 0.78 0.89 0.66 -0.712 1.00 -0.175 -0.55 -0.692 -0.58 0.428
## qsec 0.42 -0.59 -0.43 -0.71 0.091 -0.17 1.000 0.74 -0.230 -0.21 -0.656
## vs 0.66 -0.81 -0.71 -0.72 0.440 -0.55 0.745 1.00 0.168 0.21 -0.570
## am 0.60 -0.52 -0.59 -0.24 0.713 -0.69 -0.230 0.17 1.000 0.79 0.058
## gear 0.48 -0.49 -0.56 -0.13 0.700 -0.58 -0.213 0.21 0.794 1.00 0.274
## carb -0.55 0.53 0.39 0.75 -0.091 0.43 -0.656 -0.57 0.058 0.27 1.000
print(corr2$P, digits = 2)
## mpg cyl disp hp drat wt qsec vs am
## mpg NA 6.1e-10 9.4e-10 1.8e-07 1.8e-05 1.3e-10 1.7e-02 3.4e-05 2.9e-04
## cyl 6.1e-10 NA 1.8e-12 3.5e-09 8.2e-06 1.2e-07 3.7e-04 1.8e-08 2.2e-03
## disp 9.4e-10 1.8e-12 NA 7.1e-08 5.3e-06 1.2e-11 1.3e-02 5.2e-06 3.7e-04
## hp 1.8e-07 3.5e-09 7.1e-08 NA 1.0e-02 4.1e-05 5.8e-06 2.9e-06 1.8e-01
## drat 1.8e-05 8.2e-06 5.3e-06 1.0e-02 NA 4.8e-06 6.2e-01 1.2e-02 4.7e-06
## wt 1.3e-10 1.2e-07 1.2e-11 4.1e-05 4.8e-06 NA 3.4e-01 9.8e-04 1.1e-05
## qsec 1.7e-02 3.7e-04 1.3e-02 5.8e-06 6.2e-01 3.4e-01 NA 1.0e-06 2.1e-01
## vs 3.4e-05 1.8e-08 5.2e-06 2.9e-06 1.2e-02 9.8e-04 1.0e-06 NA 3.6e-01
## am 2.9e-04 2.2e-03 3.7e-04 1.8e-01 4.7e-06 1.1e-05 2.1e-01 3.6e-01 NA
## gear 5.4e-03 4.2e-03 9.6e-04 4.9e-01 8.4e-06 4.6e-04 2.4e-01 2.6e-01 5.8e-08
## carb 1.1e-03 1.9e-03 2.5e-02 7.8e-07 6.2e-01 1.5e-02 4.5e-05 6.7e-04 7.5e-01
## gear carb
## mpg 5.4e-03 1.1e-03
## cyl 4.2e-03 1.9e-03
## disp 9.6e-04 2.5e-02
## hp 4.9e-01 7.8e-07
## drat 8.4e-06 6.2e-01
## wt 4.6e-04 1.5e-02
## qsec 2.4e-01 4.5e-05
## vs 2.6e-01 6.7e-04
## am 5.8e-08 7.5e-01
## gear NA 1.3e-01
## carb 1.3e-01 NAWe can also use these functions with two separate matrices to obtain the correlations between each column of one matrix and each column of the other:
# Correlation matrix from mtcars
x <- mtcars[1:3]
y <- mtcars[4:6]
cor(x, y)
## hp drat wt
## mpg -0.7761684 0.6811719 -0.8676594
## cyl 0.8324475 -0.6999381 0.7824958
## disp 0.7909486 -0.7102139 0.8879799Relationship to slope in linear regression
The correlation coefficient and the slope in linear regression share some similarities, as both describe how \(Y\) changes with \(X\). However, in correlation, we deal with two random variables, whereas in simple linear regression we treat \(X\) as fixed and \(Y\) as a random variable, with \(Y\) regarded as a function of \(X\) (and not the other way around).
Despite their distinct interpretations, it is useful to examine the algebraic relationship between the correlation coefficient and the regression slope. Let’s simulate a dataset with two variables drawn from a multivariate normal distribution, with means of 10 and 2, respectively, and a variance-covariance matrix of:
library(MASS)
cov <- matrix(c(2.20, 0.48, 0.48, 0.20), 2, 2)
cov
## [,1] [,2]
## [1,] 2.20 0.48
## [2,] 0.48 0.20We use the mvrnomr() function to generate the dataset.
set.seed(1234)
dataset <- data.frame( mvrnorm(n=10, mu = c(10, 2), Sigma = cov) )
names(dataset) <- c("X", "Y")
dataset
## X Y
## 1 11.756647 2.547203
## 2 9.522180 2.199740
## 3 8.341254 1.862362
## 4 13.480005 2.772031
## 5 9.428296 1.573435
## 6 9.242788 1.861756
## 7 10.817449 2.343918
## 8 10.749047 2.451999
## 9 10.780400 2.436263
## 10 11.480301 1.590436The correlation coefficient and slope are as follows:
r <- with(dataset, cor(X, Y))
b1 <- coef( lm(Y ~ X, data=dataset) )[2]
r
## [1] 0.6372927
b1
## X
## 0.1785312The formula for the slope is:
\[b_1 = \frac{ \sum_{i = 1}^n \left[ ( X-\mu_X )( Y-\mu_Y )\right] }{ \sigma^2_X } \]
From this, we see that:
\[ r = b_1 \frac{\sigma_X}{ \sigma_Y }\]
and:
\[ b_1 = r \frac{\sigma_Y}{\sigma_X}\]
Indeed:
sigmaX <- with(dataset, sd(X) )
sigmaY <- with(dataset, sd(Y) )
b1 * sigmaX / sigmaY
## X
## 0.6372927
r * sigmaY / sigmaX
## [1] 0.1785312It is also easy to see that the correlation coefficient is equivalent to the regression slope of standardized \(Y\) against standardised \(X\):
Yst <- with(dataset, scale(Y, scale=T) )
summary( lm(Yst ~ I(scale(X, scale = T) ), data = dataset) )
##
## Call:
## lm(formula = Yst ~ I(scale(X, scale = T)), data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.082006 -0.067143 -0.036850 0.009214 0.237923
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.065e-18 3.478e-02 0.000 1.0000
## I(scale(X, scale = T)) 1.785e-01 7.633e-02 2.339 0.0475 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.11 on 8 degrees of freedom
## Multiple R-squared: 0.4061, Adjusted R-squared: 0.3319
## F-statistic: 5.471 on 1 and 8 DF, p-value: 0.04748Intra-class correlation (ICC)
This is a related yet distinct concept. While Pearson’s correlation coefficient is used for bivariate data—where two different variables are observed on the same independent units—intraclass correlation arises when the same variable is repeatedly observed on the same units (e.g., across different times or spaces), so that each unit acts as a grouping factor that creates clusters of measurements.
Let’s illustrate this with an example. Consider the following data table:
#|echo: false
m1 <- c(10, 12, 13, 16)
m2 <- c(21, 20, 17, 19)
dtab <- cbind(m1, m2)
row.names(dtab) <- paste0("g", 1:4)
knitr::kable(dtab)| m1 | m2 | |
|---|---|---|
| g1 | 10 | 21 |
| g2 | 12 | 20 |
| g3 | 13 | 17 |
| g4 | 16 | 19 |
When g1, g2, g3, and g4 represent four independent units (e.g., four maize plants) on which two variables, m1 and m2 (e.g., weight and height), were measured, we are in the realm of Pearson’s correlation. On the other hand, when g1, g2, g3, and g4 are entities on which multiple measurements of the same variable (e.g., yield) are taken, we enter the realm of intraclass correlation. In this second scenario, we are not assessing joint variability between two different variables; rather, we want to determine whether belonging to the same group (‘proximity’) makes observations more similar to each other than observations from different groups.
The most common intraclass correlation coefficient is:
\[ICC = {\displaystyle {\frac {\sigma _{\alpha }^{2}}{\sigma _{\alpha }^{2}+\sigma _{\varepsilon }^{2}}}.}\]
where \(\sigma _{\alpha }^{2}\) is the between-group variance component and \(\sigma _{\varepsilon }^{2}\) is the within-group variance component (or, more clearly, the variance of observations within a group). For balanced data, those two components can be obtained from one-way ANOVA, considering that there are four ‘replicates’ per group:
y <- c(10, 12, 13, 16, 21, 20, 17, 19)
group <- factor( rep(c("A", "B"), each = 4) )
nmeasures <- 4
mod <- lm(y ~ group)
anova(mod)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## group 1 84.5 84.500 18.436 0.005128 **
## Residuals 6 27.5 4.583
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sigma2_alpha <- (84.5 - 4.583)/nmeasures
sigma2_epsilon <- 4.583
icc <- sigma2_alpha/(sigma2_alpha + sigma2_epsilon)
icc
## [1] 0.8134129In general, if subjects in one group are very similar (small \(\sigma_{\varepsilon}\)), while different groups vary widely (high \(\sigma_{\alpha}\)), the ICC is high. Conversely, if measurements are just as similar across different groups as they are within the same group, the ICC is low. In theory, the ICC cannot be negative, as variance components are expected to be positive. In practice, however, the ICC calculated via the method of moments shown above can turn out negative when the group mean square is lower than the residual mean square from the ANOVA. Such negative values are usually constrained to zero or interpreted as indicating an absence of positive clustering; in both cases, intraclass correlation remains a meaningful tool.
ICC is crucial and is most often studied using linear mixed models. But… that is a story for another day…
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: andrea.onofri@unipg.it
This post was originally published on 2019-02-07