# Restricted Maximum Likelihood (REML) Estimate of Variance Component

Let, $$\mathbf y_i = \mathbf X_i\mathbf\beta + \mathbf Z_i\mathbf b_i+ \mathbf\epsilon_i,$$

where

$\mathbf y_i\sim N(\mathbf X_i\mathbf\beta, \Sigma_i=\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i'),$

$\mathbf b_i\sim N(\mathbf 0, \mathbf G),$

$\mathbf\epsilon_i\sim N(\mathbf 0, \sigma^2\mathbf I_{n_i}),$

$\mathbf y_i$ is a $n_i\times 1$ vector of response for $i^{th}$ individual at $1,2,\ldots, n_i$ time points, $\mathbf X_i$ is a $n_i\times p$ matrix, $\mathbf \beta$ is a $p\times 1$ vector of fixed effect parameters, $\mathbf Z_i$ is a $n_i\times q$ matrix, $\mathbf b_i$ is a $q\times 1$ vector of random effects, $\mathbf \epsilon_i$ is a $n_i\times 1$ vector of within errors, $\mathbf G$ is a $q\times q$ covariance matrix of between-subject, $\sigma^2$ is a scalar.

Note that, $\mathbf X_i$, $\mathbf Z_i$, and $\mathbf G$ do NOT involve $\sigma^2$.

Now I have to find out the Restricted Maximum Likelihood (REML) Estimate of $\sigma^2$, that is,

$$\hat\sigma^2_R = \frac{1}{N_0-p}\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'(\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta),\ldots (1)$$

where $N_0 = \sum_{i=1}^{N}n_i$.

So first I wrote the Restricted Maximum Log-Likelihood :

$$l_R \propto -\frac{1}{2}\sum_{i=1}^{N}\log\det(\Sigma_i)-\frac{1}{2}\sum_{i=1}^{N}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)-\frac{1}{2}\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'\Sigma_i^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta).$$

Then I have to differentiate log-likelihood, $l_R$, with respect to $\sigma^2$ and equate it to zero, i.e.,

$-\frac{1}{2}\frac{\partial}{\partial\sigma^2}\{\sum_{i=1}^{N}\log\det(\Sigma_i)+\sum_{i=1}^{N}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)+\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'\Sigma_i^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta)\}|_{\sigma^2=\hat\sigma^2_R}=0.$

But basically I can't differentiate,

$\frac{\partial}{\partial\sigma^2}\log\det(\Sigma_i)=\frac{\partial}{\partial\sigma^2}\log\det(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')$

$\frac{\partial}{\partial\sigma^2}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)=\frac{\partial}{\partial\sigma^2}\log\det(\mathbf X_i'(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}\mathbf X_i)$ and

$\frac{\partial}{\partial\sigma^2}\Sigma_i^{-1}= \frac{\partial}{\partial\sigma^2}(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}$.

Here is an answer to differentiate such derivative but I can't combine it to get the REML estimate $\hat\sigma^2_R$ in equation $(1)$.

How can I get the REML estimate $\hat\sigma^2_R$ in equation $(1)$?

• If all you need is the solution for a particular model, you can use AD Model Builders random effects module to code up this simple model and solve for the parameters estimates including $\sigma_R$. The main difference between what you are sketching out here and ADMB is that ADMB uses automatic differentiation to numerically calculate all the necessary derivatives to maximize the llikelihood. – dave fournier May 28 '16 at 13:29
• I highly recommend Variance Components by Searle et al. It goes in depth into how to obtain derivatives. – Nicholas Mancuso Feb 16 '18 at 23:11

## 2 Answers

NB. I simplify notation somewhat and do not use bold typesetting.

The following rules for matrix differentials are useful:

\begin{align} d\log \vert A\vert &= \mathrm{tr}(A^{-1}dA) \\ dA^{-1} & = -A^{-1}(dA)A^{-1}. \end{align}

A good source for such rules and how to derive them is Magnus and Neudecker. That text also explains how to go between differentials and derivatives. You can also refer to this document by Minka. The Matrix Cookbook is another standard reference for such rules. In addition to derivatives, this answer also makes repeated use of cyclic invariance and linearity of the trace operator.

Based on these rules we get the following expressions for the differentials in question: $$d\Sigma_i^{-1} = -\Sigma_i^{-1}d\Sigma_i \Sigma_i^{-1} = -\Sigma^{-2}d\sigma^2,$$ $$d\log \vert \Sigma_i \vert = \mathrm{tr}\Sigma_i^{-1}d\Sigma_i = \mathrm{tr}\Sigma_i^{-1}Id\sigma^2,$$

and

\begin{align} d\log\vert X_i'\Sigma_i^{-1}X_i\vert &= \mathrm{tr}[(X_i'\Sigma_iX)^{-1}X_i'd\Sigma_i^{-1}X_i] \\ &= -\mathrm{tr}[(X_i'\Sigma_iX)^{-1}X_i'(\Sigma_i^{-1}d\Sigma_i \Sigma_i^{-1})X_i]\\ &= -\mathrm{tr}[\Sigma_i^{-1}X_i(X_i'\Sigma_iX)^{-1}X_i'\Sigma_i^{-1}Id\sigma^2] \end{align}

Thus, the gradient of the log-restricted likelihood up to scaling and additive terms not involving $\sigma^2$ is:

\begin{align} \sum_i \frac{\partial}{\partial \sigma^2}\left[-\log \vert \Sigma_i \vert - \log \vert X_i'\Sigma_i^{-1}X_i \vert - (y_i - X_i\beta)'\Sigma_i^{-1}(y_i - X_i\beta) \right] = \\ \sum_i \mathrm{tr}\left[-\Sigma_i^{-1} + [\Sigma_i^{-1}X_i(X_i'\Sigma_i^{-1}X)^{-1}X_i'\Sigma_i^{-1}] + (y_i - X_i \beta)'\Sigma_i^{-2}(y_i - X_i\beta)\right]. \end{align}

Since the $\Sigma_i$ are (I assume) positive definite, and in particular full rank, setting the gradient to zero is equivalent to: $$\sum_i \mathrm{tr}\left[I - [X_i(X_i'\Sigma_i^{-1}X)^{-1}X_i'\Sigma_i^{-1}] - \Sigma_i^{-1}(y_i - X_i\beta)(y_i - X_i \beta)'\right] = 0.$$

Now notice that $Q_i := I - [X_i(X_i'\Sigma_i^{-1}X)^{-1}X_i'\Sigma_i^{-1}]$ is a projection matrix onto the orthogonal complement of the column space of $X_i$. Thus, its trace equals the dimension of the space onto which it is projecting: $n_i - p$. Factoring out $\sigma^2$ then gives the REML estimate of $\sigma^2$ as the solution to $$\sigma^2 = \sum_i (y_i - X\beta)'[I + Z_iGZ_i' / \sigma^2]^{-1}(y_i - X_i\beta) / (N_0 - p),$$

assuming the optimization problem is convex.

This is different from your expression. It's possible I made a mistake somewhere and, in that case, hopefully you or another reader can spot it. After reading some more though, my impression is that mixed models are commonly parameterized in terms of $\tilde{G} = G/\sigma^2$. Parameterizing the first order condition derived here in terms of $\tilde{G}$ would make it agree with your first order condition.

• These kinds of unstructured derivative calculations will work for really simple models. However as soon as the model gets a little more complicated you will experience a bit of grief. How would you for example modify this if for robustness one wanted to assume that the $\epsilon_i$ have a student t distribution with the degrees of freedom as another parameter. If this had been written in ADMB the modifications to the code would be pretty simple. – dave fournier May 28 '16 at 23:03
• With unstructured, do you mean you would like me to structure the answer better? I'd be happy to accommodate that. Anything in particular you had in mind? I wouldn't modify this for another distributional assumption; the question states the objective function of interest. – ekvall May 29 '16 at 0:45
• I am referring to the techniques for calculating derivatives known as automatic differentiation or the differentiation of algorithms. Using these ideas one sets up a general scheme for calculating derivatives which many advantages, one of which is that the user does not have to calculate the derivatives by hand. For some reason these techniques are not nearly as well known in computational statistics as they should be. – dave fournier May 29 '16 at 13:30
• That sounds interesting @davefournier. I in fact looked at your (?) software for one of my own research projects. In this particular case, the OP asked for a derivation of a stated "closed form" result, and how to differentiate a few matrix terms. I believe I answered that. Still, I'd personally be very interested in reading an answer making use of computational techniques if you want to write one up. – ekvall May 29 '16 at 13:45
• I know. You answered his question correctly, but he asked the wrong question. This is fairly common. With the vast increase in computational power what is needed in general is the ability to code up models that are much more general and ambitious than the usual glmm's the R world has saddled us with. As I said one could replace the normal $\epsilon_i$ with student t or the normal random effects with some kind of robust random effects distribution. That's just a tiny bit of what could be done. – dave fournier May 29 '16 at 14:13

This is the AD Model Builder code for the model. This version will estimate the fixed effects.

DATA_SECTION
init_number N   // number of individuals
init_int p
init_int q
init_vector n(1,N) // number of time points
init_3darray X(1,N,1,n,1,p)
init_3darray Z(1,N,1,n,1,q)
init_matrix obsy(1,N,1,n)
PARAMETER_SECTION
init_number log_sigma
!! int q1=q*(q-1)/2;
init_vector chcoff(1,q1)
init_vector logs(1,q)
init_vector beta(1,p)
random_effects_matrix eps(1,N,1,q)
objective_function_value f
PROCEDURE_SECTION
for (int i=1;i<=N;i++)
{
f+=f1(i,eps(i),beta,log_sigma,logs,chcoff);
}
SEPARABLE_FUNCTION   dvariable f1(int i,const dvar_vector& epsi,const dvar_vector beta,const prevariable & log_sigma,const dvar_vector& logs,const dvar_vector& chcoff)
dvariable sigma=exp(log_sigma);
dvar_vector s=exp(logs);
dvariable var=square(sigma);
dvar_matrix ch(1,q,1,q);
ch.initialize();
// ch is the choleski decomp of G
ch(1,1)=s(1);
int offset=0;
for (int j=2;j<=q;j++)
{
ch(j)(1,j-1)=chcoff(1+offset,j-1+offset).shift(1);
ch(j,j)=1.0;
ch(j)/=norm(ch(j));
offset+=j-1;
ch(j)*=s(j);
}
dvar_vector predy=X(i)*beta+Z(i)*(ch*epsi);
f+= 0.5*n(i)*log(var)+0.5*norm2(obsy(i)-predy)/var;
f+= 0.5*norm2(epsi);


The REML version is calculated by integrating over the fixed effects. This involves changing one line in the above code.

 init_vector beta(1,p)


is changed to

 random_effects_vector beta(1,p)


Now suppose I wanted to modify this code to use the student t distribution. The line

f+= 0.5*n(i)*log(var)+0.5*norm2(obsy(i)-predy)/var;


is changed to

f-=ld_student(obsy(i),predy,var,dof);


and the lines

init_number log_dof_coff

dvariable dof=1.0+exp(log_dof_coff);


are added to the code in the appropriate places to include the degrees of freedom parameter.
Note that while there is a closed form solution for integrating over the parameters for the normal distribution, this is not true for the student. The integration is carried out via the Laplace approximation.