# How to test for a significant difference between several dependent correlation coefficients?

I have a group of correlation coefficients (more than two). They are all dependent on one variable A in the form of r_A1, r_A2, r_A3....r_Ak, where 1, 2 ...k denote other variables; they all have the same sample size.

My question is: what statistics should I use when I want to know whether any one of the correlation coefficients is any different to any of the others? I know if there are only two dependent correlation coefficients, this can be easily compared using most statistical tools, but this is a multiple correlations test. I know that a Chi-square test can be used to compare equality of several correlation coefficients (see http://home.ubalt.edu/ntsbarsh/business-stat/otherapplets/MultiCorr.htm for an example), but to my knowledge this approach is for testing the difference between INDEPENDENT correlation coefficients. So I am wondering whether there is any approach that is equivalent to Fisher's least significant difference that can be used to make comparisons among several dependent correlation coefficients?

EDIT: thanks @russ-lenth for your answer. In general, I found that CIs computed from lsmeans are larger than those computed using Fisher's Z method. Here's an example of the CIs that I get through the lsmeans function:

 rep.meas      lsmean         SE df   lower.CL  upper.CL
M1        0.76914236 0.13325688 23  0.4934795 1.0448052
M2        0.82346705 0.11830361 23  0.5787374 1.0681967
M3        0.89294217 0.09386717 23  0.6987631 1.0871212
M4       -0.09985512 0.20747224 23 -0.5290441 0.3293339
M5        0.56183690 0.17249315 23  0.2050076 0.9186662
M6        0.79086279 0.12760947 23  0.5268825 1.0548431
M7        0.14667681 0.20625924 23 -0.2800029 0.5733566


Take M1 whose r = 0.769 as an example: the width of the CI from lsmeans is (1.0448-0.4935) 0.5513. The width of the CI computed from Fisher’s Z is (0.8948-0.5302) 0.3646, which is much smaller than the former. Is the difference between the widths of the two confidence intervals too large?

• Are you sure this is the question you want answered? What is the practical meaning if you find that one correlation is higher than two others? Is your goal in fact to select variables for a regression equation? Are in fact all these variables in the same dataset? If not, maybe they are not dependent. Anyway, this seems like a difficult technical problem, so maybe it's worth spending some time making sure it's a relevant question in the context of your research goals. – rvl Sep 1 '14 at 21:40
• Hi, thanks for your comment. Basically I am testing the performance of several models. The performance of each model is evaluated by correlating the outputs of the model with human performance data. The higher the correlation is, the better the model. Since the outputs of all the models have been correlated with the same subjective data, I think the comparison of coefficients should therefore be dependent in this case? – ksna Sep 1 '14 at 22:32

Well, here's an idea...

First standardize all of the variables (the model outputs as well as the reference variable).

Then fit a multivariate model with those standardized model outputs as a multivariate response variable (not predictors), and the common human-performance variable (standardized) as the predictor. Do not include an intercept, as it will be zero anyway. Then the regression coefficients will be equal to the correlation coefficients, due to the standardization, and their covariance matrix will be available; so you can estimate each pairwise difference and its standard error.

### R example

This is using the swiss dataset provided with R. Here is the standardized dataset

> swiss.std = as.data.frame(lapply(swiss, function(x) (x-mean(x))/sd(x)))


Note the covariances are just the correlations

> cov(swiss.std[1:4])
Fertility Agriculture Examination  Education
Fertility    1.0000000   0.3530792  -0.6458827 -0.6637889
Agriculture  0.3530792   1.0000000  -0.6865422 -0.6395225
Examination -0.6458827  -0.6865422   1.0000000  0.6984153
Education   -0.6637889  -0.6395225   0.6984153  1.0000000


Fit the multivariate model, looking at correlations with Fertility

> swiss.mlm = lm(cbind(Agriculture,Examination,Education) ~ Fertility - 1, data = swiss.std)


Here are the coefficients and variance-covariance matrix thereof

> coef(swiss.mlm)
Agriculture Examination  Education
Fertility   0.3530792  -0.6458827 -0.6637889

> vcov(swiss.mlm)
Agriculture:Fertility Examination:Fertility Education:Fertility
Agriculture:Fertility           0.019029024          -0.009967271        -0.008807663
Examination:Fertility          -0.009967271           0.012670338         0.005862729
Education:Fertility            -0.008807663           0.005862729         0.012160529


So to compare, say, the 2nd and 3rd correlation, here's the estimate

> con = c(0,1,-1)
> coef(swiss.mlm) %*% con
[,1]
Fertility 0.01790615


and its SE

> sqrt(sum(con * vcov(swiss.mlm) %*% con))
[1] 0.1144789


I can trick the lsmeans package into doing it:

> library(lsmeans)
> swiss.lsm = lsmeans(swiss.mlm, "rep.meas")


rep.meas is the default name for the levels of the multivariate response. So far, swiss.lsm is just estimating the mean, $(0,0,0,)$, but I'll change the linear function that it's using to be $1$ times each regression coefficient

> swiss.lsm@linfct = diag(c(1,1,1))


Now, here is the summary

> swiss.lsm
rep.meas        lsmean        SE df    lower.CL   upper.CL
Agriculture  0.3530792 0.1379457 46  0.07540884  0.6307495
Examination -0.6458827 0.1125626 46 -0.87245946 -0.4193060
Education   -0.6637889 0.1102748 46 -0.88576050 -0.4418172


and the pairwise comparisons:

> pairs(swiss.lsm)
contrast                    estimate        SE df t.ratio p.value
Agriculture - Examination 0.99896189 0.2272309 46   4.396  0.0002
Agriculture - Education   1.01686804 0.2209183 46   4.603  0.0001
Examination - Education   0.01790615 0.1144789 46   0.156  0.9866

P value adjustment: tukey method for a family of 3 means


If I want to compare the absolute correlations, just change the linear function

> swiss.lsm@linfct = diag(c(1,-1,-1))
> swiss.lsm
rep.meas       lsmean        SE df   lower.CL  upper.CL
Agriculture 0.3530792 0.1379457 46 0.07540884 0.6307495
Examination 0.6458827 0.1125626 46 0.41930596 0.8724595
Education   0.6637889 0.1102748 46 0.44181722 0.8857605

Confidence level used: 0.95

> pairs(swiss.lsm)
contrast                     estimate        SE df t.ratio p.value
Agriculture - Examination -0.29280352 0.1084658 46  -2.700  0.0257
Agriculture - Education   -0.31070967 0.1165085 46  -2.667  0.0279
Examination - Education   -0.01790615 0.1144789 46  -0.156  0.9866

P value adjustment: tukey method for a family of 3 means

• If you use R, I can show an example – rvl Sep 1 '14 at 23:58
• Thank you very much for your answer, I will try it out. Yes it would be really useful if you have an example of it using R, thanks! – ksna Sep 2 '14 at 9:41
• OK, I added an example. Sorry it took so long - busy morning – rvl Sep 2 '14 at 19:45
• That's perfect, thank you so much! Such a detailed example was really really useful and clear. I've spent so much time on this and haven't managed to work it out, but now it works - many thanks!! – ksna Sep 2 '14 at 21:34
• There are two issues here. First of all, the sampling distribution of correlation coefficients is not normal when the true correlation is large and $n$ is small (which is the case here). Second, when you standardize two variables and then put them into a regression model, then the standard error you get for the coefficient (i.e., correlation) is only consistent when the true correlation is zero (see: Yuan, K.-H., & Chan, W. (2011). Biases and standard errors of standardized regression coefficients. Psychometrika, 76(4), 670-690.). So, the CIs you get with this approach are not correct. – Wolfgang Sep 5 '14 at 9:42

There is a rich literature on testing dependent correlation coefficients. Here are some relevant references:

Dunn, O. J., & Clark, V. (1971). Comparison of tests of the equality of dependent correlation coefficients. Journal of the American Statistical Association, 66(336), 904-908.

Neill, J. J., & Dunn, O. J. (1975). Equality of dependent correlation coefficients. Biometrics, 31(2), 531-543.

Steiger, J. H. (1980). Tests for comparing elements of a correlation matrix. Psychological Bulletin, 87(2), 245-251.

Steiger, J. H. (1980). Testing pattern hypotheses on correlation matrices: Alternative statistics and some empirical results. Multivariate Behavioral Research, 15(3), 335-352.

Olkin, I., & Finn, J. D. (1990). Testing correlated correlations. Psychological Bulletin, 108(2), 330-333.

Meng, X., Rosenthal, R., & Rubin, D. B. (1992). Comparing correlated correlation coefficients. Psychological Bulletin, 111(1), 172-175.

Raghunathan, T. E., Rosenthal, R., & Rubin, D. B. (1996). Comparing correlated but nonoverlapping correlations. Psychological Methods, 1(1), 178-183.

Cheung, M. W.-L., & Chan, W. (2004). Testing dependent correlation coefficients via structural equation modeling. Organizational Research Methods, 7(2), 206-223.

You will find equations for the covariance of overlapping correlations in these references (for raw and r-to-z transformed correlations). You can then approach this problem from a meta-analytic perspective. In essence, you have $k$ correlations (raw or r-to-z transformed) with an approximately known variance-covariance matrix that can be computed/estimated. Let $y$ denote the (column) vector with the correlations and $V$ the corresponding var-cov matrix. Then you can test the null hypothesis that all true correlations are the same (homogeneous) by comparing the test statistic $$Q = y'Py$$ against the critical value of a chi-square distribution with $k-1$ degrees of freedom, where $P = W - WX(X'WX)^{-1}X'W$, $X$ is just a column of 1s, and $W = V^{-1}$.

You can also test more focused hypotheses, such as that the $i$th correlation is different from the rest. For this, let $X$ be a $k \times 2$ matrix with 1s in the first column and 0s in the second column, except for $i$th row, which gets a 1 in the second column. Then you can fit the model $$y = X\beta + e,$$ where $e \sim N(0, V)$. You can estimate $\beta = [\beta_0, \beta_1]'$ with $$b = (X'WX)^{-1} X'Wy$$ and the variance-covariance matrix of $b$ with $$Var(b) = (X'WX)^{-1}.$$ Now compute $z = b_1 / SE[b_1]$, where $b_1$ is the second element from $b$ and $SE[b_1]$ is the standard error of $b_1$ (i.e., the square root of the second diagonal element from $Var(b)$). You can compare $z$ against the critical values of a standard normal distribution (i.e., $\pm 1.96$ for $\alpha = .05$, two-sided) to test whether the $i$th correlation is different from the rest.

• Thank you very much for your answer - those papers were really useful! – ksna Sep 2 '14 at 21:38
• @Wolfgang, I'll look at those papers when I get a chance. I imagine one of them (older ones especially) must come pretty close to what I did in my answer. – rvl Sep 2 '14 at 22:27