Show why the estimate of variance component using REML is unbiased

Question

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!

Answer 1

I found this REML estimation of variance components lecture note extremly helpful!

I'm inspired to do eigendecomposition with $N_{\mathbf{1}_n}$ so I get

$$ N_{\mathbf{1}_n}=B \begin{bmatrix} 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0\\ 0 & 0 & \cdots & 0 & 0 \end{bmatrix} B^\top, $$

where $B^\top B=\mathbf{I}_n$. If we take first $n-1$ rows of $B^\top$ and denote them as $C^\top$, then we have

$$ C^\top C=\mathbf{I}_{n-1},CC^\top=N_{\mathbf{1}_n} $$

Note $C^\top X\sim\mathcal{N}(\mu C^\top\mathbf{1}_n,\sigma^2C^\top C)$. Since $\mathbf{0}_n=N_{\mathbf{1}_n}\mathbf{1}_n=CC^\top \mathbf{1}_n$, we can left multiply both sides by $\mathbf{1}_n^\top$ to get $\mathbf{1}_n^\top\mathbf{0}_n=0=\lVert C^\top\mathbf{1}_n\rVert$, which implies $C^\top\mathbf{1}_n=\mathbf{0}_{n-1}$!

So we have $C^\top X\sim\mathcal{N}(\mathbf{0}_{n-1},\sigma^2\mathbf{I}_{n-1})$. Writing down the log-likelihood:

$$ \text{const}-\frac{n-1}{2}\log\sigma^2-\frac{1}{2\sigma^2}X^\top N_{\mathbf{1}_n}X $$

I suppose this is a so much better way than my first attempt :)