I was recently going through the book Elements of Statistical Learning by Tibshirani et.al. In this book, while explaining the ordinary least squares model, the authors state that assume that $y_i \epsilon \mathbb{R}$ represents the observed variables, $\hat{y_i}$ represents the model output and $\mathbf{x_i} \epsilon \mathbf{R}^{p+1}$ represent the inputs. If the $y_i$s are assumed to be uncorrelated and have constant with variance $\sigma$, then the unbiased estimate of variance is $\hat{\sigma} = \frac{1}{\left (N-p-1 \right)}\sum\left( y_i - \hat{y_i} \right)^2$, summation being done from $i=1$ to $i=N$. Note that $p$ has been used here to denote the dimensionality of $\mathbf{x_i}$s. My question is why is the factor in the denominator $N-p-1$ while estimating the variance of $y_i$s i.e. $\hat{\sigma}$ ? From my understanding if the $y_s$s are real numbers that have constant variance, the factor should be equal to $N-1$.


The current accepted answer is flawed, as it implicitly assumes that the error of the model $\varepsilon$ is Gaussian (otherwise you need not have $\sum(y_i-\hat{y}_i)^2\sim\sigma^2\chi^2_{N-p-1}$).

Here's a proof with the general assumption that $\varepsilon$ has mean $0$ and variance $\sigma^2 I_N$.

First note that $\sum(y_i-\hat{y}_i)^2=\|y-X\hat\beta\|^2$.

We have $$\begin{align} y-X\hat\beta &= X\beta +\varepsilon -X(X^TX)^{-1}X^T(X\beta +\varepsilon)\\ &=X\beta +\varepsilon - X\beta -X(X^TX)^{-1}X^T\varepsilon\\ &= (I_N-H)\varepsilon\end{align}$$ where $H=X(X^TX)^{-1}X^T$ is the hat matrix. It's easy to check that $H^T=H$ and $H^2=H$ (indeed the hat matrix is merely the orthogonal projection on $\operatorname{Im}X$).

Hence $\begin{aligned}[t]E( \|y-X\hat\beta\|^2) &= E(\varepsilon^T(I_N-H)^T (I_N-H)\varepsilon)=E(\varepsilon^T(I_N-H)\varepsilon) \end{aligned}$

Note that $\varepsilon^T(I_N-H)\varepsilon=\sum_{i,j} \varepsilon_i\varepsilon_j (\delta_{ij}-H_{ij})$, thus $$E(\varepsilon^T(I_N-H)\varepsilon)=\sum_{i,j} \sigma^2\delta_{ij} (\delta_{ij}-H_{ij})=\sigma^2(N-\operatorname{tr}H)$$

Note that $\operatorname{tr}H =\operatorname{tr}(X(X^TX)^{-1}X^T)=\operatorname{tr}(X^TX(X^TX)^{-1})=\operatorname{tr}(I_{p+1})=p+1 $

Putting everything together, $E( \|y-X\hat\beta\|^2)=\sigma^2(N-p-1)$

  • $\begingroup$ Just before the variance equation in ESL it's written ".. we now assume that the observations $y_i$ are uncorrelated and have constant variance $σ^2$ , and that the $x_i$ are fixed (non random)".. I don't see where they make the claim that $Y = f(X) + \epsilon$. This was done only in chapter two. I was under the impression that the assumptions in this chapter are more relaxed. Unfortunately, without this assumption I could only get as far as this: $E[\|y-X\hat\beta\|^2] = (I-H)^T E[y y^T] (I-H)$. We could use $\sigma^2 I = E[y y^T] - E[y] E[y]^T$ but we don't know the mean. $\endgroup$
    – rka
    Mar 27 '18 at 13:25
  • $\begingroup$ I think this is the best answer. $\endgroup$
    – maplemaple
    Dec 15 '19 at 20:23

You can show that $\sum(y_i-\hat{y}_i)^2\sim\sigma^2\chi^2_{N-p-1}$. As expectation of a $\chi^2_{N-p-1}$ is $(N-p-1)$. Hence $\mathbb{E}(\frac{1}{N-p-1}\sum(y_i-\hat{y}_i)^2)=\sigma^2$.

$N-p-1$ is in the denominator to make the estimator unbiased.

  • 1
    $\begingroup$ you can convert $\sum(y_i-\hat{y}_i)^2/\sigma^2$ into sum of squares of $N-p-1$ independent standard normal variables by Gram-Schmidt orthogonalization. $\endgroup$
    – QED
    May 12 '17 at 10:42
  • 3
    $\begingroup$ $\sum(y_i-\hat{y}_i)^2\sim\sigma^2\chi^2_{N-p-1}$ requires that the error term be Gaussian, which may not be true. $\endgroup$ Jul 1 '17 at 14:30

To answer the question without using the Gaussian assumption, nor the additive model assumption as it has been done in previous answer, here is my take:

The assumption made are that the observations $y_i$ are uncorrelated and have mean zero, constant variance $\sigma^2$, and that the $x_i$ are fixed (non random).

From the previous section we have that $\hat{y} = X(X^TX)^{-1}X^Ty = Hy$ where the hat matrix can be shown easily to satisfy $H^T=H^2=H$.

Also we can rewrite $\sum_{i=1}^N(y_i - \hat{y_i})^2$ as $(y-\hat{y} )^T(y-\hat{y} )$ and thus it comes that

\begin{align*} \mathbb{E}[\sum_{i=1}^N(y_i - \hat{y_i})^2] &= \mathbb{E}[(y-Hy)^T(y-Hy)] = \mathbb{E}[((I-H)y)^T((I-H)y)] \\ &= \mathbb{E}[y^T(I-H)(I-H)y] = \mathbb{E}[y^T(I-H)y] \\ &= \mathbb{E}[\sum_{i,j}y_iy_j(\delta_{ij}-H_{i,j})] \\ &= \sum_{i,j}\mathbb{E}[y_iy_j\delta_{ij}]- \sum_{i,j}\mathbb{E}[y_iy_jH_{i,j}] \\ &= \sum_{i}\mathbb{E}[y_i^2]- \sum_{i,j}\mathbb{E}[y_iy_j]H_{i,j} \end{align*}

The hat matrix has been removed from the expectation since it is fixed as the $x_i$ are.

In addition the $y_i$ are uncorrelated of mean zero and variance $\sigma^2$ thus using the cyclic property of the trace: \begin{align*} \mathbb{E}[\sum_{i=1}^N(y_i - \hat{y_i})^2] &= \sum_{i}\sigma^2 - \sum_{i}\sigma^2H_{i,i} = \sigma^2[N - trace(H)] \\ &= \sigma^2[N - trace(X(X^TX)^{-1}X^T)] \\ &= \sigma^2[N - trace(X^TX(X^TX)^{-1})] \\ &= \sigma^2[N - p - 1] \end{align*}

So finally we have that an unbiased estimator $\hat{\sigma}^2$ of $\sigma^2$ is \begin{align} \mathbb{E}[\hat{\sigma}^2] &= \mathbb{E}[\frac{1}{N - p - 1} \cdot \sum_{i=1}^N(y_i - \hat{y_i})^2] = \sigma^2 \end{align}


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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