In practice how is the random effects covariance matrix calculated in a mixed effects model?

Basically what I'm wondering is how different covariance structures are enforced, and how the values inside these matrices are calculated. Functions like lme() allow us to chose which structure we'd like, but I'd love to know how they are estimated.

Consider the linear mixed effects model $Y=X\beta+Zu+\epsilon$.

Where $u \stackrel{d}{\sim} N(0,D)$ and $\epsilon \stackrel{d}{\sim} N(0,R)$. Furthermore:

$Var(Y|X,Z,\beta,u)=R$

$Var(Y|X,\beta)=Z'DZ+R=V$

For simplicity we'll assume $R=\sigma^2I_n$.

Basically my question is: How exactly is $D$ estimated from the data for the various parameterizations? Say if we assume $D$ is diagonal (random effects are independent) or $D$ fully parameterized (case I'm more interested in at the moment) or any of the various other parameterizations? Are there simple estimators/equations for these? (That would no doubt be iteratively estimated.)

EDIT: From the book Variance Components (Searle, Casella, McCulloch 2006) I've managed to gleam the following:

If $D=\sigma^2_uI_q$ then then the variance components are updated and calculated as follows:

$\sigma_u^{2(k+1)} = \frac{\hat{\textbf{u}}^T\hat{\textbf{u}}} {\sigma_u^{2(k)}\text{trace}(\textbf{V}^{-1}\textbf{Z}^T\textbf{Z})}$

$\sigma_e^{2(k+1)} = Y'(Y-X{\hat{\beta}}^{(k)}-{Z}\hat{{u}}^{(k)})/n$

Where $\hat{\beta}^{(k)}$ and $\hat{{u}}^{(k)}$ are the $k$th updates respectively.

Is there general formulas when $D$ is block diagonal or fully parameterized? I'm guessing in the fully parameterized case, a Cholesky decomposition is used to ensure positive definiteness and symmetry.

• arxiv.org/pdf/1406.5823 (in press at Journal of Statistical Software) might be useful ... Mar 18, 2015 at 23:51

The Goldstein .pdf @probabilityislogic linked is a great document. Here's a list of some references that discuss your particular question:

Laird and Ware, 1982: Random-effects models for longitudinal data.

McCulloch, 1997: Maximum likelihood algorithms for generalized linear mixed models.

The SAS User Guide excerpt for the MIXED procedure has some great information about covariance estimation and many more sources (starting on page 3968).

There are numerous quality textbooks on longitudinal/repeated measures data analysis, but here's one that goes into some detail about implementation in R (from the authors of lme4 and nlme):

Pinheiro and Bates, 2000: Mixed-Effects Models in S and S-PLUS.

EDIT: One more relevant paper: Lindstrom and Bates, 1988: Newton-Raphson and EM Algorithms for linear mixed-effects models for repeated-measures data.

EDIT 2: And another: Jennrich and Schluchter, 1986: Unbalanced Repeated-Measures Models with Structured Covariance Matrices.

• I've had a look at Pinheiro and Bates, specifically Chapter 2 (on theory and computation) but I didn't seem to gleam anything on how the covariance structure is enforced and estimated? I shall go over it again shortly. I've got a few of those papers just sitting here, I'll definitely have to read them again. Cheers.
– dcl
Nov 30, 2011 at 3:31
• @dcl Looking back through Chapter 2 of P&B, I do see that they may be glossing over some of the steps you are interested in (they mention optimizing the log-likelihood wrt the covariance parameters but don't say how). That being said, Section 2.2.8 may be the section that addresses your question the best.
– user5594
Nov 30, 2011 at 3:54
• @dcl Added one more source that may help.
– user5594
Nov 30, 2011 at 4:07
• thanks for the links. I've had a browse through these papers in the past, some of them get quite technical for me. I'll have another browse through them now, but on first glance I can't seem get what I want out of them.
– dcl
Nov 30, 2011 at 5:01
• @dcl Sorry for the wall of links, but your question is one that a person can spend a couple full lectures discussing (it's a very good question that is kind of swept under the rug when first learning about mixed-effects models). Apart from swimming through the literature, one thing you could do is to look at the source code for lme4 and see how it deals with this estimation.
– user5594
Nov 30, 2011 at 5:24

Harvey Goldstein isn't a bad place to start.

As with most complex estimation methods, it varies with the software package. However, often what is done is in the following steps:

1. Pick an initial value for $D$ (say $D_0$) and $R$ (say $R_0$). Set $i=1$
2. Conditional on $D=D_{i-1}$ and $R=R_{i-1}$, estimate $\beta$ and $u$ and $\epsilon$. Call the estimates $\beta_i$ and $u_i$ and $\epsilon_i$.
3. Conditional on $\beta=\beta_i$ and $u=u_i$ and $\epsilon=\epsilon_i$, estimate $D$ and $R$. Call the estimates $D_i$ and $R_i$
4. Check for convergence. If not converged, set $i=i+1$ and return to step 2

One simple, and fast method is IGLS, which is based on iterating between two least squares procedures, and is described in detail in chapter two. Downside is that it doesn't work well for variance components close to zero.

• I know that this is the general method, but how are D and R estimated, what equations are used for the various structures? What are good initial values? I will check out the pdf now, cheers.
– dcl
Nov 30, 2011 at 3:24

The following article gives a closed form solution for D:

two more references that could be useful Variance Components by Searle Et al and Lynch and Walsh Genetics and Analysis of Quantitative Traits . The Lynch and Walsh book gives a step by step algorithm if I recall right