Are the coefficients in `lme4` correct with nested effects?

Question

The coefficients I'm seeing from a merMod model produced by lmer are not what I'd expect. I made a reproducible example to illustrate.

I assumed that the Days coefficient reported for a nested effect (in this case, subject "308" within "Cohort1") would be the sum of the Days fixed effect plus the random effect coefficient for Days within Cohort and within Subject.

I can produce this coefficient myself by adding exactly this way (18.7) and confirm with predictions (18.7). But the reported coefficient is different (21.15).

Am I interpreting something wrong? Appreciate any help!

library(tidyverse)
library(lme4)

# add a categorical variable "Cohort" that Subject can be nested in
df <- sleepstudy %>%
  mutate(Cohort = if_else(
    as.numeric(as.character(Subject)) < 336,
    "Cohort1",
    "Cohort2"
  ))

# fit
sleep_model <- lmer(Reaction ~ Days + (Days | Cohort/Subject), df)
#> boundary (singular) fit: see help('isSingular')

# the report coef
coef(sleep_model)[["Subject:Cohort"]]["308:Cohort1", "Days"]
#> [1] 21.15208

# what I think the coef would be
fixef(sleep_model)["Days"] +
  ranef(sleep_model)[["Cohort"]]["Cohort1", "Days"] +
  ranef(sleep_model)[["Subject:Cohort"]]["308:Cohort1", "Days"]
#>     Days 
#> 18.70102

# predict and divide delta y / delta x 
predict(
  sleep_model, 
  tibble(Cohort = "Cohort1", Subject = "308", Days = c(1, 2))
) %>%
  diff()
#>        2 
#> 18.70102

^{Created on 2023-03-26 with reprex v2.0.2}

Answer 1

coef() gives you the adjusted coefficients (fixed + random effects) for each random-effect term separately; the key thing to know here is that Cohort/Subject is counted as two separate random effect terms by the model. That is, (Days | Cohort/Subject) is internally expanded to (Days | Cohort) + (Days | Cohort:Subject), and dealt with thereafter (for ranef()/coef() purposes) as two separate terms, not as a combined term. So what you are getting when you ask for coef(sleep_model)[["Subject:Cohort"]]["308:Cohort1", "Days"] is indeed just the adjustment for the subject-by-cohort combination, not the combination of the cohort and subject effects.

I realize that may be a nuisance. It would be possible to augment/expand the coef method so that it collapsed these terms together (although a bit of a pain, as the fitted model object doesn't keep track of which RE terms were expanded from nested components in the original formula ...)

I think this should do what you want, provided the grouping variable names/levels within them are not pathologically set up so that grep() gets weird results.

object <- sleep_model
## nested terms to evaluate
nests <- c("Cohort", "Subject:Cohort")
ref <- ranef(object)
maxnest <- ref[nests][[length(nests)]] ## innermost nested RE
## start with FE, shape matching innermost RE
vals <- lapply(fixef(object), rep, nrow(maxnest))
vals <- data.frame(vals, check.names = FALSE)
rownames(vals) <- rownames(maxnest)
## 
for (n in nests) {
    for (nm in rownames(ref[[n]])) {
        ## add relevant terms to running total (yikes)
        matches <- grep(nm, rownames(vals))
        vals[matches,] <- sweep(vals[matches,], 
                 MARGIN = 2, ## columnwise
                 STATS = unlist(ref[[n]][nm,]), ## cond modes for this group
                 FUN = "+")
    }
}