Inside the Subspace Where Spurious Correlations Are Born

between cholesterol and blood pressure is measured in 10 patients. The observed correlation is 0.62.

Is this enough evidence that a relationship exists?

Studies that measure 20,000 gene expression levels in 10 mice routinely apply multiple-testing corrections and stringent correlation thresholds. Researchers know that the more correlations they examine, the more likely they are to observe large values purely by chance, even when the genes are completely unrelated.

With small studies, researchers often assume they are safer from spurious correlations. Are they?

The distribution of a sample correlation is the same regardless of the total number of variables measured. It depends primarily on the number of subjects. High-dimensional datasets do not create spurious correlations—they simply make them easier to encounter.

This article answers the question: what correlation values should we expect when the variables are independent and the true population correlation is zero? It also explains why such correlations occur.

Using the geometry of Pearson’s correlation coefficient, the article visualizes the effects of centering and normalization. It then discusses the role of rotational invariance. Together, these ideas build intuition for the exact distribution of Pearson’s correlation coefficient and its asymptotic normal approximation.

Finally, the article offers practical guidance on how to use this intuition when reading and conducting research.

The code to reproduce all figures and experiments is available in the accompanying notebook.

The Gaussian Experiment

Imagine measuring d variables on n individuals.

To study correlations that arise purely by chance, consider a simulated dataset whose columns are independent Gaussian random vectors of dimension n. The resulting feature smatrix, where rows correspond to subjects and columns to variables, is

with columns satisfying

What is the sample correlation between two vectors?

Even though the population correlation between any pair of variables is exactly zero, the sample correlations we observe are almost always nonzero. This raises a question: when the columns are truly independent, what values should we expect for the sample correlation?

The sample (Pearson) correlation between columns k and l is

where

is the sample mean of the k-th column. The vector of all column means is

Because the columns are random vectors, the sample correlation $r\,(X_{\cdot k},\ X_{\cdot l})$ is itself a random variable taking values in $\left[-1, 1\right]$.

The intuition behind its derivation becomes clearer by understanding how Pearson’s correlation coefficient transforms the data.

Step 1: centering moves vectors onto a hyperplane

The numerator is effectively centering the vectors.

Centering subtracts the sample mean from each column

where $\mathbf{1}_n$ is the column vector of ones

After centering, the coordinates of every column sum to 0

This means every centered column now lies in the (n-1)-dimensional subspace

perpendicular to the vector of ones​. The vectors remain n-dimensional, but after centering they are restricted to an (n-1)-dimensional hyperplane. They no longer vary in the direction of the all-ones vector.

Sample size n = 3, centering constrains all vectors to a plane

For n = 3, the hyperplane is an ordinary plane in three-dimensional space.

Rotating the plane to align with the xy-plane makes it visible: there is no variability in the z-direction, and the vectors are restricted to a two-dimensional subspace.

Sample size n = 4: centering constrains all vectors to a three-dimensional hyperplane

After choosing a coordinate system orthogonal to the all-ones vector, the centered vectors can be visualized in three dimensions. They form a cloud centered at the origin.

Step 2: normalizing moves vectors onto the unit sphere

Notice that Pearson’s correlation is exactly the cosine of the angle between two centered vectors

The distribution of the sample correlation between two independent random vectors is therefore the same as the distribution of the cosine of the angle between their centered versions.

When the centered columns are normalized

the correlation matrix simplifies to

Centering moved the vectors onto an (n−1)-dimensional hyperplane, removing information about location. Normalizing removes differences in length

removing information about scale. The vectors are now constrained to the unit sphere embedded in the (n-1)-dimensional hyperplane, an (n−2)-dimensional manifold

n = 3: normalizing moves vectors onto the unit circle

For n = 3, centering moved the vectors onto a plane. Normalizing fixes their length at one, restricting them to the one-dimensional sphere within that plane—the unit circle.

n = 4: normalizing moves vectors onto the unit sphere

For n = 4, centering moved the vectors into a cloud centered at the origin within a three-dimensional subspace. Normalizing fixes their length at one, placing them on the unit sphere.

Notice that Gaussianity has not been used so far. Every random vector becomes a point on the same sphere after centering and normalization.

Step 3: Gaussianity allows the derivation of the exact distribution

The standard multivariate normal distribution is rotationally invariant. This means that the probability density of a random vector depends only on its distance from the origin, not on its direction.

Conditional on a radius, directions are uniformly distributed on a sphere. Therefore, after normalization to unit length, the vectors become uniformly distributed on the unit sphere.

Assuming the vectors are independent and Gaussian, the exact sampling distribution of the correlation C can be derived by fixing one direction and measuring the surface area of all points on the sphere that produce a given correlation value. Rotational invariance guarantees that a result obtained for one fixed direction applies equally to every direction. After obtaining the result and normalizing it to 1, the exact sampling distribution of the correlation is

where B is the Beta function. The numerator reflects the geometry of the sphere, while the denominator is the normalizing constant.

Since the variance decreases as $\frac{1}{n−1}$, the typical magnitude of the correlation decreases as $(\sqrt{n})^{-1}$. So under the independent Gaussian null the typical correlations are of order

Equivalently, the squared correlation under the independent Gaussian model follows a Beta distribution

which makes exact quantiles and significance thresholds easy to compute.

n = 3: correlation distribution is U-shaped

This distribution is a shifted and rescaled version of the classical arcsine distribution.

For a study with only 3 subjects, the two-sided 5% rejection threshold for the null hypothesis of zero population correlation is ∣r∣>0.997. Smaller correlations are not unusual enough to reject the hypothesis of no relationship.

Notice that with only 3 subjects, it is more likely to observe a large correlation than a correlation around 0, even though the variables are unrelated.

n = 4: correlation follows a uniform distribution

Because the numerator in the exact sampling distribution formula simplifies to 1 when n = 4

For a study with 4 subjects, the rejection threshold remains very high: $∣r∣>0.95$. Smaller correlations provide insufficient evidence to reject the hypothesis of no relationship.

For this sample size, all correlations are equally likely to be observed.

n = 5: correlation distribution is semicircle-shaped

It is a Wigner semicircle distribution.

For n = 6 and higher the distribution becomes bell-shaped. As the sample size increases, it approaches the asymptotic normal distribution.

Step 4: when n goes to infinity

As n grows, the dimension of the sphere increases and almost all pairs of random directions become nearly orthogonal. As a result, the distribution of the correlation concentrates around zero. More precisely,

Therefore, for sufficiently large n

What about non-normal sample distributions?

Surprisingly, the asymptotic result also holds for non-normal data, provided the underlying distribution has finite variance.

More surprisingly, the exact sampling distribution often remains an excellent approximation even when the data is non-normal and the sample size is small.

Exponential distribution

The derivation of the exact sampling distribution relies on directions being uniformly distributed on the unit sphere and therefore requires rotational invariance. The exponential distribution does not possess this property.

This sphere is obtained by centering and normalizing independent exponential random vectors. Unlike the Gaussian case, the points are not uniformly distributed on the sphere and some directions are preferred. The correlation distribution no longer follows the exact law derived earlier.

For small n, the empirical sampling distribution of the correlation is very different from the exact sampling distribution. For this simulated exponential example, empirical quantiles would provide a more appropriate reference for hypothesis testing than either the exact Gaussian distribution or the asymptotic normal approximation.

As n increases, the empirical sampling distribution of the correlation starts resembling the asymptotic normal distribution.

For n = 10, the empirical sampling distribution is still neither close to the asymptotic normal distribution nor to the exact sampling distribution of the correlation. The bell shape is visibly skewed. This has practical consequences: the rejection region based on the exact distribution is no longer correct, even though the observed Type I error rate may appear approximately correct or even slightly conservative.

For a sample size of 100, the empirical distribution is fairly close to the asymptotic normal distribution and nearly symmetric.

Although some residual skewness remains, the approximation can be considered good enough.

Symmetric distributions

A very interesting thing happens for simulated symmetric distributions. As long as the variance is finite, the empirical distribution is very close to the exact distribution derived under the Gaussian independence assumption.

This is true for sample sizes as small as n = 3 for the Laplace and uniform distributions.

A transformed U-shaped $\mathrm{Beta}(0.2, 0.2)$ distribution, symmetric around 0, behaves differently. Its empirical sampling distribution of Pearson’s correlation differs from the exact one.

No wonder. The vectors are not uniformly distributed on the unit sphere.

Surprisingly, it still converges to the exact sampling distribution of the correlation.

For n as small as 8, the empirical distribution looks very similar to the exact one, even though the latter was derived under the assumption of rotational invariance.

In these simulations, skewness appears to be the dominant source of deviation from the exact sampling distribution derived under the normality assumption. The Cauchy distribution, which has infinite variance, does not show the same convergence. But the simulated symmetric non-normal distributions with finite variance do move toward the exact sampling distribution of the correlation.

Notice that the convergence happens long before the asymptotic normal approximation becomes accurate.

What does it all mean for a practitioner?

Small studies are susceptible to spurious correlations

A common misunderstanding is that spurious correlations only arise in large (high-dimensional) datasets. But in a small study with the same number of subjects, each sample correlation comes from exactly the same sampling distribution.

It’s like rolling three fair dice. The probability of getting a total greater than 16 is close to 2%. This extreme result can occur on the first try. But as the number of rolls increases, the chance of observing at least one such outcome also increases. If you roll the dice 100 times, the probability of observing at least one sum greater than 16 is about 84%.

This is why high-dimensional datasets contain extreme correlations more often, even when the variables are unrelated. But a small study with the same number of subjects draws its sample correlations from exactly the same sampling distribution. The same kind of extreme correlation can therefore occur purely by chance.

Correlations should never be interpreted in isolation from the sample size

In many applied fields, correlations of 0.4 or greater are often described as strong.

Consider a study with only 3 subjects. Even when the variables have no true relationship, correlations are more likely to be close to ±1 than to 0. The probability that $|r| > 0.4$ is close to 74% when the variables are independent.

With 10 observations, there is still a 25% chance of getting a value $|r| > 0.4$.

This value cannot be considered evidence of a relationship in small-sample studies.

In a study of 10 mice measuring the expression levels of about 20,000 genes against one unrelated outcome, we expect over 1300 pairs with $|r| > 0.6$ even when no biological relationship exists. This is why multiple-testing corrections are routinely used in high-dimensional datasets.

Does it mean small-sample studies shouldn’t be done or published? Not at all. The experiment can be replicated, and eventually, when several laboratories repeat it, a meta-analysis can place each result in context.

For n = 100 subjects, the picture is different. The probability that $∣r∣>0.4$ is only 0.00359% (about 36 per million).

For this sample size, absolute correlations of 0.2 and larger are considered statistically significant. The upper 2.5% critical value for various sample sizes is shown in the table below.

Take note, however, that statistical significance does not guarantee that a relationship is real. The null distribution has support on the entire interval $[-1,1]$, so every possible sample correlation can arise even when the true population correlation is zero. Correlations in the rejection region are simply improbable under the null, not impossible. They are the source of type I errors in hypothesis testing and can still lead to false discoveries.

The test based on the exact distribution is surprisingly robust

Many software implementations, including scipy.stats.pearsonr, use the exact null distribution of Pearson’s correlation coefficient. Because it is derived for independent Gaussian variables, researchers might think they should use other methods for non-normal data. It appears that the exact sampling distribution of Pearson’s correlation coefficient provides a surprisingly good approximation, long before the sampling distribution approaches its large-sample normal limit—sometimes even for n = 3 or n = 4.

This empirical robustness is remarkably strong for simulated symmetric distributions with finite variance and holds in cases far from asymptotic conditions. Theoretical explanations or references from readers are welcome.

When the underlying distribution is known to be skewed, one can estimate empirical null quantiles by simulating the distribution of the correlation under the null hypothesis of independence. For the simulated exponential distribution, the estimated quantiles are:

These are Monte Carlo estimates from simulations (one run); more stable values can be obtained with more simulation runs.

For larger sample sizes, the midpoint moves toward 0, and the empirical quantiles move closer to both the exact Gaussian null distribution and the asymptotic normal approximation.

The one takeaway

If there is only one sentence I could recommend you remember, it is this: investigate the sample size of any study that boasts a remarkable result.

It could be a genuine discovery, but it could also be two vectors that happened to land unusually close together on a sphere by chance.

Feedback welcome

If you’d like to replicate the simulations, please visit the GitHub repository. All figures were created by me, and the simulation-based figures can be reproduced from the code there.

If you enjoyed the article or have a suggestion, don’t hesitate to contact me on LinkedIn. All thoughtful feedback is greatly appreciated.