How can I extract and report the random slopes and confidence from a glmer in R?

I would like to report the random slopes from a binomial lme4::glmer model along with their confidence or deviations. I am adding the fixed effect to each random effect to obtain slopes, but how do I combine condsd from ranef and 2.5 %, 97.5 % from confint?

mod <- lme4::glmer(dependent~pred1+pred2+pred1+pred2|ID), data = df, family = "binomial")

fixed_effects <- confint(mod, method = "Wald") %>% 
  as.data.frame() %>% 
  rownames_to_column(var = "term") %>%
  left_join(lme4::fixef(mod) %>% 
              as.data.frame() %>% 
              rename(fixef = ".") %>% 
              rownames_to_column(var = "term"), 
            by = join_by(term)) 

random_effects <- as.data.frame(lme4::ranef(mod))

effects <- random_effects %>% 
  left_join(fixed_effects) %>% 
  mutate(rnd_slope = fixef+condval)

So far, I have found many helpful posts here and on CrossValidated to address the earlier steps of modeling and then extracting model coefficients, but I have found none that addressed handling confidence when reporting random slopes (rather than random effects).

Solution

This is possible, but currently much harder than it ought to be. Getting the joint covariance matrix of random and fixed effects in lme4 has only been possible relatively recently thanks to work by Emi Tanaka; it's used in predict.merMod(..., se.fit = TRUE) but the machinery isn't exposed in a way that's more generally useful (i.e. the user should be able to retrieve the full joint (conditional) covariance matrix, and coef.merMod() should also optionally be able to return conditional SDs ...

So here we go (key parts of the code are taken from predict.merMod()):

set up examples

library(lme4)
object <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)

derive joint conditional covariance matrix of b, beta

s <- sigma(object)
L <- getME(object, "L")
RX <- getME(object, "RX")
RZX <- getME(object, "RZX")
Lambdat <- getME(object, "Lambdat")
RXtinv <- solve(t(RX))
LinvLambdat <- solve(L, Lambdat, system = "L")
Minv <- s * rbind(
                cbind(LinvLambdat,
                      Matrix(0, nrow = nrow(L), ncol = ncol(RX))),
                cbind(-RXtinv %*% t(RZX) %*% LinvLambdat, RXtinv)
            )
Cmat <- crossprod(Minv)

helper function for efficient diagonal of cross product

quad.tdiag <- function(M, x) {
    ## only real-valued, so drop Conj
    rowSums(tcrossprod(x, M) * x)
}
## more brute-force equivalent: diag(x %*% M %*% t(x))

set up 'contrast' matrix

This is the ugliest part (finding the joint covariance matrix is more magical, but this part is ugly because I haven't thought about a general way to do it, instead hacking something that will work OK in this example). We want to set up a matrix such that each row when multiplied by a (b, beta) vector gives the slope for a particular subject. We need to know that (1) the b vector is stored as an (intercept, slope) offset pair for each subject; and (2) the beta vector (intercept, slope) is the last two elements in the vector. Thus e.g. if we take a matrix row

0 1 0 0 0 ...... 0 1

and multiply it by (b, beta) we'll get b_slope(1) + beta_slope, i.e. the slope for subject 1.

This code constructs a matrix of such rows (not necessarily in the most efficient way possible).

nsubj <- length(levels(sleepstudy$Subject))
v <- Matrix(0, nrow = nsubj, ncol = nrow(Cmat))
v[,ncol(Cmat)] <- 1
v[col(v) == 2*row(v)] <- 1

Now the quadratic form v Cmat v^T gives us the covariance matrix of the elements of v, but we only need the diagonal of this matrix. Take the square root to get SDs/SEs:

sqrt(quad.tdiag(Cmat, v))
##  [1] 2.336278 2.336278 2.336278 2.336278 2.336278 2.336278 2.336278 2.336278
##  [9] 2.336278 2.336278 2.336278 2.336278 2.336278 2.336278 2.336278 2.336278
## [17] 2.336278 2.336278

(all the same because this is a balanced design)

Compare to naive (assuming uncorrelated) answer:

fixed_slope_var <- vcov(object)[2, 2]
ranef_var <- attr(ranef(object)$Subject, "postVar")[2,2,]
sqrt(fixed_slope_var + ranef_var)

 [1] 2.775202 2.775202 2.775202 2.775202 2.775202 2.775202 2.775202 2.775202
 [9] 2.775202 2.775202 2.775202 2.775202 2.775202 2.775202 2.775202 2.775202
[17] 2.775202 2.775202