1
$\begingroup$

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

```
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – rvl Oct 10 at 18:45
  • $\begingroup$ Oh, and you have to add sigma too, for the total SD of the random effects. $\endgroup$ – rvl Oct 10 at 18:48
  • $\begingroup$ 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. $\endgroup$ – JWalker Oct 11 at 23:52
  • $\begingroup$ Also, I had to update my emmeans to get the bias.adj=TRUE to work. $\endgroup$ – JWalker Oct 12 at 12:41
2
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Excellent. Thanks for the time posting this. $\endgroup$ – JWalker Oct 10 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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