# Why is the sum of individual Spearman's rho squared less than 1 as opposed to Pearson's r in a synthetic example?

A relatively low number of iid random vectors of a relatively high dimension (10,000) are added up together element wise: $$\sum_{i=1}^{n}X_i=Y$$ where $$dim(X_i)=dim(X_j)=dim(Y),\forall i,j$$ and $$dim(X_i)\gg n$$

My modeling suggests that $$E \left[\sum_{i=1}^{n}\rho_{X_i,Y}^2\right]<1$$ for the described case.

More generally $$E \left[\sum_{i=1}^{n}\rho_{X_i,Y}^2\right] and $$E \left[\sum_{i=1}^{n}r_{X_i,Y}^2\right]\ge1$$ where $$\rho$$ is Spearman's rank correlation coefficient, $$r$$ is Pearson correlation coefficient.

# Computational assumptions

For computing $$\rho_{X,Y}^2$$ and $$r_{X,Y}^2$$ in R I use

cor(X, Y, method='spearman')^2


and

cor(X, Y, method='pearson')^2


respectively.

I also assume $$X_i \sim \mathcal{N}(0, 1)$$ but the following works for some other common distributions as well.

# Fixed number of variables, fixed dimension, varied correlation function

Given just 3 of 10,000-vectors the distributions of Spearman's $$\rho_{X_1,Y}^2+\rho_{X_2,Y}^2+\rho_{X_3,Y}^2$$ and Pearson's $$r_{X_1,Y}^2+r_{X_2,Y}^2+r_{X_3,Y}^2$$ are as follows:

# Fixed number of variables, varied dimension

Here's how the distributions of $$\rho_{X_1,Y}^2+\rho_{X_2,Y}^2+\rho_{X_3,Y}^2$$ and $$r_{X_1,Y}^2+r_{X_2,Y}^2+r_{X_3,Y}^2$$ look like next to each other for various dimensions. It appears as with the increase of the number of dimensions both distributions travel to the left but Pearson's stops at 1.0 while Spearman's continues to move below 1.0.

# Fixed dimension, varied number of variables

When the number of variables increases, the Spearman's distribution seems to drift even further below 1.0 in the beginning and then comes back and exceeds 1.0 whereas Pearson's doesn't go below 1.0 at all and bounces back at the high number of variables just as Spearman's does:

# Varied dimension, varied number of variables

A 3D plot to corroborate the above.

# Questions

1. Is my modeling correct?
2. Why is the Pearson sum of squares always greater than Spearman's on average?
3. Why does the Pearson sum of squares never get below 1.0 on average whereas Spearman's does?

Here's my code:

library(ggplot2)
library(dplyr)
library(purrr)
library(plotly)

r_squared_sum <- function(num_vars, distrib_func, rsq_func, iterations=1000) {
R_sq_sums <- c()
for (i in 1:iterations) {
data <- c()
for (j in 1:num_vars) {
data <- append(data, list(distrib_func()))
}
Y = Reduce(+, data)
R_sq_sum = Reduce(+, lapply(data, partial(rsq_func, y=Y)))

R_sq_sums <- append(R_sq_sums, R_sq_sum)
}
return(data.frame(r_squared_sum=R_sq_sums))
}

# ---------------------------------------------
# Fixed number of variables, fixed dimensions,
# varied correlation function, histogram
# ---------------------------------------------
data_spearman <- r_squared_sum(
num_vars=3,
distrib_func=partial(rnorm, 10000, 0, 1),
rsq_func=function (x, y) cor(x, y, method='spearman')^2
)
data_spearman$method <- 'spearman' data_pearson <- r_squared_sum( num_vars=3, distrib_func=partial(rnorm, 10000, 0, 1), rsq_func=function (x, y) cor(x, y, method='pearson')^2 ) data_pearson$method <- 'pearson'

ggplot(
bind_rows(data_pearson, data_spearman),
aes(x=r_squared_sum, fill=method, color=method)
) + geom_histogram(position="identity", alpha=0.8)

# -------------------------------------------------------
# Fixed number of variables, varied dimension, boxplots
# -------------------------------------------------------
dimension_range <- c(10, 100, 1000, 10000)
distrib_func_range <- lapply(
dimension_range,
partial,
.f=partial(rnorm, mean=0, sd=1)
)
data_frames <- lapply(
distrib_func_range,
partial(
r_squared_sum,
num_vars=3,
rsq_func=function (x, y) cor(x, y, method='spearman')^2
)
)
for (i in 1:length(dimension_range)) {
data_frames[[i]]$$dimension <- as.character(dimension_range[i]) data_frames[[i]]$$method <- 'spearman'
}
data_spearman <- bind_rows(data_frames)

data_frames <- lapply(
distrib_func_range,
partial(
r_squared_sum,
num_vars=3,
rsq_func=function (x, y) cor(x, y, method='pearson')^2
)
)
for (i in 1:length(dimension_range)) {
data_frames[[i]]$$dimension <- as.character(dimension_range[i]) data_frames[[i]]$$method <- 'pearson'
}
data_pearson <- bind_rows(data_frames)

ggplot(
bind_rows(data_pearson, data_spearman),
aes(y=dimension, x=r_squared_sum, fill=method, color=method)
) + geom_boxplot()

# -------------------------------------------------------
# Fixed dimension, varied number of variables, boxplots
# -------------------------------------------------------
num_vars_range = c(10, 100, 1000)
data_frames <- lapply(
num_vars_range,
partial(
r_squared_sum,
distrib_func=partial(rnorm, 10000, 0, 1),
rsq_func=function (x, y) cor(x, y, method='spearman')^2
)
)
for (i in 1:length(num_vars_range)) {
data_frames[[i]]$$num_vars = as.character(num_vars_range[i]) data_frames[[i]]$$method <- 'spearman'
}
data_spearman <- bind_rows(data_frames)

data_frames <- lapply(
num_vars_range,
partial(
r_squared_sum,
distrib_func=partial(rnorm, 10000, 0, 1),
rsq_func=function (x, y) cor(x, y, method='pearson')^2
)
)
for (i in 1:length(num_vars_range)) {
data_frames[[i]]$$num_vars = as.character(num_vars_range[i]) data_frames[[i]]$$method <- 'pearson'
}
data_pearson <- bind_rows(data_frames)

ggplot(
bind_rows(data_pearson, data_spearman),
aes(y=num_vars, x=r_squared_sum, fill=method, color=method)
) + geom_boxplot()

# ---------------------------------------------------------
# Varied dimension, varied number of variables, 3D surfaces
# ---------------------------------------------------------
num_vars_range = seq(10, 100, by=10)
dimension_range <- seq(100, 1000, by=100)
distrib_func_range <- lapply(
dimension_range,
partial,
.f=partial(rnorm, mean=0, sd=1)
)
nrow <- length(num_vars_range)
ncol <- length(distrib_func_range)

data_pearson <- matrix(data=NA, nrow=nrow, ncol=ncol)
for (i in 1:nrow) {
for (j in 1:ncol) {
r_squared_sums <- r_squared_sum(
num_vars_range[[i]],
distrib_func_range[[j]],
function (x, y) cor(x, y, method='pearson')^2,
iterations=100
)
data_pearson[i, j] <- mean(r_squared_sums$r_squared_sum) } } data_spearman <- matrix(data=NA, nrow=nrow, ncol=ncol) for (i in 1:nrow) { for (j in 1:ncol) { r_squared_sums <- r_squared_sum( num_vars_range[[i]], distrib_func_range[[j]], function (x, y) cor(x, y, method='spearman')^2, iterations=100 ) data_spearman[i, j] <- mean(r_squared_sums$r_squared_sum)
}
}

fig <- plot_ly(
x=dimension_range,
y=num_vars_range,
showscale=FALSE,
showlegend=TRUE,
name='sums of r^2',
width=700,
height=500
) %>% layout(
scene=list(
xaxis=list(title="vectors"),
yaxis=list(title="dimensions"),
zaxis=list(title="mean sum r_sq")
)
)
z=data_spearman,
name='spearman',
colorscale=list(c(0, 1), c("#00BFC4","#00BFC4")),
hovertemplate=paste(' vectors: %{x}<br>', 'dimensions: %{y}<br>', 'mean sum r_sq: %{z}')
)
z=data_pearson,
name='pearson',
colorscale=list(c(0, 1), c("#F8776D","#F8776D")),
hovertemplate=paste(' vectors: %{x}<br>', 'dimensions: %{y}<br>', 'mean sum r_sq: %{z}')
)
z=matrix(data=1, nrow=nrow, ncol=ncol),
name='level',
colorscale=list(c(0, 1), c("#696969","#696969")),
hovertemplate=paste(' vectors: %{x}<br>', 'dimensions: %{y}<br>', 'mean sum r_sq: %{z}')
)
fig


I've implemented similar modeling in Python (my main language) with the same results.

• Your code does not appear to study the situation you have described mathematically. Your equations concern the correlations between individual $X_i$ and the final $Y,$ which we ought to more explicitly write as $Y_n = X_1+X_2+\cdots+X_n.$ Is this what you really intend??
– whuber
Commented Jun 30, 2022 at 13:04
• I'm not sure I completely understand your concern @whuber but I'm adding up a finite number (say $n$) of iid $X_i$s to get $Y$ and then try to measure Spearman's rho squared between each $X_i$ and $Y$ which ideally should all be the same and equal to $\frac{1}{n}$ each. Does it make it clearer? Commented Jun 30, 2022 at 13:15
• That helps a little, thanks. But how does your code assess this? It seems you have hard-coded the value of $n=3.$ Unless I am mistaken in my reading of the code, all you're doing is using Monte-Carlo to estimate a quantity you could compute analytically; and you aren't studying what happens as $n$ grows large, as suggested by the language you use.
– whuber
Commented Jun 30, 2022 at 14:37
• Thanks for pointing out the caveats @whuber. I've rewritten the post. Commented Jul 2, 2022 at 16:03
• That code differs from what you have expressed mathematically. You aren't evaluating the expectation of the squared correlation of two $3$-vectors: you are evaluating the squared correlation of two $10000$-vectors. When "10000" is made less, the sum of squared correlations often exceeds $1.0.$
– whuber
Commented Jul 3, 2022 at 18:08

TLDR; the setting (some sum of variables) where you observe a different Spearman's rank correlation and Pearson correlation can be greatly simplified.

If you use $$X$$ as normal distributed, then your question is equivalent to the case of bivariate normal distribution for which the Spearman's rank correlation is lower than the Pearson correlation.

For that case the Spearman's rank correlation can be computed as

$$r_{x,y} = \frac{6}{\pi} \sin^{-1}(\rho/2)$$

See the question: Finding correlation between CDF of two normal distributions

### Simplifying the expression

The expectation of a sum equals the expectation of the components.

So we can also express this as a comparison for only one of the components. Your problem is equivalent to:

$$E[\rho_{X_1,Y}^2] < 1/n$$ $$E[r_{X_1,Y}^2] = 1/n$$

We can simplify this further. This pair of variables $${X_1,Y}$$, with $$Y$$ as a sum of $$n$$ terms is effectively the correlation with a sum of two variables instead of $$n$$ variables.

$$Y = \sum_{i=1}^n X_i = X_1 + Z \\ \qquad \text{with } Z = \sum_{i=2}^n X_i \sim N(0,\sqrt{n-1})$$

### The correlation of the distribution

The Pearson correlation of the distribution for the population we can compute exactly

$$r_{X_1,X_1+Z}^2 = \frac{\text{Cov}(X_1,X_1+Z)^2}{{\text{Var}(X_1)\text{Var}(X_1+Z)}}$$

and we can express the covariance of a sum as $$\text{Cov}(A,B+C) = \text{Cov}(A,C) + \text{Cov}(A,B)$$ which in your specific case becomes $$\text{Cov}(X_1,X_1+Z) = \text{Cov}(X_1,Z) + \text{Cov}(X_1,X_1) = 0 + \text{Var}(X_1)$$ If we fill that in then we get to

$$r_{X_1,X_1+Z}^2 = \frac{\text{Var}(X_1)^2}{{\text{Var}(X_1)\text{Var}(X_1+Z)}}$$

We also have

$$\text{Var}(X_1+Z) = \sum_{i=1}^{n}\text{Var}(X_i) = n \text{Var}(X_1)$$

where the latter equality is because all $$\text{Var}(X_i)$$ have the same value of the $$X_i$$ are similarly distributed.

And

$$r_{X_1,X_1+Z}^2 = \left(\frac{\text{Var}(X_1)}{\sqrt{\text{Var}(X_1)n \text{Var}(X_1)}}\right)^2 = \frac{1}{{n}}$$

This is the correlation for the population of $$X_1$$ and $$Y$$. The expectation of the sample correlation will approach towards this value of we increase the sample size $$n$$. It will not be exactly the same.

In the case that the $$X_i$$ are normal distributed then the exact distribution for the sample correlation is this monster. I could not easily find computations of the expectation value of the observed squared correlation coefficient, so I did a simple computation and it might be that we have the following relationship

$$E[\hat{r}^2] = \frac{n}{n-1} r^2$$

where $$n$$ is the sample size. The computation below is a simulation that verifies that this relationship might be true.

    set.seed(1)

p = 3
n = 10

corsample = function() {
X = matrix(rnorm(p*n),n)
X1 = X[,1]
Y = rowSums(X)
cor(X1,Y)^2
}

rho = replicate(10^3,corsample())
mean(rho)
(1/3)*n/(n-1)


### A way to approximate the expectation value of the sample correlation for large samples.

The sample Pearson correlation is the square of the following sum

$$\hat{r}_{x_1,y}^2 = \left(\frac{\frac{1}{n}\sum_{i=1}^{n} (x_i-\bar{x})(y_i-\bar{y}) }{s_{x_1}s_z} \right)^2$$

It is not very rigorous but when $$n$$ is very large then we might approximate the expectation of this product as a product of expectations and also the sample deviation and sample mean can be approximated by the population deviation and population mean, $$s_{x_1} \approx \sigma_{x_1}$$, $$s_{y}^2 \approx \sigma_{y}^2$$, $$\bar{x} = \mu_x$$ and $$\bar{y} = \mu_y$$. For large sample $$n$$ this seems reasonable and we may expect any bias from this approximation to reduce to zero. Let's also assume without loss of generality that $$\mu_x=\mu_y=0$$

So

$$E[\hat{r}_{x_1,y}^2] \approx \left(\frac{\frac{1}{n}\sum_{i=1}^{n} E[x_iy_i] }{\sigma_{x_1}\sigma_y} \right)^2 = \frac{E[x^2]^2}{n\sigma_x^4} = \frac{1}{n}$$

In a similar fashion we should be able to obtain the expectation value for the Spearman correlation, but now we will work with a transformed variable. We are not looking for the correlation of the variables $$x$$ and $$y$$ but instead of the variables $$rank(x)$$ and $$rank(y)$$. For large $$n$$ these ranks will get closer to the quantiles of the distribution and we use these to approximate the limit for large $$n$$.

In the code below is a computation that uses the distribution of the ranks and a computational integration. For three variables the integration results in

### Spearman
[1] 0.3127959
### Pearson
[1] 0.3333071


And your expectation for the sum of the case with three variables is expected to approach $$0.9384$$.

We could say intuitively that the difference between the two correlations is in this transformation. With the Pearson correlation we compute correlation of $$x$$ and $$y$$ with the Spearman correlation we compute the correlation of transformed versions of $$x$$ and $$y$$.

So even if you would be able to increase your sample to a very large amount, then the expectation value for the Spearman's correlation will be lower in the case of a normal distribution of $$X$$.

set.seed(1)
n = 3

x_rank = function(x) {
p = pnorm(x,0,1)-0.5
}

y_rank = function(x) {
p = pnorm(x,0,sqrt(n))-0.5
}

Sigma = matrix(c(1,1,1,n),2)
m = solve(Sigma)/2
const = det(2*pi*Sigma)^(-0.5)

mvdnorm = function(x,y, Sigma) {
const*exp(-m[1,1]*x^2 - m[2,2]*y^2 - 2*m[1,2]*x*y)
}

dx = 0.1
range = seq(-8,8,dx)

xy = 0
rxry = 0
for (x in range) {
for (y in range) {
rxry = rxry + mvdnorm(x,y,Sigma) *x_rank(x)*y_rank(y)*dx^2
xy = xy + mvdnorm(x,y,Sigma)*x*y*dx^2
}
}

var_rx = sum(dnorm(range,0,1)*x_rank(range)^2)*dx
var_ry = sum(dnorm(range,0,sqrt(n))*y_rank(range)^2)*dx
var_x = sum(dnorm(range,0,1)*range^2)*dx
var_y = sum(dnorm(range,0,sqrt(n))*range^2)*dx

rxry^2/var_rx/var_ry
xy^2/var_x/var_y

rxry^2/var_rx/var_ry*n

v = MASS::mvrnorm(10^5,c(0,0),Sigma)

plot(v[,1],v[,2], pch = 19, col = rgb(0,0,0,0.03), cex = 0.5, xlab = "X", ylab = "Y", main = "r^2 = 0.3333")
plot(x_rank(v[,1]),y_rank(v[,2]), pch = 19, col = rgb(0,0,0,0.03), cex = 0.5, xlab = "quantile X", ylab = "quantile Y", main = "rho^2 = 0.3128")

cor(x_rank(v[,1]),y_rank(v[,2]))^2
cor((v[,1]),(v[,2]))^2

• Thanks for taking the effort of answering this mess of a question. Couple notes. The linked post seems to be about Pearson correlation of CDFs, not ranks like Spearman. The expression $E[\hat{r}^2] = \frac{n}{n-1} r^2$ is known as the Bessel's correction and your assumption is correct. Commented Jul 6, 2022 at 12:44
• What really bugs me about this is that the above inequality between Spearman's and Pearson's seem to persist for some other distributions like uniform and exponential. I'm wondering whether there's a more general phenomenon going on between the two. Commented Jul 6, 2022 at 12:49
• It doesn't seem that there can be a nice formula for the difference (factor?) between the two correlations especially one that works across all (can this be proven?) possible distributions. So I'm curious whether there's a way for the Spearman's to be brought in line with the behavior of Pearson's numerically. E.g. can we simply scale the whatever value of $\sum_{i=1}^{n}\rho_{X_i,Y}^2$ to 1.0 and apply the same factor to all the $\rho_{X_i,Y}^2$ and then claim the resulting $C\cdot \rho_{X_i,Y}^2$ to be the respective shares of the total variance? Would this be theoretically sound? Commented Jul 6, 2022 at 12:57
• @ayorgo I don't think that there is a general formula. Take for instance a sum of three heavy tailed distributions like a t-distribution with $\nu=3$ set.seed(1); n = 10^6; d = 3; X = matrix(rt(3*10^6,3),n); Y = rowSums(X); X1 = X[,1]; cor(X1,Y)^2*3; cor(X1,Y, method ="spearman")^2*3 now the squared Spearman's rank correlation is around 0.86, a bit lower. Commented Jul 6, 2022 at 13:25
• @ayorgo Responses to your first comment (1) In the limit of a large sample size, the Spearman's rank correlation will approach the Pearson correlation of the quantiles (the CDFs). (2) Indeed it looks like some sort of Bessel's correction. But, I would not consider it the same. Bessel's correction is for the estimate of the variance (and could be generalized to covariance) and it is not a correction for the correlation. Commented Jul 6, 2022 at 15:02