Bias in the MLE of variance component in a multivariate Gaussian? Given an $n$-vector $y$ (responses) and a design matrix $X$, I wish to fit them with a simple linear regression model
$$y=X\beta+e,$$
where $e\sim\mathcal{N}(0, \sigma^2I)$. Then, we have 
$$y\sim\mathcal{N}(X\beta, \sigma^2I).$$
Then the maximum likelihood estimations (MLE) of $\beta$ and $\sigma^2$ are just
$$\hat\beta=(X^TX)^{-1}X^Ty$$ and 
$$\hat{\sigma^2}=\frac{(y-X\hat\beta)^T(y-X\hat\beta)}{n}.$$
I understand that $\sigma^2$'s estimator is biased in that $E\{\hat{\sigma^2}\}\neq\sigma^2$.
On Page 4 of this lecture notes, the author claims $E\{\hat{\sigma^2}\}=\frac{n-r}{n}\sigma^2$ without specifying what $r$ is. I guess $r$ is the DoF loss while estimating $\beta$, i.e., the dimension of $\beta$. But how do we derive it?
I tried to prove it myself but got confused at "over which distribution are we taking the expectation?"
 A: Yes, $r$ is indeed the length/dimension of $\beta$. Define the 'hat' matrix $H=X(X^TX)^{-1}X^T$ so that $Hy=X\hat\beta=\hat y$. Then $E((y-X\hat\beta)^T(y-X\hat\beta))=E((y-Hy)^T(y-Hy))=E(y^T(I-H)^T(I-H)y)$. Now $I-H$ is an orthogonal projection so $(I-H)^T(I-H)=I-H$ and hence the expectation above reduces to $E(y^Ty)-E((Hy)^T(Hy))$. Given that $y\sim\mathcal{N}(X\beta, \sigma^2I)$, it follows that $E(y^Ty)=\sigma^2trace(I)+(X\beta)^T(X\beta)=n\sigma^2+(X\beta)^T(X\beta)$. It also follows that the transformed variable $Hy\sim\mathcal{N}(HX\beta,\sigma^2HH^T)=\mathcal{N}(X\beta,\sigma^2H)$. Hence $E((Hy)^T(Hy))=\sigma^2trace(H)+(X\beta)^T(X\beta)$. Because $H$ is an orthogonal projection, its trace is equal to its rank, which equals the rank of $X$, which is just $dim(\beta)$. So, combining the obtained equations for $E(y^Ty)$ and $E((Hy)^T(Hy))$ leads to $E((y-X\hat\beta)^T(y-X\hat\beta))=(n-dim(\beta))\sigma^2$.
A: Just make our life easier, suppose $X$ is just one random variable and your model has no intercept then the MLE of $\sigma^2$ is $\hat{\sigma^2}=\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n}$ (This is just a simpler case for $\hat{\sigma^2}=\frac{(y-X\hat\beta)^T(y-X\hat\beta)}{n}$ and $\hat\mu=X\hat\beta=\bar{X}$)
We know  $\hat{\sigma^2}=\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n}$ is a biased estimator of $\sigma^2$, the unbiased one should be  $\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n-1}$
Therefore, $E(\frac{n}{n-1} \hat{\sigma^2})=E(\frac{n}{n-1}*\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n})=\sigma^2 $
$\Rightarrow E(\hat{\sigma^2)}=\frac{n-1}{n}\sigma^2$
you can see $r=1$ here.
For multivariate situation I guess it should be the same, hope someone else can help.
