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

Question

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.

Answer 1

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

Harville, 1976: Extension of the Gauss-Markov Theorem to include the estimation of random effects.

Harville, 1977: Maximum likelihood approaches to variance component estimation and to related problems.

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.

Answer 2

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.

Answer 3

The following article gives a closed form solution for D:

• J. Shao, Mari Palta & Roger Qu, (1998), "Least squares estimation of regression parameters in mixed effects models with unmeasured covariates", Communications in Statistics - Theory and Methods, 27:6, 1487-1501, DOI:10.1080/03610929808832172

Answer 4

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