--- title: "Covariance Structures" output: rmarkdown::html_vignette bibliography: lmer.bib vignette: > %\VignetteIndexEntry{Covariance Structures} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} %\VignetteDepends{glmmTMB} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## 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](https://glmmtmb.github.io/glmmTMB/articles/covstruct.html)), 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](https://cran.r-project.org/package=lme4/vignettes/lmer.pdf) [@Bates_JSS]. 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`; $\sigma^2_2, \dots $\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`. ```{r cs-fit} library(lme4) fm1.cs <- lmer(Reaction ~ Days + cs(1 + Days | Subject), sleepstudy) ``` Extracting the covariance structure: ```{r getReCovs} print(fm1.cs_cov <- getReCovs(fm1.cs)) ``` 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: ```{r getInfo} getME(fm1.cs, "par") getME(fm1.cs, "theta") ``` The ${\mathbf \Lambda}$ matrix is large, so we'll view it instead of printing: ```{r lambda-mat} 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): ```{r varcorr} vc_mat <- VarCorr(fm1.cs) vc_mat$Subject ``` For detailed comparisons with `glmmTMB`, see `browseURL(system.file("extra_docs", "covariance_structure_comparison.html", package = "lme4"))`. ## References