Finding the unbiased variance estimator in high dimensional spaces

Question

The problem comes from linear regression. Assume the regression function is linear, i.e. $$ f(X) = \beta_0+\sum_{j=1}^pX_j\beta_j $$ .Given a set of training data $(x_1, y_1),\ldots,(x_N,y_N)$,we try to estimate the parameters $\beta$ by minimizing the residual sum of squares: $$\text{RSS}(\beta) = (\sum_{i=1}^N(y_i-\beta_0-\sum_{j=1}^px_{ij}\beta_j)^2$$ Here each $x_i=(x_{i1},\ldots,x_{ip})^T$ is vector in $\mathbb{R}^p$,$y_i\in\mathbb{R}$,and $\beta=(\beta_0,\ldots,\beta_{p+1})$. Let $\mathbf{X}$ denote matrix $(x_1,\ldots,x_N)^T$, $\mathbf{y}=(y_1,\ldots,y_N)$ and assume it has full column rank, it's easy to get $$\hat{\beta}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$$ Suppose $\mathbf{\hat{y}}=\mathbf{X}\hat{\beta}$, and $\mathbf{y}=\mathbf{X}\beta$ for some $\beta$, to pin down the sampling properties of $\hat{\beta}$, now assume that the observations $y_i$ are uncorrelated and have constant variance $\sigma^2$, and that the $x_i$ are fixed(non random). It's easy to see that $$\text{Var}(\hat\beta) = (\mathbf{X}^T\mathbf{X})^{-1}\sigma^2$$ ,so we turn to estimating the $\sigma^2$, and now comes the question: Let $\hat\sigma$ defined as below: $$\hat\sigma^2=\frac{1}{K}\sum_{i=1}^N(y_i-\hat{y_i})^2$$ It is said that only when $K=N-p-1$ did $\hat\sigma$ be an unbiased estimator of $\sigma$. But how do $p$(the dimension of $x$) be introduced in?

Answer 1

Well it might be an overkill but I think this proof is OK. I would use basic linear algebra tools. Starting with slight change in your notations:

Let $X$ denote the matrix $(x_1\; x_2 \;...\; x_N)^T$ where $x_i=(1 \; x_{i2} \; x_{i3} \; ... \; x_{ip})$. So now we have $p-1$ covariates. Our model is $$ y=X\beta+\epsilon $$ with $E(\epsilon)=0$ and $Var(\epsilon)=\sigma^2I$ Now, we can write $$ \hat{y}=X\hat{\beta}=P_Xy $$ where $P_X$ is the projection matrix onto the space spanned by the columns of $X$. That is, $$ P_X=X(X^TX)^{-1}X^T. $$ Note that $P_X$ is $n\times n$ matrix. Now, note that we can actually write $\sum_{i=1}^{N}(y_i-\hat{y}_i)^2$ as $||y-P_Xy||^2 $ with $||\cdot||$ being the $\ell_2$ norm of a vector. Next, $$ ||y-P_Xy||^2=||(I-P_X)y||^2=y^T(I-P_X)y $$ where the last equality holds since $I-P_X$ is a projection matrix. Note also that $(I-P_X)X=0$ and hence $y^T(I-P_X)y=\epsilon^T(I-P_X)\epsilon$ (using our model definition).

If we assume normality of the error, we could continue to derive the appropriate $\chi^2$ distribution. But you wrote nothing regarding the distribution so I am now turning to the expectation Recall that $I-P_X$ is $N\times N$ projection matrix. We can use the eigendecomposition of a matrix to write $$ (I-P_X)=U\Lambda U^T $$ where $UU^T=I$ and $\Lambda$ is a diagonal matrix with the diagonal being the eigenvalues of $I-P_X$.

Since $I-P_X$ is a projection matrix, all of its eigenvalues are equal to zero or one. How many "ones" there are? as the dimension of the subspace that $I-P_X$ projecting onto. Since this is the complementary of the subspace spanned by $P_X$, we get $n-p$ ones in the diagonal, and by using the fact that $UU^T=I$ we get the final answer which is $$ E(||(I-P_X)y||^2)= E(\epsilon^TU\Lambda U^T\epsilon)=\sum_{j=1}^n\lambda_j\sigma^2=(n-p)\sigma^2 $$ and as written before there are $n-p \; \lambda_j$'s that are equal to one and $p \; \lambda_j$'s that equals to zero ($\lambda_j, i=1,...,n$ are the eigenvalues of $I-P_X$). Some relevant Wikipedia pages for the linear algebra:

http://en.wikipedia.org/wiki/Projection_(linear_algebra) http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix http://en.wikipedia.org/wiki/Symmetric_matrix