# Should modeled means in glmm be averaged on link or response scale?

I'm using emmeans to estimate means of 4 treatment levels. The response is a count with poisson distribution. The design is randomized complete block with 14 blocks, each with single replicate of all four levels. N = 14*4=56. The model is glmm using poisson family. If the means of each treatment are averaged over the 14 blocks on the link scale and then backtransformed to the response scale, the treatment means are substantially less than the means taken on the raw scale (or the true means). So, the modeled mean doesn't predict very well the true mean. Do we want this behavior? Up to now, I've agreed with Russell Lenth (author of emmeans) that means of glm fits should be averaged on the link scale and then backtransformed. But this behavior questions this for glmm.

This code replicates the question with fake data. I include modeled means from a lm, glm, lmm, and the glmm. The final column is the means taken over the prediction on the response scale.

library(lme4)
library(emmeans)
library(data.table)
set.seed(10)
n_sites <- 14
treatments <- c("cn", "tr-a", "tr-b", "tr-c")
n_treatments <- length(treatments)
fake_data <- data.table(site=letters[1:n_sites], treatment=rep(treatments, each=n_sites))
X <- model.matrix(~treatment, data=fake_data)
sigma_i <- 1.000001
b0 <- log(4.6)
b0_i <- rnorm(n_sites, sd=sigma_i)
b <- c(b0, 0, 0, 0)
fake_data[, mu:=exp((X%*%b)[,1] + rep(b0_i, n_treatments))]
fake_data[, count:=rpois(n_sites*n_treatments, lambda=mu)]
m0 <- lm(count ~ treatment, data=fake_data)
m1 <- lmer(count ~ treatment + (1|site), data=fake_data)
m2 <- glm(count ~ treatment, family=poisson(link="log"), data=fake_data)
m3 <- glmer(count ~ treatment + (1|site), family=poisson(link="log"), data=fake_data)

fake_data[, glmm:=predict(m3, type="response")]

# get modeled means
m0_emm <- data.table(summary(emmeans(m1, specs="treatment", type="response")))
setnames(m0_emm, old="emmean", new="lm")
part1 <- m0_emm[, .SD, .SDcols=c("treatment", "lm")]
m1_emm <- data.table(summary(emmeans(m1, specs="treatment", type="response")))
setnames(m1_emm, old="emmean", new="lmm")
part2 <- m1_emm[, .SD, .SDcols=c("treatment", "lmm")]
m2_emm <- data.table(summary(emmeans(m2, specs="treatment", type="response")))
setnames(m2_emm, old="rate", new="glm")
part3 <- m2_emm[, .SD, .SDcols=c("treatment", "glm")]
m3_emm <- data.table(summary(emmeans(m3, specs="treatment", type="response")))
setnames(m3_emm, old="rate", new="glmm")
part4 <- m3_emm[, .SD, .SDcols=c("treatment", "glmm")]
m3_resp <- fake_data[, .(glmm_resp=mean(glmm)), by=treatment]
# put it all together
merge(merge(merge(merge(part1, part2, by="treatment"), part3, by="treatment"), part4, by="treatment"), m3_resp, by="treatment")

   treatment       lm      lmm      glm     glmm glmm_resp
1:        cn 5.500000 5.500000 5.500000 3.955373  5.487603
2:      tr-a 4.357143 4.357143 4.357143 3.133485  4.347332
3:      tr-b 4.785714 4.785714 4.785714 3.441701  4.774945
4:      tr-c 5.142857 5.142857 5.142857 3.698570  5.131319

$$$$

• One way to think about this is that you are really estimating the geometric means. If you want things on the expected response scale, add bias.adj = TRUE to the emmeans calls. See the vignette on transformations. – Russ Lenth Oct 10 '19 at 18:45
• Oh, and you have to add sigma too, for the total SD of the random effects. – Russ Lenth Oct 10 '19 at 18:48
• Yes. Not having had any formal course/training in glm or glmm (or even lmm), I've found it helpful to explore the behavior of the conditional vs. marginal means by changing my sigma_1, so controlling the among:within variance. The whole exercise has been really mind-expanding for me, and it helps alot to combine my exploration with your math. – JWalker Oct 11 '19 at 23:52
• Also, I had to update my emmeans to get the bias.adj=TRUE to work. – JWalker Oct 12 '19 at 12:41

What you are observing here is an example of the marginal versus conditional interpretation of the fixed effects coefficients from generalized linear mixed-effects models (GLMMs). Namely, in GLMMs the fixed effects have an interpretation conditional on the random effects.

For your particular model, and because you have random intercepts only the fixed effect intercept has an interpretation conditional on the random effects. The rest of the coefficients have the usual population-averaged interpretation. You can obtain the marginal intercept by using the transformation: $$\beta_0^M = \beta_0^C + \frac{\sigma_b^2}{2},$$ where $$\beta_0^C$$ is the intercept you obtain from fitting the GLMM, and $$\sigma_b^2$$ the variance of the random intercepts. With the new intercept the results are much closer, i.e.,

betas <- fixef(m3)
betas[1] <- betas[1] + 0.5 * VarCorr(m3)[[1]][1, 1]

X <- model.matrix(~ treatment,
data = data.frame(treatment = unique(fake_data\$treatment)))
exp(c(X %*% betas))
[1] 5.692687 4.509802 4.953396 5.323088
`
• Excellent. Thanks for the time posting this. – JWalker Oct 10 '19 at 14:52