How to Get Covariance Matrix for Random Effects (Blups/Conditional Modes) from Lme4

How to get covariance matrix for random effects (BLUPs/conditional modes) from lme4

From ?ranef:

If ‘condVar’ is ‘TRUE’ each of the data frames has an attribute
called ‘"postVar"’ which is a three-dimensional array with symmetric
faces; each face contains the variance-covariance matrix for a
particular level of the grouping factor. (The name of this attribute
is a historical artifact, and may be changed to ‘condVar’ at some
point in the future.)

Set up an example:

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
rr <- ranef(fm1,condVar=TRUE)

Get the variance-covariance matrix among the b values for the intercept

pv <- attr(rr[[1]],"postVar")
str(pv)
##num [1:2, 1:2, 1:18] 145.71 -21.44 -21.44 5.31 145.71 ...

So this is a 2x2x18 array; each slice is the variance-covariance matrix among the conditional intercept and slope for a particular subject (by definition, the intercepts and slopes for each subject are independent of the intercepts and slopes for all other subjects).

To convert this to a variance-covariance matrix (see getMethod("image",sig="dgTMatrix") ...)

library(Matrix)
vc <- bdiag(  ## make a block-diagonal matrix
            lapply(
                ## split 3d array into a list of sub-matrices
                split(pv,slice.index(pv,3)),
                   ## ... put them back into 2x2 matrices
                      matrix,2)) 
image(vc,sub="",xlab="",ylab="",useRaster=TRUE)

Extract raw model matrix of random effects from lmer objects (lme4, R)

If mmList is there, then it's not going away (however poorly documented it may be -- feel free to suggest documentation improvements ...). How about

do.call(cbind,getME(m2,"mmList"))

(which would seem to generalize correctly for multi-term models)?

I agree that it's a bit of a pain that Zt doesn't distinguish correctly between structural and non-structural zeroes -- it might be possible to change the underlying code to make this work if it were sufficiently important, but I think it would be hard enough that we'd need a pretty compelling use case ...

Extract the confidence intervals of lmer random effects; plotted with dotplot(ranef())

rr <- ranef(fm1)  ## condVar = TRUE has been the default for a while

With as.data.frame: gives the conditional mode and SD, from which you can calculate the intervals (technically, these are not "confidence intervals" because the values of the BLUPs/conditional modes are not parameters ...)

dd <- as.data.frame(rr)
transform(dd, lwr = condval - 1.96*condsd, upr = condval + 1.96*condsd)

Or with broom.mixed::tidy:

broom.mixed::tidy(m1, effects = "ran_vals", conf.int = TRUE)

broom.mixed::tidy() uses as.data.frame.ranef.mer() (the method called by as.data.frame) internally: this function takes the rather complicated data structure described in ?lme4::ranef and extracts the conditional modes and standard deviations in a more user-friendly format:

If ‘condVar’ is ‘TRUE’ the ‘"postVar"’
attribute is an array of dimension j by j by k (or a list of such
arrays). The kth face of this array is a positive definite
symmetric j by j matrix. If there is only one grouping factor in
the model the variance-covariance matrix for the entire random
effects vector, conditional on the estimates of the model
parameters and on the data, will be block diagonal; this j by j
matrix is the kth diagonal block. With multiple grouping factors
the faces of the ‘"postVar"’ attributes are still the diagonal
blocks of this conditional variance-covariance matrix but the
matrix itself is no longer block diagonal.

In this particular case, here's what you need to do to replicate the condsd column of as.data.frame():

## get the 'postVar' attribute of the first (and only) RE term
aa <- attr(rr$Subject, "postVar")
## for each slice of the array, extract the diagonal;
##  transpose and drop dimensions;
##  take the square root
sqrt(c(t(apply(aa, 3, diag))))

How does one compute the normalized model residuals based via lme4/merMod in R?

nlme::lme allows for modeling correlation and heteroscedasticity in the residual error term, while lme4::lmer doesn't (these are called "R-side" structures in SAS terminology). Specifying type = "normalized" provides residuals that account for/correct for any modeled structure in the residuals; since lme4::lmer doesn't have those structures, normalizing the residuals wouldn't do anything.

If a linear mixed model (with R-side structures) is expressed as

y ~ MVN(mu, Sigma_r(phi))
mu = X beta + Z b
b ~ MVN(0, Sigma_g(theta))

where X is the fixed-effect model matrix, beta is the FE parameter vector, Z is the random-effect model matrix, b is the vector of BLUPs/conditional modes, Sigma_r and Sigma_g are the covariance matrices of the residual error and the random effects, and phi and theta are the parameter vectors defining these matrices ...

... then Sigma_r is the error correlation (covariance) matrix. In lmer, Sigma_r is always a homogeneous diagonal matrix (phi^2*I).

Variance-covariance matrix for random effect from GAM with mgcv package

library(lme4)
data(sleepstudy)

fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
fm1.rr <- ranef(fm1,condVar=TRUE)$Subject[,1]
fm1.pv <- sqrt(attr(ranef(fm1,condVar=TRUE) [['Subject']],"postVar")[1,1,])

library(mgcv)
fm2 <- gam(Reaction ~ Days + s (Subject, bs="re"), 
data = sleepstudy,   method = "REML")

To extract random effect for each Subject

idx <-grep("Subject", names(coef(fm2)))
fm2.rr<-coef(fm2)[idx]
attributes(fm2.rr)<-NULL

We can see that random effects in both models are identical as expected.

To extract variance-covariance matrix for random effect and calculate an error we use parameter Vp which is a Bayesian posterior covariance matrix:

fm2.pv <-sqrt(diag(fm2$Vp))[idx]

Or frequentist estimated covariance matrix Ve

fm2.pv <-sqrt(diag(fm2$Ve))[idx]

We can see that random effect errors estimated with mgcv slightly differ that those estimated with lme4 model. Errors based on a Bayesian posterior covariance matrix are larger, whereas based on a frequentist matrix are smaller.

---

Related Topics