# Show why the estimate of variance component using REML is unbiased

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!

• I read it on Extending the Linear Model with R. Generalized Linear, Mixed Effects and Nonparametric Regression Models(Faraway). Supposing $y\sim\mathcal{N}(X\beta,\sigma^2K)$, from the book, it says "The idea (of REML) is to ﬁnd all independent linear combinations of the response, $k$, such that $k^\top X = 0$. Form matrix $K$ with columns $k$, so that: $K^\top y\sim \mathcal{N}(0, K^\top VK)$ " – Chris Cloud Oct 29 '20 at 8:45
• So I guess my choosing $K^\top$ to be $E_{n-1,n}N_{\mathbf{1}_n}$ is the right thing to do here? – Chris Cloud Oct 29 '20 at 8:52

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