18
$\begingroup$

I am wondering how to derive the formula for the standard error of Pearson's correlation coefficient which is given in Zar for example as

$$ \newcommand{\cov}{{\rm Cov}} \newcommand{\var}{{\rm Var}} \newcommand{\sd}{{\rm SD}} SE_r =\sqrt{\frac{1-r^2}{n-2}}$$

I tried to get it from estimating the variance of r when

$$r =\frac{\cov(x,y)}{\sd(x)\sd(y)}$$

and $V(X) = E(X^2) - E(X)^2$ so we get $Var(r) = E\bigg(\frac{\cov(x,y)^2}{\var(x)\var(y)}\bigg) - r^2$. But from here I don't know how to continue since $E\bigg(\frac{\cov(x,y)^2}{\var(x)\var(y)}\bigg)$ would have to be $\frac{1-(n-3)r^2}{n-2}$ to get finally to

$$\var(r) =\frac{1-r^2}{n-2}$$

Any suggestions or references where I could look this up?

$\endgroup$
10
  • $\begingroup$ Maybe that helps en.youscribe.com/catalogue/tous/knowledge/formal-sciences/… $\endgroup$
    – Robert
    Jul 30, 2016 at 1:41
  • $\begingroup$ Thanks for the link but it just lists lots of formulas for SE of particular estimates but doesn't give any information of how these formulas are obtained. The fact that this cannot be found by using google somehow tells me that people just use it without knowing where it comes from which i find a bit scary... $\endgroup$
    – R. Carlos
    Jul 30, 2016 at 12:57
  • 1
    $\begingroup$ That's a really good question. Contrary to many textbooks,, that is not the standard error of r. For starters, the standard error of a statistic cannot be defined in terms of a statistic. It should be defined in terms of parameters. I believe you use rho in the formula and divide by sqrt(N-1) rather than 2 but I don't have a reference handy. $\endgroup$
    – David Lane
    Mar 16, 2017 at 2:11
  • 2
    $\begingroup$ One more reference: tandfonline.com/doi/abs/10.1080/… $\endgroup$
    – David Lane
    Mar 16, 2017 at 2:53
  • 4
    $\begingroup$ From all the references above and below, it seems that the formula you want to recover is wrong, anyway... $\endgroup$
    – Pascal
    Mar 13, 2018 at 12:03

3 Answers 3

8
$\begingroup$

There are two equations here for computing the statistical significance of the correlation coefficient.  The first is the variance of the true correlation coefficient $\rho$ of two bivariate normal random variables:\begin{equation}\text{var}\left(r\right)=\frac{\left(1-\rho^2\right)^2}{n},\end{equation} and the second is a t-statistic associated with the hypothesis that in the linear regression of $Y$ on $X$, the main effect of $X$ is zero:\begin{equation}t=r\sqrt{\frac{n-2}{1-r^2}}.\end{equation}Whence the standard error of $r$ mentioned by the OP: $\text{se}\left(r\right)=\sqrt{\frac{1-r^2}{n-2}}$.

Both of these expression can derived from the principle of maximum-likelihood.  That is, if we assume a parameter $\theta$ should be distributed normally,\begin{equation}\mathcal{L}\left(\theta\right)\sim\exp{\left(-\frac{\theta^2}{2\sigma^2_{\theta}}\right)},\end{equation}then the standard error of the parameter can be estimated from the curvature of the log-likelihood $\ell=\log{\mathcal{L}}$, function via\begin{equation}\sigma^2_{\theta}=\frac{-1}{\frac{\partial^2\ell}{\partial\theta^2}\bigr|_{\theta=\hat{\theta}}},\end{equation}where $\hat{\theta}$ is the maximum-likelihood estimate of $\theta$, got from the condition\begin{equation}\frac{\partial\ell}{\partial\theta}\bigr|_{\theta=\hat{\theta}}=0.\end{equation}

Now, Pearson derived the first expression \begin{equation}\text{var}\left(r\right)=\frac{\left(1-\rho^2\right)^2}{n\left(1+\rho^2\right)}\end{equation} in VII. Mathematical contributions to the theory of evolution.-III. Regression, heredity, and panmixia and https://royalsocietypublishing.org/doi/10.1098/rspl.1897.0091 by expanding the joint distribution of $n$ pairs of bivariate normal variables about the true value of $\rho$.  We can summarize his method here.  If we let $f$ be the bivariate normal density of two zero-mean random variables, i.e.,\begin{equation}f\left(X,Y\right)=\frac{1}{2\pi\sqrt{1-\rho^2}\sigma_X\sigma_Y}\exp{\left(-\frac{X^2}{2\sigma_X^2\left(1-\rho^2\right)}+-\frac{\rho XY}{2\sigma_X\sigma_Y\left(1-\rho^2\right)}-\frac{Y^2}{2\sigma_Y^2\left(1-\rho^2\right)}\right)},\end{equation}then we can get the variance $\sigma^2_{\rho}$ of the correlation coefficient by evaluating $\frac{-\partial^2\log{f}}{\partial \rho^2}\bigr|_{\rho=\hat{\rho}}$.  The first derivative of $\log{f}$ is\begin{align}\frac{\partial\log{f}}{\partial \rho}&=\frac{\rho}{1-\rho^2}+\frac{2\rho}{1-\rho^2}\cdot\frac{-X^2}{2\left(1-\rho^2\right)\sigma_X^2}+\left(\frac{1}{1-\rho^2}+\frac{2\rho^2}{\left(1-\rho^2\right)^2}\right)\frac{XY}{\sigma_X\sigma_Y}+\frac{2\rho}{1-\rho^2}\cdot\frac{-X^2}{2\left(1-\rho^2\right)\sigma_Y^2}\nonumber\\&=\frac{\rho}{1-\rho^2}+\left(\frac{2\rho}{1-\rho^2}\right)\left(\log{f}-\log{\frac{1}{2\pi\sigma_X\sigma_Y\sqrt{1-\rho^2}}}\right)+\frac{1}{1-\rho^2}\frac{XY}{\sigma_X\sigma_Y},\end{align}where at the maximum-likelihood solution $\hat{\rho}=\frac{\mathbb{E}\left(XY\right)}{\sigma_X\sigma_Y}$ the middle term becomes\begin{equation}\left[\log{f}-\log{\frac{1}{2\pi\sigma_X\sigma_Y\sqrt{1-\rho^2}}}\right]\bigr|_{\rho=\hat{\rho}}=-\frac{\mathbb{E}\left(X^2\right)}{\sigma_X^2\left(1-\rho^2\right)}+\frac{\rho\mathbb{E}\left(XY\right)}{\sigma_X\sigma_Y\left(1-\rho^2\right)}-\frac{\mathbb{E}\left(Y^2\right)}{\sigma_X^2\left(1-\rho^2\right)}=-1.\end{equation}Whence upon taking the second derivative and evaluating at $\rho=\hat{\rho}$, we get\begin{align}\frac{\partial^2\log{f}}{\partial\rho^2}\bigr|_{\rho=\hat{\rho^2}}&=\left(\frac{1}{1-\rho^2}+\frac{2\rho^2}{\left(1-\rho^2\right)^2}\right)+\frac{2\rho}{\left(1-\rho^2\right)^2}\frac{\mathbb{E}\left(XY\right)}{\sigma_X\sigma_Y}+\left(\frac{2}{1-\rho^2}+{4\rho^2}{\left(1-\rho^2\right)}\right)\left[\log{f}-\log{\frac{1}{2\pi\sigma_X\sigma_Y\sqrt{1-\rho^2}}}\right]\bigr|_{\rho=\hat{\rho}}+\frac{2\rho}{1-\rho^2}\frac{\partial}{\partial\rho}\left[\log{f}-\log{\frac{1}{2\pi\sigma_X\sigma_Y\sqrt{1-\rho^2}}}\right]\bigr|_{\rho=\hat{\rho}}\nonumber\\&=\frac{1+\rho^2}{\left(1-\rho^2\right)^2}+\frac{2\rho^2}{\left(1-\rho^2\right)^2}-2\left(\frac{1+\rho^2}{\left(1-\rho^2\right)^2}\right)-\frac{2\rho^2}{\left(1-\rho\right)^2},\end{align}so that\begin{equation}\sigma^2_{\rho}=\frac{-1}{\frac{\partial^2\log{f}}{\partial\rho^2}\bigr|_{\rho=\hat{\rho}}}=\frac{\left(1-\rho^2\right)^2}{1+\rho^2}.\end{equation}Then by the central limit theorem, the sampling variance of $\rho$ is $\frac{\left(1-\rho^2\right)^2}{n\left(1+\rho^2\right)}$.

For the second form of the statistic, let's drop the assumption of bivariate normality and consider the regression of $Y$ on $X$ with normally-distributed error $\varepsilon$: if\begin{equation}Y=\alpha+\beta X+\varepsilon\end{equation}and\begin{equation}\sigma^2_Y=\beta^2\sigma^2_X+\sigma^2,\end{equation} then according to the relationship $\beta=\rho\frac{\sigma_Y}{\sigma_X}$, it must be the case that error variance is $\sigma^2=\left(1-\rho^2\right)\sigma^2_Y$.  Then the distribution of $\varepsilon$ is:\begin{equation}\mathcal{L}\left(\varepsilon\right)\sim\Pi_i\exp{\left(-\frac{\left(Y_i-\alpha-\beta X_i\right)^2}{2\sigma^2}\right)},\end{equation}so that the log-likelihood is\begin{equation}\ell=-\sum_i\frac{\left(Y_i-\alpha-\beta X_i\right)^2}{2\left(1-\rho^2\right)\sigma_Y^2}.\end{equation}The first two derivatives are\begin{equation}\frac{\partial\ell}{\partial \beta}=\sum_i\frac{\left(Y_i-\overline{Y}-\beta\left(X_i-\overline{X}\right)\right)\left(X_i-\overline{X}\right)}{\left(1-\rho^2\right)\sigma_Y^2}\end{equation}and\begin{equation}\frac{\partial^2\ell}{\partial\beta^2}=-\sum_i\frac{\left(X_i-\overline{X}\right)^2}{\left(1-\rho^2\right)\sigma_Y^2}.\end{equation}Now, making the substitution $\beta=\rho\frac{\sigma_Y}{\sigma_X}$ and evaluating at $\rho=\hat{\rho}=r$ gives\begin{equation}\frac{-1}{\frac{\partial^2\ell}{\partial\beta^2}\bigr|_{\beta=\hat{\beta}}}=\frac{-1}{\frac{\partial^2\ell}{\partial\rho^2}\bigr|_{\rho=\hat{\rho}}}\frac{\sigma^2_Y}{\sigma^2_X}=\frac{\sigma^2_Y\left(1-r^2\right)}{\sigma^2_X\left(n-2\right)},\end{equation}whence the sampling variance of the measured correlation coefficient $\hat{\rho}=r$ is\begin{equation}\sigma^2_r=\frac{1-r^2}{n-2},\end{equation} in which we lose two degrees of freedom from the estimation of two parameters $\alpha$ and $\beta$.  Finally, we can form a t-statistic to test the hypothesis that $r=0$ using \begin{equation}t=\frac{r}{\sigma_r}=r\sqrt{\frac{n-2}{1-r^2}}.\end{equation}For an alternate derivation, see The Analysis of Physical Measurements, pp. 193-199, by Pugh and Winslow cited in A brief note on the standard error of the Pearson correlation.

A comparison of the two formulas shows\begin{equation}\frac{\sigma^2_{\rho}}{n}=\frac{\left(1-\rho^2\right)^2}{n\left(1+\rho^2\right)}<\frac{1-r^2}{n-2}=\sigma^2_r.\end{equation}In other words, there is less variance of the true parameter $\rho$ than in the value estimated from linear regression.  It should also be pointed out that Pearson's formula is only true for bivariate normal variables, while the standard error of $r$ is valid for any linear regression.  However, we see that the test of whether $r=0$ is equivalent to the test of whether $\beta=0$ and does not really tell us anything new.

$\endgroup$
2
  • $\begingroup$ I didn't follow everything but the derivation looks well done, at least it is in accordance with my simulations and on the fact that the OP has to square the numerator of his formula. But why did you cross out some sentences, do you think they are wrong or is it just a style that you apply? $\endgroup$
    – Jean Paul
    Feb 18, 2023 at 14:00
  • $\begingroup$ I had a change of heart about my first derivation and decided it must be this way instead. However, I'd still like to see how we could derive the variance directly from E(X^2Y^2) - E(X)^2E(Y)^2, because I was getting nonsensical results (that had to be crossed out) from that approach. $\endgroup$ May 2, 2023 at 21:22
7
$\begingroup$

After looking for a long time for an answer to this same question, I found a couple interesting links:

$\bullet$ The Standard Deviation of the Correlation Coefficient, where we can only see the first page but that's where the derivation is. The "standard deviation by dr Sheppard" is given by something called the Asymptotic distribution of moments, of which you can see a bit in the following source.

$\bullet$ A History of Parametric Statistical Inference from Bernoulli to Fisher, 1713-1935.

The reason for the "n-2" instead of "n" in the root, is that your formula assumes a t-distribution with n-2 degrees of freedom, while the one in the links assumes a normal distribution.

$\endgroup$
2
  • 1
    $\begingroup$ It looks like the equal sign in the formula should be written as an approximately sign then cause it reflect asymptotic behaviour rather than a closed form solution. Further, the SE of the correlation coefficient is also then heavily dependent on the sample size as far as I read it from Pearson original papers. Conclusion is that the formula for the SE of r cannot be simply obtained as it is the case for the SE of the mean or the variance. $\endgroup$
    – R. Carlos
    Apr 24, 2017 at 2:22
  • $\begingroup$ Although it makes sense to have n-2 in the denominator since we use r as an estimate of $\rho$ there still must be somewhere a mathematical derivation of this which I would like to see. With the variance estimated from a sample we also get n-1 in the denominator which is mathematically derived as well and not just put because we reduce the degrees of freedom by one when using the samples expected value. $\endgroup$
    – R. Carlos
    Jul 10, 2022 at 8:29
-2
$\begingroup$

I do not have the answer, but for me there is an error in the formula of the question.

It is: $$SE_r =\frac{1-r^2}{\sqrt{n-2}}$$ and $$\newcommand{\Var}{\operatorname{Var}} \Var(r)=\frac{(1-r^2)^2}{n-2}$$

I will try to check this by simulation:

library(MASS)

N = 100000
r = 0.8
n = 100

Sigma = matrix(c(1, r, r, 1), nrow=2)
r_obs =  replicate(N, cor(mvrnorm(n, c(0,0), Sigma))[2,1])

> mean(r_obs)
[1] 0.7984783
> sd(r_obs)
[1] 0.03690896

So the standard error for r=0.8 with n=100 is approximately 0.037.

If I use the formula of the question I get:

> sqrt((1-r^2)/(n-2))
[1] 0.06060915

And with the formula I gave:

> (1-r^2)/sqrt((n-2))
[1] 0.03636549

The second formula seems to be much closer to the true value than the first.

Edit:

To try to explain why there are two different formulas for the standard error which are circulating, I found that it depends on how you compute it.

In my first simulation, I used Pearson formula to compute the correlation, but one can also use the least square regression coefficient. I can confirm that the latter method has the standard error proposed in the question:

r_obs <- replicate(N, {
    M<-mvrnorm(n, c(0,0), Sigma)
    c(pearson=cor(M, method="pearson")[2,1],
      regression=lm(M[,2]~M[,1])$coef[[2]])
    })

> apply(r_obs, 1, mean)
   pearson regression 
 0.7981580  0.7998433 
> apply(r_obs, 1, sd)
   pearson regression 
0.03707184 0.06094964 

These are two estimators of the correlation which do not have the same variance.

Edit 2:

That try to reconcile the two formulas did not work, because I forgot to normalize the regression coefficient. The formula to compute the correlation from the regression coefficient is: $$r=\frac{SD(x)}{SD(y)}b$$ with $$E(y|x)=a+b\cdot x$$

By redoing the correct computation, I obtain in fact the same result:

r_obs <- replicate(N, {
    M<-mvrnorm(n, c(0,0), Sigma)
    c(pearson=cor(M, method="pearson")[2,1],
      regression=lm(M[,2]~M[,1])$coef[[2]]*sd(M[,1])/sd(M[,2]))
    })
> apply(r_obs, 1, sd)
   pearson regression 
0.03676248 0.03676248 
> apply(r_obs, 1, mean)
   pearson regression 
 0.7992615  0.7992615 

Which is reassuring in some way. So I maintain that the standard error formula in the question is incorrect, but maybe I explained where does the error come from.

$\endgroup$
6
  • 2
    $\begingroup$ There is no error in the original formula: see en.wikipedia.org/wiki/Pearson_correlation_coefficient. You will discover your formula doesn't seem to work, either, if you set n <- 3 instead of 100 and r <- 0 instead of 0.8. There are several problems, one of which is that the formula of the question is a reasonable standard error only when $r$ is close to zero. Another is that the sampling distribution can be highly skewed, rendering any SE less than useful: that's why better analyses re-express $r$ with the Fisher transformation. $\endgroup$
    – whuber
    Nov 6, 2018 at 16:02
  • 2
    $\begingroup$ I think that n=3 is a too special case to judge. I don't know, maybe n-2 is not the good divisor, maybe it is n-3. Beyond that, the two formulae give the same result for r=0. $\endgroup$
    – Jean Paul
    Nov 6, 2018 at 17:05
  • 1
    $\begingroup$ Please explain the downvotes. In addition, I would like to know where I can find the formula of the standard error of the raw correlation coefficient (not the Fisher transform) in the Wikipedia page. Logically it should be in this section: en.wikipedia.org/wiki/Pearson_correlation_coefficient#Inference, but I am unable to see it. $\endgroup$
    – Jean Paul
    Oct 10, 2019 at 15:36
  • $\begingroup$ First of all: 'Checking' by simulation is not a mathematical proof. Second, just because it is not on wikipedia it does not mean it is not existent. It just means no one has yet copy/paste information to it. $\endgroup$
    – R. Carlos
    Jul 10, 2022 at 8:19
  • $\begingroup$ Yes checking by simulation is not a mathematical proof, but if simulations do not match a formula but match another there are reasons to think that the other formula is better. For Wikipedia it is just the other guy who suggested me to check the formula there, I just answered that I was unable to find it, I did not deduce anything more. But now the formula of the standard error has been added, and guess what, that's not the formula of the question but the one of my answer which is written there. $\endgroup$
    – Jean Paul
    Jul 22, 2022 at 19:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.