Comparing correlation coefficients with R

An empirical bootstrap approach may be useful to understand the underlying concepts. This approach proceeds through the following steps:

1. Resample \(n\) individuals with replacement (\(n\) being the original sample size)
2. Compute the two correlation coefficients \(\rho^*_{X, Y_{20}}\) and \(\rho^*_{X, Y_{50}}\) using the bootstrap samples
3. Compute their difference \(d^* = \rho^*_{X, Y_{20}} - \rho^*_{X, Y_{50}}\)
4. Execute steps 1 to 3 \(B\) times (\(B > 500\) or \(1000\) may seem a reasonable choice, depending on the original sample size \(n\))
5. Use the empirical distribution of the \(B\) values of \(d^*\) to compute a 95% confidence interval for the difference between the two correlation coefficients.

This is a standard bootstrap procedure, recalled in many textbooks [2,8]. Here is an example of R implementation (I don’t use the R package boot here, so as to make the algorithm more readable and transparent):

## Set some constants:
set.seed(2020)
B <- 999 # number of bootstrap replicates
n <- nrow(dat) # total sample size
## Initialize an empty vector for bootstrap statistics:
dstar <- rep(NA, B)

## Run the bootstrap procedure:
for (b in 1:B) {
    indices <- sample(x = 1:n, size = n, replace = TRUE)
    bootsample <- dat[indices, ]
    rho1 <- cor(bootsample$X, bootsample$Y20)
    rho2 <- cor(bootsample$X, bootsample$Y50)
    dstar[b] <- rho1 - rho2
}

## Plot histogram and display confidence interval:
hist(dstar, main = "Difference in correlation coefficients",
     xlab = "d*")
dstar <- sort(dstar) # sort to compute empirical CI
abline(v = dstar[c(0.025*B, 0.975*B)], lty = 2, col = "blue")

As shown on figure 1, the 95% empirical confidence interval includes the null value. This bootstrap procedure does not provide sufficient evidence to discard the null hypothesis here, and we cannot conclude (or generalize) that \(Y_{50}\) is really a better candidate than \(Y_{20}\) to be correlated with the canine values \(X\).

There are several exact methods to provide an alternative to the resampling approach developed in the previous section. The R package cocor (whose the usage is extensively, and very clearly exposed in [1]) implements many of them. Furthermore, it is really easy-to-use, either through the command-line interface or through a dedicated website.

The main function of the package, also named cocor, can handle R formulas, as shown below. Here I ask for the results returned by two modern methods: Zou’s method [9] and Hittner’s method [4].

library(cocor)

cocor(~X + Y20 | X + Y50, data = dat,
      test = c("hittner2003", "zou2007"))

  Results of a comparison of two overlapping correlations based on dependent groups

Comparison between r.jk (X, Y20) = 0.9424 and r.jh (X, Y50) = 0.8463
Difference: r.jk - r.jh = 0.0961
Related correlation: r.kh = 0.782
Data: dat: j = X, k = Y20, h = Y50
Group size: n = 15
Null hypothesis: r.jk is equal to r.jh
Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
Alpha: 0.05

hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
  z = 1.6186, p-value = 0.1055
  Null hypothesis retained

zou2007: Zou's (2007) confidence interval
  95% confidence interval for r.jk - r.jh: -0.0217 0.3404
  Null hypothesis retained (Interval includes 0)

None of those methods can help to discard the null hypothesis at the specified alpha level (0.05), which is consistent with the bootstrap simulations.

Well, this question has received quite enlightening answers on a thread on Researchgate, and I could probably not give a better one!

A straightforward modification of the previous case apply: the same bootstrap procedure remain efficient here. Once again, we will compute a 95% confidence interval for the difference between the two correlation coefficients under study. But this time, let’s use the R package boot for a greater efficiency.

library(boot)
set.seed(2020)
## First define a function to work out the difference of coefs:
diff_corr <- function(data, indices) {
    data <- data[indices, ]
    cor1 <- cor(data$Sepal.Length, data$Sepal.Width)
    cor2 <- cor(data$Petal.Length, data$Petal.Width)
    return(cor1 - cor2)
}
## Then apply a bootstrap procedure with 999 draws:
res_boot <- boot(data = setosa,
                 R = 999,
                 statistic = diff_corr,
                 stype = "i")
## Retrieve the empirical 95% confidence interval:
boot.ci(res_boot, type = "perc", conf = 0.95)

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = res_boot, conf = 0.95, type = "perc")

Intervals : 
Level     Percentile     
95%   ( 0.1877,  0.6805 )  
Calculations and Intervals on Original Scale

This time, the null value is not a plausible value for the difference between the two correlation coefficients, thus bringing a strong evidence against the null hypothesis of no difference.

The interface formula remains very easy-to-use in this case. Let’s try with the methods zou2007 [9] and silver2004 [5]:

cocor(~Sepal.Length + Sepal.Width | Petal.Length + Petal.Width,
      data = setosa, test = c("silver2004", "zou2007"))

  Results of a comparison of two nonoverlapping correlations based on dependent groups

Comparison between r.jk (Sepal.Length, Sepal.Width) = 0.7425 and r.hm (Petal.Length, Petal.Width) = 0.3316
Difference: r.jk - r.hm = 0.4109
Related correlations: r.jh = 0.2672, r.jm = 0.2781, r.kh = 0.1777, r.km = 0.2328
Data: setosa: j = Sepal.Length, k = Sepal.Width, h = Petal.Length, m = Petal.Width
Group size: n = 50
Null hypothesis: r.jk is equal to r.hm
Alternative hypothesis: r.jk is not equal to r.hm (two-sided)
Alpha: 0.05

silver2004: Silver, Hittner, and May's (2004) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
  z = 3.0306, p-value = 0.0024
  Null hypothesis rejected

zou2007: Zou's (2007) confidence interval
  95% confidence interval for r.jk - r.hm: 0.1405 0.6986
  Null hypothesis rejected (Interval does not include 0)

The results from those two methods are consistent with the conclusion (and confidence interval) derived from our home-made bootstrap procedure, and brings a strong evidence against the null hypothesis of no difference.