How are calculations done for REML?

Question

I've read a few questions on this site (e.g., https://stats.stackexchange.com/a/48676/46427), but they rarely go beyond intuitive explanations. I am particularly interested with how to calculate an estimate via REML, not the intuition behind REML.

I know what maximum likelihood estimation is; basically, take your likelihood function $L(\boldsymbol{\theta} \mid \mathbf{y}) = f(\mathbf{y})$ and find the values of $\boldsymbol{\theta}$ which maximize $L$; that is your maximum-likelihood estimator.

Now we're talking about REML in a linear models class... something about how you take $n - r(\mathbf{X})$ linearly independent row vectors ($\mathbf{X}$ the design matrix, I suppose...) such that when multiplied by $\mathbf{X}$, you get the vector $\mathbf{0}^{T}$, use that as your data, compute a likelihood (or loglikelihood) based on that data, and do MLE on this (log)likelihood?

Searching online hasn't helped me very much. Could someone explain REML to someone who's very familiar with MLE? I cannot take an intuitive explanation here, like in the link I have in this post; I actually need to be able to do it.

Answer 1

Here is a simple example with calculations that shows the idea. We work in the linear model $Y = X\beta + e, e\sim N(0, \Sigma(\theta))$, where $Y$ is the $N \times 1$ response vector, $X$ the $n\times p$ design matrix, and $\theta$ parametrizes the covariance matrix. Suppose interest lies in estimating $\theta$.

Assume $X$ has full column rank and let $A$ be an (any!) $n \times (n-p)$matrix with orthonormal columns such that $A'X = 0$, then Harville, 1974 showed that the likelihood for $A'Y$ is:

\begin{align} RL(\theta; Y)=(2\pi)^{-(n-p)/2}\vert X'X\vert^{1/2}\vert X' \Sigma^{-1}X\vert^{-1/2}\vert \Sigma \vert^{-1/2}\exp\left(-\frac{1}{2}r'\Sigma^{-1}r \right), \end{align}

where $r = QY:=\left(I - X(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}\right)Y$. REML entails maximization of this likelihood over $\theta$ only, instead of maximizing the usual likelihood:

\begin{align} L(\theta, \beta; Y) = (2\pi)^{-n/2}\vert \Sigma\vert^{-1/2}\exp\left(-\frac{1}{2}[Y - X\beta]'\Sigma^{-1}[Y-X\beta]\right), \end{align}

which gives MLEs for $(\theta, \beta)$.

To get some easy calculations, suppose $X = 1_n$, a vector of ones, and $\Sigma = \sigma^2I$. The log-restricted likelihood, up to terms that does not depend on $\sigma^2$, is in this particular case:

\begin{align} RL(\theta ; Y) &\propto \log\vert \sigma^{2} n\vert^{1/2}+\log\vert\sigma^2I\vert^{-1/2} -\frac{1}{2\sigma^2}SSE \\ &\propto -\frac{n-1}{2}\log(\sigma^2) -\frac{1}{2\sigma^2}SSE, \end{align}

which is maximized at $\sigma^2 = SSE/(n-1)$. Note the following which was used to get the restricted log-likelihood:

• When $\Sigma = \sigma^2I$, $Q = I - X(X'X)^{-1}X'$, which does not depend on $\sigma^2$.
• When $X = 1_n$, $X'X = n$
• $SSE = (Y - 1_n\bar{Y})'(Y - 1_n\bar{Y})$, since $Q = I - 1_n1_n'/n$.

It is conceptually easy to generalize to a general linear mean function $X\beta$ and general covariance matrix $\Sigma(\theta)$; $\hat{\theta} = \arg \max RL(\theta; Y)$. However, the actual computations will generally require a little more work.