# How are calculations done for REML?

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.

• Possible duplicate of What is "restricted maximum likelihood" and when should it be used? – Tim Apr 4 '16 at 20:12
• @Tim I've already read that, as you can clearly see from my link. Can you explain to me how any of these answers explain how to do REML? – Clarinetist Apr 4 '16 at 20:13
• I suggest you edit your question to make it more clearly distinct from the one you intend it not to be a duplicate of. This would be helpful for you (it will make it less likely that it is closed as a duplicate) and it be helpful for future readers (at the moment, if they search for the topic, how would they distinguish between this thread and the other)? I'm voting to leave open, but I strongly suggest you make some edits - particularly to the title, which is completely indistinguishable from the other question. Perhaps something like "How is REML estimation calculated?" – Silverfish Apr 4 '16 at 20:45
• @Silverfish See the edit. – Clarinetist Apr 4 '16 at 22:37
• Thanks for your edit. The new title is definitely more distinct. – Silverfish Apr 4 '16 at 22:50

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.