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

Question

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]<E \left[\sum_{i=1}^{n}r_{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")
    )
)
fig <- fig %>% add_surface(
    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}')
)
fig <- fig %>% add_surface(
    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}')
)
fig <- fig %>% add_surface(
    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.

Answer 1

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