I'm trying to use a very simple example to illustrate how REML makes the estimate of variance component unbiased:
Consider $X_1,\dots,X_n\overset{i.i.d.}{\sim}\mathcal{N}(\mu,\sigma^2)$, we denote one realization as $x_1,\dots,x_n$. We are interested in estimating $\sigma^2$. If we choose ML approach, we will have $\hat{\mu}_{\text{ML}}=\bar{x}=\sum_{i=1}^nx_i,\hat{\sigma}^2_{\text{ML}}=\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2$. Let $X^\top=[X_1,\dots,X_n]^\top$, so $X\sim\mathcal{N}(\mu\mathbf{1}_n,\sigma^2\mathbf{I}_n)$. We can transform $X$ into $N_{\mathbf{1}_n}X$, where $N_{\mathbf{1}_n}=\mathbf{I_n}-\frac{\mathbf{1}_n\mathbf{1}_n^\top}{n}$, now $N_{\mathbf{1}_n}X\sim\mathcal{N}(\mathbf{0}_n,\sigma^2N_{\mathbf{1}_n})$. However, we cannot write the p.d.f. of $N_{\mathbf{1}_n}X$ since $N_{\mathbf{1}_n}$ is not invertible.
Let
$$ E_{n-1,n}= \begin{bmatrix} 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix}_{n-1\times n} $$
Then $E_{n-1,n}N_{\mathbf{1}_n}X\sim\mathcal{N}(\mathbf{0}_{n-1},\sigma^2E_{n-1,n}N_{\mathbf{1}_n}E_{n-1,n}^\top)$, with
$$ (E_{n-1,n}N_{\mathbf{1}_n}E_{n-1,n}^\top)^{-1}= \begin{bmatrix} 2 & 1 & \cdots & 1 \\ 1 & 2 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 2 \end{bmatrix}_{n-1\times n-1} $$
Hence we can write down the likelihood for $E_{n-1,n}N_{\mathbf{1}_n}X$:
$$ \frac{1}{(2\pi)^{\frac{n-1}{2}}(\sigma^2)^{\frac{n-1}{2}}[\det(E_{n-1,n}N_{\mathbf{1}_n}E_{n-1,n}^\top)]^{\frac{1}{2}}}\exp(-\frac{1}{2 \sigma^2}X^\top N_{\mathbf{1}_n}E_{n-1,n}^\top (E_{n-1,n}N_{\mathbf{1}_n}E_{n-1,n}^\top)^{-1}E_{n-1,n}N_{\mathbf{1}_n}X) $$
Take log:
$$ \text{const}-\frac{n-1}{2}\log\sigma^2-\frac{1}{2\sigma^2}\color{red}{X^\top N_{\mathbf{1}_n}E_{n-1,n}^\top (E_{n-1,n}N_{\mathbf{1}_n}E_{n-1,n}^\top)^{-1}E_{n-1,n}N_{\mathbf{1}_n}X} $$
I'm stuck with showing $\color{red}{X^\top N_{\mathbf{1}_n}E_{n-1,n}^\top (E_{n-1,n}N_{\mathbf{1}_n}E_{n-1,n}^\top)^{-1}E_{n-1,n}N_{\mathbf{1}_n}X}=SSE=\sum_{i=1}^n(x_i-\bar{x})^2$, which doesn't seem very obvious to me.
Or could someone suggest me a more elegant way to show $\hat{\sigma}^2_{\text{REML}}=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2$ in this example?
Any help appreciated!