Extracting the (i) estimated variance-covariance matrix of random effects and/or (ii) mixed model equation solutions from lme4 as matrices?

0
votes

As the title says I am trying to extract matrices from an lme4 (or other packages?) object. To make clear what I want precisely I think it is easiest to refer to the SAS documentation: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect022.htm

Variance-covariance matrix of random effects

In SAS notation this matrix is called G and is the variance-covariance matrix of the random effect parameter gamma. By using the option "G" in PROC MIXED and the Output Delivery System you obtain G as a matrix.

I am aware that it is relatively simple to construct this matrix manually once I have the variance components and dimensions of gamma. I nevertheless expected there to be an even simpler way.

Mixed model equations solution

In SAS notation these are called C. By using the option "MMEQSOL" in PROC MIXED and the Output Delivery System you request that a solution to the mixed model equations be produced, as well as the inverted coefficients matrix. It is the latter that I am interested in.

Thanks in advance!

1
rsaslme4mixed-modelsnlme

1 Answers

3
votes

Not a very sensible model (see ?lme4::cake), but reasonable for illustration:

library(lme4)
fm1 <- lmer(angle ~ temperature + 
           (1|recipe)+(1|replicate), cake)

The VarCorr() method gives a list of variance-covariance matrices for each term (in this case each one is 1x1), with its own print method:

v <- VarCorr(fm1)

You can combine these into a single matrix by using the bdiag() (block-diagonal) function from Matrix (as.matrix() converts from a sparse matrix to a standard (dense) R matrix object).

as.matrix(Matrix::bdiag(v))
##          [,1]      [,2]
## [1,] 39.21541 0.0000000
## [2,]  0.00000 0.4949681

The C matrix is unfortunately not so easy to get. As discussed in vignette("lmer",package="lme4"), lme4 doesn't use the Henderson-equation formulation. The upper block of C (variance-covariance matrix of fixed effects) is accessible via vcov(), but the variance-covariance matrix of the variances is not so easy: see e.g. here.