Marginal variance-covariance matrix from an LMM

Calculates the marginal variance-covariance matrix from an lme4::lmer or nlme::lme object.

Usage

lmm_marginal_varcov(x)

Arguments

x: An lme4::lmer or nlme::lme model object.

Value

dgCMatrix

Details

For linear mixed model \(\boldsymbol{y}_i = \boldsymbol{X}_i \boldsymbol{\beta} + \boldsymbol{Z}_i \boldsymbol{b}_i + \boldsymbol{\epsilon}_i\), where \(i\) is the group index, \(\boldsymbol{y}_i\) is an \(n_i\)-dimensional vector of observed responses, \(\boldsymbol{X}_i\) and \(\boldsymbol{Z}_i\) are known \(n_i \times p\) and \(n_i \times q\) regression matrices corresponding to the \(p\)-dimensional fixed effects vector \(\boldsymbol{\beta}\) and the \(q\)-dimensional random effects vector \(b_i\), respectively, and \(\boldsymbol{\epsilon}_i\) is an \(n_i\)-dimensional vector of within-group errors.

The \(\boldsymbol{b}_i\) are assumed to be independently distributed as \(N(\boldsymbol{0}, \boldsymbol{\Psi})\) and the \(\boldsymbol{\epsilon}_i\) are assumed to be independently distributed as \(N(\boldsymbol{0}, \boldsymbol{\Lambda}_i)\), independent of the \(\boldsymbol{b}_i\). The \(\boldsymbol{\Psi}\) covariance matrix may be unstructured or structured (commonly a compound symmetric or general symmetric positive definite matrix parameterized by \(\boldsymbol{\theta}\)). The \(\boldsymbol{\Lambda}\) matrices are assumed to depend on \(i\) only through their dimensions \(\boldsymbol{\lambda}\) (commonly a diagonal, compound/general symmetric, or AR(1) matrix parameterized by \(\boldsymbol{\lambda}\)).

We are interested in jointly estimating the regression model parameters (relationship between response and covariates) and the covariance parameters (the correlation and variance structure of the response). The model approach using random effects is called the conditional model. The model approach which directly models the within-group error covariance is the called the marginal model.

Likelihood estimation is based on the marginal distribution of the observed response vectors \(\boldsymbol{y}_i\). Since the random effects and within-group error term enter the model linearly, it can be shown that the \(\boldsymbol{y}_i\) are marginally distributed as independent \(N(\boldsymbol{X}_i\boldsymbol{\beta}, \boldsymbol{\Sigma}_i)\) random vectors. The marginal covariance matrix is given by \(var(\boldsymbol{y}_i) = \boldsymbol{\Sigma}_i = \boldsymbol{Z}_i \boldsymbol{\Psi} \boldsymbol{Z}_i^{\prime} + \boldsymbol{\Lambda}_i\), where \(\boldsymbol{Z}_i \boldsymbol{\Psi} \boldsymbol{Z}_i^{\prime}\) is the random effects variance component and \(\boldsymbol{\Lambda}_i\) is the within-group error variance component.

Linear mixed effects models can account for within-group correlation and heteroscedasticity through the random effects \(\boldsymbol{b}_i\) and the within-group errors \(\boldsymbol{\epsilon}_i\). Since the \(\boldsymbol{b}_i\) are fixed by group, the within-group observations share the same random effects and are correlated. The \(\boldsymbol{Z}_i \boldsymbol{\Psi} \boldsymbol{Z}_i^{\prime}\) variance component describes this correlation. The diagonal of \(\boldsymbol{Z}_i \boldsymbol{\Psi} \boldsymbol{Z}_i^{\prime}\) can also accommodate heteroscedasticity. If fitting a conditional model, \(\boldsymbol{Z}_i \boldsymbol{\Psi} \boldsymbol{Z}_i^{\prime}\) is usually favored for describing the structure of the data and \(\boldsymbol{\Lambda}_i\) typically assumes a simple form such as \(\boldsymbol{\Lambda}_i = \sigma^2 \boldsymbol{I}_i\). If fitting a marginal model, \(\boldsymbol{\Sigma}_i = \boldsymbol{\Lambda}_i\) so it must accommodate correlation (non-diagonal elements) and heteroscedasticity.

In practice of fitting LMMs, one can use either variance component to account for correlation and heteroscedasticity. However, the two components may compete in explaining the marginal covariance and lead to overparameterized models.

References

https://stackoverflow.com/a/45655597

Pinheiro JC (2006). “Conditional versus Marginal Covariance Representation for Linear and Nonlinear Models.” Austrian Journal of Statistics, 35(1). ISSN 1026-597X, doi:10.17713/ajs.v35i1.346 .

Venables WN, Ripley BD (2002). Modern Applied Statistics with S, Statistics and Computing. Springer, New York, NY. ISBN 9781441930088, doi:10.1007/978-0-387-21706-2 .

Pinheiro JC, Bates DM (2000). Mixed-Effects Models in S and S-PLUS. Springer, New York, NY. ISBN 9781475781441, doi:10.1007/b98882 .

Examples


#----------------------------------------------------------------------------
# lmm_marginal_varcov() examples
#----------------------------------------------------------------------------
library(bkstat)

data(Orthodont, package="nlme")

mod_lmer <- lme4::lmer(
  formula = distance ~ age + Sex + (1 + age | Subject),
  data = Orthodont
)
mod_lme <- nlme::lme(
  fixed = distance ~ age + Sex,
  random = ~ 1 + age | Subject,
  #correlation = nlme::corCompSymm(form = ~ 1 | Subject),
  data = Orthodont
)

vcov_lmer <- lmm_marginal_varcov(mod_lmer)
vcov_lme <- lmm_marginal_varcov(mod_lme)
vcov_lme_check <- nlme::getVarCov(
  obj = mod_lme,
  individuals = levels(mod_lme$groups$Subject),
  type = "marginal"
)
vcov_lme_check <- Matrix::bdiag(vcov_lme_check)

Matrix::image(vcov_lmer)

Matrix::image(vcov_lme)

Matrix::image(vcov_lme_check)


vcov_lmer[1:8, 1:8]
#> 8 x 8 sparse Matrix of class "dgCMatrix"
#>          1        2        3        4        5        6        7        8
#> 1 5.060294 3.194403 3.044722 2.895042 .        .        .        .       
#> 2 3.194403 4.965992 3.305161 3.360540 .        .        .        .       
#> 3 3.044722 3.305161 5.281810 3.826039 .        .        .        .       
#> 4 2.895042 3.360540 3.826039 6.007747 .        .        .        .       
#> 5 .        .        .        .        5.060294 3.194403 3.044722 2.895042
#> 6 .        .        .        .        3.194403 4.965992 3.305161 3.360540
#> 7 .        .        .        .        3.044722 3.305161 5.281810 3.826039
#> 8 .        .        .        .        2.895042 3.360540 3.826039 6.007747
vcov_lme[1:8, 1:8]
#> 8 x 8 sparse Matrix of class "dsCMatrix"
#>                                                                             
#> [1,] 5.060332 3.194382 3.044636 2.894890 .        .        .        .       
#> [2,] 3.194382 4.965917 3.305045 3.360376 .        .        .        .       
#> [3,] 3.044636 3.305045 5.281658 3.825862 .        .        .        .       
#> [4,] 2.894890 3.360376 3.825862 6.007553 .        .        .        .       
#> [5,] .        .        .        .        5.060332 3.194382 3.044636 2.894890
#> [6,] .        .        .        .        3.194382 4.965917 3.305045 3.360376
#> [7,] .        .        .        .        3.044636 3.305045 5.281658 3.825862
#> [8,] .        .        .        .        2.894890 3.360376 3.825862 6.007553
vcov_lme_check[1:8, 1:8]
#> 8 x 8 sparse Matrix of class "dsCMatrix"
#>                                                                             
#> [1,] 5.060332 3.194382 3.044636 2.894890 .        .        .        .       
#> [2,] 3.194382 4.965917 3.305045 3.360376 .        .        .        .       
#> [3,] 3.044636 3.305045 5.281658 3.825862 .        .        .        .       
#> [4,] 2.894890 3.360376 3.825862 6.007553 .        .        .        .       
#> [5,] .        .        .        .        5.060332 3.194382 3.044636 2.894890
#> [6,] .        .        .        .        3.194382 4.965917 3.305045 3.360376
#> [7,] .        .        .        .        3.044636 3.305045 5.281658 3.825862
#> [8,] .        .        .        .        2.894890 3.360376 3.825862 6.007553