# Why is the denominator $N-p-1$ in estimation of variance?

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)$

• 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.
– rka
Mar 27 '18 at 13:25
• I think this is the best answer. 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.

• 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.
– QED
May 12 '17 at 10:42
• $\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. 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}