The equivalence of sample correlation and R statistic for simple linear regression

Question

It is often stated that the square of the sample correlation $r^2$ is equivalent to the $R^2$ coefficient of determination for simple linear regression. I have been unable to demonstrate this myself and would appreciate a full proof of this fact.

Answer 1

There seems to be some variation in notation: in a simple linear regression, I've usually seen the phrase "sample correlation coefficient" with symbol $r$ as a reference to the correlation between observed $x$ and $y$ values. This is the notation I have adopted for this answer. I have also seen the same phrase and symbol used to refer to the correlation between observed $y$ and fitted $\hat y$; in my answer I have referred to this as the "multiple correlation coefficient" and used the symbol $R$. This answer addresses why the coefficient of determination is both the square of $r$ and also the square of $R$, so it shouldn't matter which usage was intended.

The $r^2$ result follows in one line of algebra once some straightforward facts about correlation and the meaning of $R$ are established, so you may prefer to skip down to the boxed equation. I assume we don't have to prove basic properties of covariance and variance, in particular:

$$\text{Cov}(aX+b, Y) = a\text{Cov}(X,Y)$$ $$\text{Var}(aX+b) = a^2\text{Var}(X)$$

Note that the latter can be derived from the former, once we know that covariance is symmetric and that $\text{Var}(X)= \text{Cov}(X,X)$. From here we derive another basic fact, about correlation. For $a \neq 0$, and so long as $X$ and $Y$ have non-zero variances,

$$\begin{align} \text{Cor}(aX+b, Y) &= \frac{\text{Cov}(aX+b, Y)}{\sqrt{\text{Var}(aX+b) \text{Var} (Y)}} \\ &= \frac{a}{\sqrt{a^2}} \times \frac{\text{Cov}(X, Y)}{\sqrt{\text{Var}(X) \text{Var} (Y)}} \\ \text{Cor}(aX+b, Y) &= \text{sgn}(a) \, \text{Cor}(X,Y) \end{align} $$

Here $\text{sgn}(a)$ is the signum or sign function: its value is $\text{sgn}(a) = +1$ if $a>0$ and $\text{sgn}(a) = -1$ if $a<0$. It's also true that $\text{sgn}(a) = 0$ if $a=0$, but that case doesn't concern us: $aX+b$ would be a constant, so $\text{Var}(aX+b) = 0$ in the denominator and we can't calculate the correlation. Symmetry arguments let us generalise this result, for $a, \, c \neq 0$:

$$\text{Cor}(aX+b, \, cY+d) = \text{sgn}(a) \, \text{sgn}(c) \, \text{Cor}(X,Y)$$

We won't need this more general formula to answer the current question, but I include it to emphasise the geometry of the situation: it simply states that correlation is unchanged when either variable is scaled or translated, but reverses in sign when a variable is reflected.

We need one more fact: for a linear model including a constant term, the coefficient of determination $R^2$ is the square of the multiple correlation coefficient $R$, which is the correlation between the observed responses $Y$ and the model's fitted values $\hat Y$. This applies for both multiple and simple regressions, but let us restrict our attention to the simple linear model $\hat Y = \hat \beta_0 + \hat \beta_1 X$. The result follows from the observation that $\hat Y$ is a scaled, possibly reflected, and translated version of $X$:

$$\boxed{R = \text{Cor}(\hat Y, Y) = \text{Cor}(\hat \beta_0 + \hat \beta_1 X, \, Y) = \text{sgn}(\hat \beta_1) \, \text{Cor}(X, Y) = \text{sgn}(\hat \beta_1) \, r}$$

So $R = \pm r$ where the sign matches the sign of the estimated slope, which guarantees $R$ not to be negative. Clearly $R^2 = r^2$.

The preceding argument was made simpler by not having to consider sums of squares. To achieve this, I skipped over the details of the relationship between $R^2$, which we normally think of in terms of sums of squares, and $R$, for which we think about correlations of fitted and observed responses. The symbols make the relationship $R^2 = (R)^2$ seem tautological but this is not the case, and the relationship breaks down if there is no intercept term in the model! I'll give a brief sketch of a geometric argument about the relationship between $R$ and $R^2$ taken from a different question: the diagram is drawn in $n$-dimensional subject space, so each axis (not shown) represents a single unit of observation, and variables are shown as vectors. The columns of the design matrix $\mathbf{X}$ are the vector $\mathbf{1_n}$ (for the constant term) and the vector of observations of the explanatory variable, so the column space is a two-dimensional flat.

The fitted $\mathbf{\hat{Y}}$ is the orthogonal projection of the observed $\mathbf{Y}$ onto the column space of $\mathbf{X}$. This means the vector of residuals $\mathbf{e} = \mathbf{y} - \mathbf{\hat{y}}$ is perpendicular to the flat, and hence to $\mathbf{1_n}$. The dot product is $0 = \mathbf{1_n} \cdot \mathbf{e} = \sum_{i=1}^n e_i$. As the residuals sum to zero and $Y_i = \hat{Y_i} + e_i$, then $\sum_{i=1}^n Y_i = \sum_{i=1}^n \hat{Y_i}$ so that both fitted and observed responses have mean $\bar{Y}$. The dashed lines in the diagram, $\mathbf{Y} - \bar{Y}\mathbf{1_n}$ and $\mathbf{\hat{Y}} - \bar{Y}\mathbf{1_n}$, are therefore the centered vectors for the observed and fitted responses, and the cosine of the angle $\theta$ between them is their correlation $R$.

The triangle these vectors form with the vector of residuals is right-angled since $\mathbf{\hat{Y}} - \bar{Y}\mathbf{1_n}$ lies in the flat but $\mathbf{e}$ is orthogonal to it. Applying Pythagoras:

$$\|\mathbf{Y} - \bar{Y}\mathbf{1_n}\|^2 = \|\mathbf{Y} - \mathbf{\hat{Y}}\|^2 + \|\mathbf{\hat{Y}} - \bar{Y}\mathbf{1_n}\|^2 $$

This is just the decomposition of the sums of squares, $SS_{\text{total}} = SS_{\text{residual}} + SS_{\text{regression}}$. The conventional formula for the coefficient of determination is $1 - \frac{SS_{\text{residual}}}{SS_{\text{total}}}$ which in this triangle is $1 - \sin^2 \theta = \cos^2 \theta$ so is indeed the square of $R$. You may be more familiar with the formula $R^2 = \frac{SS_{\text{regression}}}{SS_{\text{total}}}$, which immediately gives $\cos^2 \theta$, but note that $1 - \frac{SS_{\text{residual}}}{SS_{\text{total}}}$ is more general, and will (as we've just seen) reduce to $\frac{SS_{\text{regression}}}{SS_{\text{total}}}$ if a constant term is included in the model.

Answer 2

The $R^2$ is defined as $$R^2=\frac{\hat{V}(\hat{y}_i)}{\hat{V}(y_i)} =\frac{1/(N-1)\sum_{i=1}^N(\hat{y}_i-\bar{y})^2}{1/(N-1)\sum_{i=1}^N(y_i-\bar{y})^2}=\frac{ESS}{TSS} $$ The squared sample correlation coefficient: $$r^2(y_i,\hat{y}_i)=\frac{\left(\sum_{i=1}^N(y_i-\bar{y})(\hat{y}_i-\bar{y})\right)^2}{\left(\sum_{i=1}^N(y_i-\bar{y})^2\right)\left(\sum_{i=1}^N(\hat y_i-\bar{y})^2\right)} $$ is equivalent, as it is easily verified using: $$\hat V(y_i)=\hat V(\hat y_i)+\hat V(e_i) $$ (see Verbeek, §2.4)

Answer 3

We consider here the problem of multiple linear regression, that is we are given $n\in\mathbb N$ data points (i.e. vectors in $\mathbb R^{p+1}$ for some $p\in\mathbb N\cup\{0\}$) $$\mathbf x_1 = \begin{pmatrix} x_{1,0}, x_{1,1},\dots, x_{1,p}\end{pmatrix},\ \mathbf x_2 = \begin{pmatrix} x_{2,0}, x_{2,1},\dots, x_{2,p}\end{pmatrix},\dots,\mathbf x_n = \begin{pmatrix} x_{n,0}, x_{n,1},\dots, x_{n,p}\end{pmatrix}.$$ We are also given the corresponding labels $$y_1,\dots,y_n\in\mathbb R.$$

Of course, in practice, $n$ denotes the number of data points we have available while $p$ denotes the number of features we use for prediction.

I will write $$\mathbf y = \begin{pmatrix}y_1\\ \vdots\\y_n\end{pmatrix}\in\mathbb R^n$$ as well as $$\mathbf X = \begin{pmatrix}\mathbf x_1\\ \vdots\\ \mathbf x_n\end{pmatrix}\in\mathbb R^{n\times (p+1)}.$$

Our goal is to predict the labels based on the data points using linear regression while minimizing the Gaussian squared loss [Footnote 1], that is, we are looking for a $\beta\in\mathbb R^{p+1}$ which minimizes

$$L(\mathbf X, \mathbf y, \beta)=\lVert\mathbf y-\mathbf X\beta\rVert^2=\sum_{k=1}^n \left(y_k - \langle\mathbf x_k, \beta\rangle\right)^2.$$

A direct computation of the gradient of $L$ gives

$$\nabla_\beta L(\mathbf X, \mathbf y, \beta) = 2(\mathbf X^\top\mathbf X\beta - \mathbf X^\top\mathbf y).$$

The Hessian of $L$ satisfies

$$\text{Hess}_\beta L(\mathbf X, \mathbf y, \beta) = \mathbf X^\top\mathbf X.$$

In particular, since $\mathbf X^\top\mathbf X$ is positive semi-definite, $\beta\mapsto L(\mathbf X,\mathbf y, \beta)$ is convex for all fixed $\mathbf X, \mathbf y$, and thus any solution $\beta$ of $\nabla_\beta L(\mathbf X,\mathbf y,\beta)=0$ is a global minimizer of $\beta\mapsto L(\mathbf X,\mathbf y, \beta)$.

$\nabla_\beta L(\mathbf X,\mathbf y,\beta)=0$ is equivalent to $$\mathbf X^\top\mathbf X\beta - \mathbf X^\top\mathbf y=0,$$ so that, assuming $\mathbf X^\top\mathbf X$ is invertible, as we will do from now on, (this is the case if the rank of $\mathbf X$ is equal to $\min(n,p)$), we have that $$\beta = (\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\mathbf y$$ is a minimizer of our quadratic loss.

Define now

$$\mathbf H = \mathbf X(\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\in\mathbb R^{n\times n}.$$

Then the predictions of our linear model are $$\hat{\mathbf y}\overset{\text{Def.}}=\mathbf H\mathbf y.$$

We define the "coefficient of determination" [Footnote 2] $R^2$ to be equal to $$ R^2 \overset{\text{Def.}}= \frac{\sum_{k=1}^n (\hat y_k-\bar y)^2}{\sum_{k=1}^n (y_k-\bar y)^2}, $$ where $\hat y_k$ denotes the $k$-th coordinate of $\hat{\mathbf y}$ and $\bar y\overset{\text{Def.}}=\frac{y_1+\dots y_n}{n}$. In fact, this Definition is common also if the predictions $\hat{\mathbf y}$ were obtained by methods different than linear regression. (Note that if $\sum_{k=1}^n (y_k-\bar y)^2=0$, then $R^2$ is undefined.)

---

The goal of this first section will be establish basic relationships between $\hat{\mathbf y}$ and $\mathbf y$. As a Corollary of these relationships we will obtain that if $R^2$ is well-defined, then $$R^2 = 1-\frac{\sum_{k=1}^n (y_k-\hat y_k)^2}{\sum_{k=1}^n (y_k-\bar y_k)^2}.$$

Lemma 1a. We have $\mathbf H^\top = \mathbf H$ and $\mathbf H^2 = \mathbf H$.

Proof. Straightforward calculation. $\square$

Lemma 1b. Let $\text{Id}_n$ denote the identity matrix in $\mathbb R^n$, then $(\text{Id}_n-\mathbf H)^2 = \text{Id}_n-\mathbf H$.

Proof. Follows directly from $\mathbf H^2 = \mathbf H$. $\square$

Definition. The linear regression is done with intercept if we impose that $x_{1,0}=\dots=x_{n,0}=1$. In other words, our prediction $\hat{\mathbf y}$ is an affine function of $(x_{\cdot, 1},\dots, x_{\cdot, p})$ and not necessarily a linear one (we are free to choose the intercept, i.e. the value of $\hat y$ when $(x_{\cdot, 1},\dots, x_{\cdot, p})=0$).

Lemma 2. If the linear regression is done with intercept, then $(1,\dots, 1)^\top\in\mathbb R^n$ is an eigenvector of $\mathbf H$ with eigenvalue $1$.

Proof. Since $\mathbf H\mathbf X=\mathbf X$, every column of $\mathbf X$ is an eigenvector of $\mathbf H$ with eigenvalue $1$. By Definition, the first column of $\mathbf X$ equals $(1,\dots, 1)^\top$ if fitting is done with intercept. $\square$

Remark. Lemma 2 is wrong in general if fitting is done without intercept.

Disclaimer. We assume from now on that fitting is done with intercept.

Proposition 1. Let $\mathbf 1_n$ denote the matrix in $\mathbb R^{n\times n}$ with every entry equal to $1$. Then $\left(\mathbf H-\frac 1n\mathbf 1_n\right)^2 = \mathbf H-\frac 1n\mathbf 1_n$.

Proof. We have $$\left(\mathbf H-\frac 1n\mathbf 1_n\right)^2 = \mathbf H^2 - \frac 1n\mathbf 1_n\mathbf H - \frac 1n\mathbf H\mathbf 1_n +\frac 1{n^2}\mathbf 1_n^2.$$ We have $\mathbf 1_n^2 = n \mathbf 1_n$. Since $\mathbf H=\mathbf H^\top$, by Lemma 2, $\mathbf 1_n\mathbf H = \mathbf 1_n = \mathbf H\mathbf 1_n$. This ends the proof. $\square$

Proposition 2. We have $$\sum_{k=1}^n (y_k-\bar y)^2 = \sum_{k=1}^n (\hat y_k-\bar y)^2+\sum_{k=1}^n(y_k-\hat y_k)^2.$$

Proof. Note that $$\sum_{k=1}^n (\hat y_k-\bar y)^2 = \left\lVert \hat{\mathbf y}-\frac 1n\mathbf 1_n \mathbf y\right\rVert^2 = \left\langle\left(\mathbf H-\frac 1n\mathbf 1_n\right)\mathbf y, \left(\mathbf H-\frac 1n\mathbf 1_n\right)\mathbf y\right\rangle.$$

Using $\mathbf H^\top=\mathbf H$ and Proposition 2, we get $$\left\langle\left(\mathbf H-\frac 1n\mathbf 1_n\right)\mathbf y, \left(\mathbf H-\frac 1n\mathbf 1_n\right)\mathbf y\right\rangle = \left\langle\mathbf y, \left(\mathbf H-\frac 1n\mathbf 1_n\right)\mathbf y\right\rangle.$$

Using an analogous argument together with Lemma 1a, we get $$\sum_{k=1}^n (y_k-\hat y_k)^2=\left\lVert(\operatorname{Id}_n-\mathbf H)\mathbf y\right\rVert^2 = \left\langle\mathbf y, (\operatorname{Id}_n-\mathbf H)\mathbf y\right\rangle.$$

In sum,

$$\sum_{k=1}^n (\hat y_k-\bar y)^2 + \sum_{k=1}^n (y_k-\hat y_k)^2 = \left\langle \mathbf y, \left(\operatorname{Id}_n-\frac 1n\mathbf 1_n\right)\mathbf y\right\rangle.$$

Since $\left(\operatorname{Id}_n-\frac 1n\mathbf 1_n\right)$ is symmetrical and $$\left(\operatorname{Id}_n-\frac 1n\mathbf 1_n\right)^2 =\operatorname{Id}_n-\frac 1n\mathbf 1_n,$$ we thus have $$\sum_{k=1}^n (\hat y_k-\bar y)^2 + \sum_{k=1}^n (y_k-\hat y_k)^2 = \left\lVert \left(\operatorname{Id}_n-\frac 1n\mathbf 1_n\right)\mathbf y\right\rVert^2 = \sum_{k=1}^n (y_k - \bar y)^2,$$ thus achieving a proof. $\square$

Corollary. If $R^2$ is well-defined, then $$R^2 = 1-\frac{\sum_{k=1}^n (y_k-\hat y_k)^2}{\sum_{k=1}^n (y_k-\bar y_k)^2}.$$

---

We will now prove the equality of $R^2$ and the correlation between $\mathbf y$ and $\hat{\mathbf y}$.

Notation. Let $$\text{Cov}(\mathbf y,\hat{\mathbf y})\overset{\text{Def.}}=\frac 1n\sum_{k=1}^n (y_k-\bar y)(\hat y_k-\bar{\hat y}),$$ as well as $$V(\mathbf y) \overset{\text{Def.}}=\frac 1n\sum_{k=1}^n (y_k-\bar y)^2$$ and $$V(\hat{\mathbf y}) \overset{\text{Def.}}=\frac 1n\sum_{k=1}^n (y_k-\bar{\hat y})^2.$$

Theorem. If $V(\hat{\mathbf y})\neq 0$, then $V(\mathbf y)\neq 0$, $R^2$ is well-defined and $$R^2 = \left(\frac{\text{Cov}(\mathbf y,\hat{\mathbf y})}{\sqrt{V(\mathbf y)V(\hat{\mathbf y})}}\right)^2.$$

Proof. We first note the following identities that can easily be verified using the machinery provided above:

• $V(\mathbf y) = \frac 1n\left\lVert\left(\text{Id}_n-\frac 1n\mathbf 1_n\right)\mathbf y\right\rVert^2 = \frac 1n\left\langle\mathbf y,\left(\text{Id}_n-\frac 1n\mathbf 1_n\right)\mathbf y\right\rangle $,
• $V(\hat{\mathbf y}) = \frac 1n\left\lVert\left(\mathbf H-\frac 1n\mathbf 1_n\right)\mathbf y\right\rVert^2= \frac 1n\left\langle\mathbf y,\left(\mathbf H-\frac 1n\mathbf 1_n\right)\mathbf y\right\rangle$,
• $\text{Cov}(\mathbf y,\hat{\mathbf y}) = \frac 1n\left\langle \left(\operatorname{Id}_n-\frac 1n\mathbf 1_n\right)\mathbf y,\left(H-\frac 1n\mathbf 1_n\right)\mathbf y\right\rangle=\frac 1n\left\langle\mathbf y,\left(H-\frac 1n\mathbf 1_n\right)\mathbf y\right\rangle=V(\hat{\mathbf y})$.

Assume from now on that $V(\hat{\mathbf y})\neq 0$. This is equivalent to $\mathbf H \left(\operatorname{Id}_n-\frac 1n\mathbf 1_n\right)\mathbf y\neq 0$, which implies in particular that $\left(\operatorname{Id}_n-\frac 1n\mathbf 1_n\right)\mathbf y\neq 0$ and thus $V(\mathbf y)\neq 0$.

From the above identities, we also have

$$\frac{\text{Cov}(\mathbf y,\hat{\mathbf y})}{\sqrt{V(\mathbf y)V(\hat{\mathbf y})}} = \sqrt{\frac{V(\hat y)}{V(y)}} = \sqrt{\frac{\sum_{k=1}^n (\hat y_k-\bar y)^2}{\sum_{k=1}^n (y_k-\bar y)^2}}=\sqrt{R^2}.$$ $\square$

Remark. Note that $V(\mathbf y)\neq 0$ does not always imply $V(\hat{\mathbf y})\neq 0$, consider for example $\mathbf x_1 = (-1, 0)$, $y_1 = 1$, $\mathbf x_2 = (1, 0)$, $y_2 = 0$ and $\mathbf x_3 = (1,1), y_3 = 1$, then $\beta=(2/3, 0)$, so that $V(\mathbf y)\neq 0$ but $V(\hat{\mathbf y})=0$.

Corollary. If $p=1$, (still working with intercept), the data of interest is the second column of $\mathbf X$, denoted by $\mathbf X_{:,1}$, which is a vector in $\mathbb R^n$, and in particular if $V(\mathbf X_{:,1})\neq 0$ and $\beta_1\neq 0$ (writing $\beta=(\beta_0,\beta_1)$), we have that $R^2$ is well-defined and the identity $$R^2 = \left(\frac{\text{Cov}(\mathbf y, \mathbf X_{:,1})}{\sqrt{V(\mathbf y)V(\mathbf X_{:,1})}}\right)^2$$ holds.

Proof. This follows directly from the Theorem and the identity $\hat{\mathbf y}= \vec\beta_0+\mathbf X_{:,1} \beta_1$, when writing $\beta=(\beta_0,\beta_1)$ and $\vec\beta_0=(\beta_0,\dots,\beta_0)^\top\in\mathbb R^n$. $\square$

---

[Footnote 1]. Other choices of a loss, such as the loss in absolute value or a cubed loss etc. can of course also be very useful, however quadratic loss is particularly pleasant to deal with mathematically thanks to its convexity and differentiability.
[Footnote 2] The name might be somewhat misleading, see for instance https://en.wikipedia.org/wiki/Explained_variation#Criticism.