3
$\begingroup$

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!

$\endgroup$

1 Answer 1

3
$\begingroup$

Quick responses, subject to someone else pointing out my dumb oversights...

  1. Yes, I think that is a reasonable interpretation.

  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.

  3. There are a lot of moving parts here and I need to ponder this a bit more. But my instincts are that the covariance of the estimates - which is used to construct the CIs -- have you covered for population-level inferences on proportions, once the bias adjustment is performed.

$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$ Dec 5, 2020 at 12:58
  • $\begingroup$ In addition, the odds ratio reported by contrasts needs to be corrected too? $\endgroup$ Dec 5, 2020 at 13:22

Your Answer

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

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