Covariance Structures

Introduction

The current version of lme4 offers four covariance classes/structures: Covariance.us (unstructured), Covariance.diag (diagonal), Covariance.cs (compound symmetry), and Covariance.ar1 (autoregressive order 1). The syntax for use is similar to the glmmTMB package (see covariance structures in glmmTMB), although the results are slightly different.

This vignette provides a detailed mathematical exposition of the new machinery, along with brief notes on the appropriate syntax. In a separate document, we provide comparisons of results between glmmTMB and lme4 (the computations are too slow to satisfy CRAN timing requirements; if we tried to store and re-load saved fits, they would be too big to satisfy CRAN package size requirements).

Background

For exact details of this structure, refer to the lmer vignette (Bates et al. 2015). We provide a quick summary in this vignette.

In matrix notation, a linear mixed model can be represented as: $$\mathbf{y} = \mathbf{X} \boldsymbol \beta + \mathbf Z \mathbf{b} + \boldsymbol{\epsilon}$$ Where $\mathbf{b}$ represents an unknown vector of random effects, with $\mathbf b \sim \textrm{MVN}(\mathbf 0, \boldsymbol \Sigma)$.

We create the relative co-factor $\mathbf \Lambda_{\mathbf \theta}$ which is a $q \times q$ block diagonal matrix that depends on the variance-component parameter vector $\boldsymbol \theta$. Let $\sigma$ be the scale parameter of the variance of a linear mixed model. In lme4, the variance-covariance matrix is constructed by: $${\mathbf \Sigma}_{\mathbf{\theta}} = \sigma^2 {\mathbf \Lambda}_{\mathbf \theta} {\mathbf \Lambda}_{\mathbf \theta}^{\top},$$

For generalized linear mixed models, $\mathbf \Lambda_\theta$ instead represents the unscaled Cholesky factor; that is, the scaling term $\sigma^2$ is omitted from the equation above.

The major difference between the four covariance classes (Covariance.us (unstructured), Covariance.diag (diagonal), Covariance.cs (compound symmetry), and Covariance.ar1 (autoregressive order 1)) is in the construction of the the relative Cholesky factor $\mathbf \Lambda_{\mathbf \theta}$.

Covariance Structures

Suppose that a particular random-effects term has $p_i$ varying effects (i.e., the number of columns of the model matrix constructed from the left-hand side f of the random-effects formula (f|g). If the grouping variable (g) has $\ell_i$ levels and $q_i = \ell_i p_i$. In general the covariance matrix for the term, $\mathbf \Sigma_i$, is a $q \times q$ block diagonal matrix composed of homogeneous $p_i \times p_i$ blocks; the (relative) Cholesky factor $\boldsymbol \Lambda_i$ has a similar structure. In what follows we will describe the construction of these template matrices, for either the covariance matrix or the Cholesky factor.

The unstructured covariance, which is the default in lme4, of size $p \times p$ has the following form: $$\mathbf{\Sigma} = \begin{bmatrix} \sigma^{2}_{1} & \sigma_{12} & \dots & \sigma_{1p} \\ \sigma_{21} & \sigma^{2}_{2} & \dots & \sigma_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{p1} & \sigma_{p2} & \dots & \sigma^{2}_{p} \\ \end{bmatrix}$$

To specify unstructured (general positive-semidefinite) covariances for a particular term, specify either (f|g) (as in lme4 < 2.0) or us(f|g).

The next three covariance structures can either be heterogeneous or homogeneous. If we have a homogeneous covariance structure (hom = TRUE), then we assume $\sigma_{1} = \sigma_{2} = \dots = \sigma_{p}$.

The diagonal covariance has the following form: $$\mathbf{\Sigma} = \begin{bmatrix} \sigma^{2}_{1} & 0 & \dots & 0 \\ 0 & \sigma^{2}_{2} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \sigma^{2}_{p} \\ \end{bmatrix}$$ By default, we assume a heterogeneous diagonal covariance structure.

To specify a diagonal covariance structure, use diag(f|g) (or diag(f|g, hom = FALSE)). Unlike the unstructured covariance matrix, which is invariant to the parameterization of the terms in the varying effects model f, the diagonal covariance matrix (and the others described below) depend on the contrasts assigned to factor variables and on the inclusion of an intercept term. For example, suppose that f represents a factor variable with $n$ levels, and that the standard treatment contrasts are used. The template matrices will be $n \times n$.

• If the intercept is not included (i.e. the term is specified as (0+f|g) or (f-1|g)), then the variances $\sigma^2_1, \dots, \sigma^2_n$ will represent variances of observations in each level of the factor around the population mean;
• if the intercept is included (the term is specified as (f|g) or (1+f|g)), $\sigma^2_1$ will represent the among-cluster variance of observations from the first (baseline) level of f; $^2_2, $\sigma^2_n$ will represent variances of differences in expected values between observations from levels 2 to $n$ and those in the baseline level.

If a different, non-default contrast is used, then the variances describe variation in the corresponding parameters. For example, with successive-difference contrasts (MASS::contr.sdif()) $\sigma^2_1$ would estimate among-cluster variance for observations at the first level, $\sigma^2_2$ would estimate variance in the difference between level 2 and level 1, and so forth.

Similar points apply to the compound symmetric and autoregressive covariance structures below.

The compound symmetric covariance has the following form: $$\mathbf{\Sigma} = \begin{bmatrix} \sigma^{2}_{1} & \sigma_{1}\sigma_{2}\rho & \sigma_{1}\sigma_{3}\rho & \dots & \sigma_{1}\sigma_{p}\rho \\ \sigma_{2}\sigma_{1} \rho & \sigma^{2}_{2} & \sigma_{2}\sigma_{3}\rho & \dots & \sigma_{2}\sigma_{p}\rho \\ \sigma_{3}\sigma_{1} \rho & \sigma_{3}\sigma_{2}\rho & \sigma^{2}_{3} & \dots & \sigma_{3}\sigma_{p}\rho \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \sigma_{p}\sigma_{1} \rho & \sigma_{p}\sigma_{2}\rho & \sigma_{p}\sigma_{3}\rho & \dots & \sigma^{2}_{p} \end{bmatrix}$$ By default, we assume a heterogeneous compound symmetric covariance structure.

Use cs(f|g) or cs(f|g, hom = TRUE) to specify a compound symmetric structure.

The AR1 (auto-regressive order 1) covariance has the following form: $$\mathbf{\Sigma} = \begin{bmatrix} \sigma^{2}_{1} & \sigma_{1}\sigma_{2} \rho & \sigma_{1}\sigma_{3}\rho^{2} & \dots & \sigma_{1}\sigma_{p}\rho^{p-1} \\ \sigma_{2}\sigma_{1} \rho & \sigma^{2}_{2} & \sigma_{2}\sigma_{3}\rho & \dots & \sigma_{2}\sigma_{p}\rho^{p-2} \\ \sigma_{3}\sigma_{1} \rho^{2} & \sigma_{3}\sigma_{2}\rho & \sigma^{2}_{3} & \dots & \sigma_{3}\sigma_{p}\rho^{p-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \sigma_{p}\sigma_{1} \rho^{p-1} & \sigma_{p}\sigma_{2}\rho^{p-2} & \sigma_{p}\sigma_{3}\rho^{p-3} & \dots & \sigma^{2}_{p} \end{bmatrix}$$ Unlike the diagonal and compound symmetric structures, by default we assume a homogeneous ar1 covariance structure.

Use ar1(f|g) or ar1(f|g, het = TRUE) to specify an AR1 structure.

Construction of the Relative Co-factor

For the unstructured covariance matrix, lme4 estimates the following parameters in par: $\mathbf{\theta} = (\theta_{1}, \theta_{2}, \dots, \theta_{p(p+1)/2})$ to construct the relative co-factor ${\mathbf \Lambda}_{\mathbf{\theta}}$ (this procedure is the same as in pre-2.0 versions of lme4): $$\Lambda_{\mathbf{\theta}} = \begin{bmatrix} \theta_1 & 0 & 0 & \dots & 0 \\ \theta_2 & \theta_3 & 0 & \dots & 0 \\ \theta_4 & \theta_5 & \theta_6 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \theta_{...} & \theta_{...} & \theta_{...} & \dots & \theta_{p(p+1)/2} \end{bmatrix}$$

The definition of the parameter vectors differs for the other covariance structures. In the diagonal covariance matrix case, $\mathbf{\theta}$ (or par) only contains the standard deviations. The relative co-factor ${\mathbf \Lambda}_{\mathbf{\theta}}$ is: $${\mathbf \Lambda}_{\mathbf{\theta}} = \begin{bmatrix} \theta_1 & 0 & 0 & \dots & 0 \\ 0 & \theta_2 & 0 & \dots & 0 \\ 0 & 0 & \theta_3 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \theta_{p} \end{bmatrix}$$

For the compound symmetry covariance structure, the parameter vector par contains the $p$ standard deviations $(\sigma_1, \sigma_2, \ldots, \sigma_p)$ and the common correlation $\rho$. In contrast to glmmTMB, the correlation is estimated on its original scale (bounded between -1 and 1), rather than on an unconstrained, transformed scale.

The relative co-factor ${\mathbf \Lambda}_{\mathbf{\theta}}$ is a lower triangular $p \times p$ matrix. Consider the form: $${\mathbf \Lambda}_{\mathbf{\theta}} = \begin{bmatrix} \theta_{11} & 0 & 0 & \cdots & 0 \\ \theta_{21} & \theta_{22} & 0 & \cdots & 0 \\ \theta_{31} & \theta_{32} & \theta_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \theta_{p1} & \theta_{p2} & \theta_{p3} & \cdots & \theta_{pp} \end{bmatrix}$$

Its elements $\theta_{ij}$ are constructed as follows. First, define the sequence $\{a_j\}$ recursively: $$a_1 = 0, \quad a_{j+1} = a_j + \frac{(1 - \rho \cdot a_j)^2}{1 - \rho^2 \cdot a_j} \quad \text{for } j = 1, \ldots, p-1$$

Then the elements of ${\mathbf \Lambda}$ are given by: $$\theta_{ij} = \begin{cases} \sqrt{1 - \rho^2 a_j} & \text{if } i = j \text{ (diagonal)} \\[0.5em] \dfrac{\rho - \rho^2 a_j}{\sqrt{1 - \rho^2 a_j}} & \text{if } i > j \text{ (below diagonal)} \\[0.5em] 0 & \text{if } i < j \text{ (above diagonal)} \end{cases}$$ Users can extract the values $\theta_{ij}$ via the getTheta() function of the Covariance.cs object, in which theta will be a vector in the column-wise elements of ${\mathbf \Lambda}$.

The setup is similar for the autoregressive order 1 (AR1) covariance structure. Again, the parameter vector contains the $p$ standard deviations $(\sigma_1, \sigma_2, \ldots, \sigma_p)$ and the autocorrelation parameter $\rho$. The relative co-factor ${\mathbf \Lambda}_{\mathbf{\theta}}$ is a lower triangular $p \times p$ matrix whose form is similar to the compound symmetric case.

The elements $\theta_{ij}$ are given by: $$\theta_{ij} = \begin{cases} \rho^{i-1} & \text{if } j = 1 \text{ (first column)} \\[0.5em] \rho^{i-j} \sqrt{1 - \rho^2} & \text{if } 1 < j \leq i \text{ (below diagonal)} \\[0.5em] 0 & \text{if } i < j \text{ (above diagonal)} \end{cases}$$ Again, these values can be extracted using getTheta() function of the Covariance.ar1 object, in which theta will be still be a vector in the column-wise elements of ${\mathbf \Lambda}$.

Extracting model components

This section illustrates how to extract par, theta, Lambda, as described in the previous section, as well as the variance covariance matrices of a model, from a merMod object.

We’ll fit the standard sleepstudy example, except that we will use a model with a compound symmetric covariance structure. Because this model has only two varying effects (intercept and slope with respect to day) per subject, and hence the covariance matrix is $2 \times 2$, there is no difference in overall model fit between the compound-symmetric and the unstructured covariance matrices. However, the models are parameterized differently, so this example will highlight the differences between par and theta.

library(lme4)
#> Loading required package: Matrix
fm1.cs <- lmer(Reaction ~ Days + cs(1 + Days | Subject), sleepstudy)

Extracting the covariance structure:

print(fm1.cs_cov <- getReCovs(fm1.cs))
#> [[1]]
#> An object of class "Covariance.cs"
#> Slot "hom":
#> [1] FALSE
#> 
#> Slot "nc":
#> [1] 2
#> 
#> Slot "par":
#> [1] 0.96679232 0.23140420 0.06561148

The result is a list with only one element as we only have one random-effects term (cs(1 + Days | Subject)). To see the values of par and theta for this object, we can call:

getME(fm1.cs, "par")
#> Subject.(Intercept)        Subject.Days         Subject.rho 
#>          0.96679232          0.23140420          0.06561148
getME(fm1.cs, "theta")
#>      Subject.(Intercept) Subject.Days.(Intercept)             Subject.Days 
#>               0.96679232               0.01518277               0.23090558

The ${\mathbf \Lambda}$ matrix is large, so we’ll view it instead of printing:

library(Matrix)
image(getME(fm1.cs, "Lambda"))

To most users, the most crucial information is simply the variance-covariance matrices. Extract these via VarCorr.merMod() (the list has one element per random-effect term in the model — in this case, only one):

vc_mat <- VarCorr(fm1.cs)
vc_mat$Subject
#>             (Intercept)      Days
#> (Intercept)  612.160273  9.613533
#> Days           9.613533 35.070443
#> attr(,"class")
#> [1] "vcmat_cs" "matrix"   "array"   
#> attr(,"stddev")
#> (Intercept)        Days 
#>    24.74187     5.92203 
#> attr(,"correlation")
#>             (Intercept)       Days
#> (Intercept)  1.00000000 0.06561148
#> Days         0.06561148 1.00000000

For detailed comparisons with glmmTMB, see browseURL(system.file("extra_docs", "covariance_structure_comparison.html", package = "lme4")).

References

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.