emmeans output interpretation of a glmer fit with nesting

I have read that the interpretation of generalized linear mixed models (GLMM) at the response level is more complex because the back transformation is nonlinear and the random terms do not play a strictly additive role. In this sense, I would like to know what would be the interpretation of the emmeans result of a glmer fit.

I simulated a data table with two fixed and two random variables (one is nested)

  rm(list=ls(all=TRUE))
set.seed(1)
df1 <- data.frame(period = rep(c(-1,0,0.5,1),2), treat = rep(c(0,1),each = 4))
df2 <- data.frame(id = round(rnorm(20),2), id2 = rep(round(rnorm(4),2),each = 5))
df <- merge(df1,df2)
df$$link <- rowSums(df); df$$size = 10
df$$pos <- rbinom(160,10,exp(df$$link)/(1+exp(df$$link))) df$$period <- factor(df$$period); df$$treat = factor(df$treat) str(df) mod <- glmer(cbind(pos, size - pos) ~ period + treat + (1 | id2/id), data = df, family = binomial) fixef(mod) (Intercept) period0 period0.5 period1 treat1 -0.8248334 0.9193305 1.5729769 2.0988980 0.9039252 exp(fixef(mod)[5]) treat1 2.469277 #Descriptive proportions agg <-aggregate(cbind(pos, size) ~ period+treat, sum, data = df) agg$$prop <- agg$$pos/agg$size;agg

period treat pos size  prop
1     -1     0  79  200 0.395
2      0     0 103  200 0.515
3    0.5     0 124  200 0.620
4      1     0 142  200 0.710
5     -1     1  98  200 0.490
6      0     1 137  200 0.685
7    0.5     1 157  200 0.785
8      1     1 168  200 0.840

#Estimated probability
emm <- emmeans(mod, ~ period+treat, type = "response");emm

period treat  prob     SE  df asymp.LCL asymp.UCL
-1     0     0.305 0.1104 Inf     0.136     0.549
0      0     0.524 0.1299 Inf     0.284     0.753
0.5    0     0.679 0.1139 Inf     0.432     0.855
1      0     0.781 0.0897 Inf     0.561     0.909
-1     1     0.520 0.1299 Inf     0.281     0.750
0      1     0.731 0.1027 Inf     0.494     0.883
0.5    1     0.839 0.0709 Inf     0.651     0.936
1      1     0.898 0.0484 Inf     0.758     0.961

Confidence level used: 0.95
Intervals are back-transformed from the logit scale


Included after Prof. Lenth response

  var.mod <- as.data.frame(lme4::VarCorr(mod))[,4]
total.SD = sqrt(sum(var.mod)) ; total.SD
[1] 1.371339

#Bias corrected estimation
summ <- summary(emm, type = "response", bias.adjust = TRUE, sigma = total.SD) ; summ
period treat  prob     SE  df asymp.LCL asymp.UCL
-1     0     0.383 0.0823 Inf     0.217     0.526
0      0     0.513 0.0692 Inf     0.366     0.665
0.5    0     0.605 0.0809 Inf     0.463     0.772
1      0     0.691 0.0876 Inf     0.533     0.846
-1     1     0.510 0.0691 Inf     0.364     0.662
0      1     0.645 0.0853 Inf     0.497     0.809
0.5    1     0.753 0.0836 Inf     0.586     0.887
1      1     0.830 0.0690 Inf     0.669     0.929

Confidence level used: 0.95
Intervals are back-transformed from the logit scale
Bias adjustment applied based on sigma = 1.3713

#Comparing estimations
comparison = data.frame(desc = agg\$prop, biased = as.data.frame(emm)[,3],
unbiased = as.data.frame(summ)[,3]) ; comparison
desc    biased  unbiased
1 0.395 0.3047386 0.3825390
2 0.515 0.5236067 0.5125329
3 0.620 0.6787740 0.6054697
4 0.710 0.7814378 0.6910434
5 0.490 0.5197627 0.5104859
6 0.685 0.7307483 0.6453686
7 0.785 0.8391704 0.7530862
8 0.840 0.8982555 0.8298073

#Distance from estimations
dist(t(comparison))
desc     biased
biased   0.16151096
unbiased 0.06200439 0.19796862



I have three questions:

1) Regarding the fixed coefficient treat, is it correct to say that, keeping all the rest fixed including the random effects, the expected odds of pos increases 2.47 times in the treatment group (treat1)?

2) Regarding the emmeans output, can I interpret the estimated probability in the population level? For instance "With 95% confidence, the population probability of treat0 in the period -1 is within the range [0.136 ; 0.549]"?

3) Do the presence of nesting affects the interpretation, or do I need to choose at which hierarchical level I will make the interpretation (i.e. at id or id2 level)?

Thank you very much!

2. No. The estimate and CI are biased due to the issue of back-transforming on a nonlinear scale. You need todo a bias adjustment based on the total SD of the random effects. To do this, obtain the random-effects estimates via VarCorr(), and refer to the last section of the emmeans transformations vignette -- https://cran.r-project.org/web/packages/emmeans/vignettes/transformations.html#bias-adj.
• Thank you, Prof. Lenth! I edited my question including the bias correct estimation, as described in your vignette. If I need to report point estimation, can we say that the biased corrected is a better choice than the descriptive proportions calculated in agg? Dec 5, 2020 at 12:58
• In addition, the odds ratio reported by contrasts needs to be corrected too? Dec 5, 2020 at 13:22