0
$\begingroup$

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!

$\endgroup$
2
  • 1
    $\begingroup$ 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 find 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)$ " $\endgroup$ – Chris Cloud Oct 29 '20 at 8:45
  • $\begingroup$ So I guess my choosing $K^\top$ to be $E_{n-1,n}N_{\mathbf{1}_n}$ is the right thing to do here? $\endgroup$ – Chris Cloud Oct 29 '20 at 8:52
0
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.