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!


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

# 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 
  tibble(Cohort = "Cohort1", Subject = "308", Days = c(1, 2))
) %>%
#>        2 
#> 18.70102

Created on 2023-03-26 with reprex v2.0.2


1 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 = "+")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.