I am trying to understand how to analyze data for a generalized mixed model (GLMM) with a binary response. I am interested in visualizing the predicted probabilities, as well as a measure of effect size (e.g., odds ratios). Synthetic data are presented below. (I will also note that my actual data are skewed and have very few observations for Size >100.)


  1. Is the visualization of estimated marginal means appropriate (see below) given the structure of my data and model (i.e., ID as a random effect)? I am not interested in the conditional within-group effects, but I want to account for the repetition of ID in the data. I don't know of a way to account for ID with emmeans (or using predict(), etc.) when ID is a random effect.

  2. To obtain odds ratios for a few contrasts, I am a bit confused about the interpretation, particularly because I have used log(Size) and there are various transformations involved. Do the contrasts presented below compare odds ratios for Size 50, 100, and 250, or log(50), log(100), and log(250)? What is a simple, plain-language interpretation?

  3. I think the asymptotic CIs obtained from emmeans and confint are Wald CIs. Based on some reading, it seems like many statisticians prefer profile or bootstrapped CIs. Is there a simple way to do that for both the predicted probabilities and ORs?

The reprex:

Creating the synthetic data
# Set seed

# Create a unique ID list
unique_ids <- 1:500

# Initialize empty lists to store data
ID_list <- c()
Size_list <- c()
Practice_list <- c()
Use_list <- c()

# Generate synthetic data
for (id in unique_ids) {
  size <- runif(1, 0.5, 300)  # random Size for each ID
  practices <- c("A", "B", "C", "D", "E")  # 5 practices
  use <- sample(0:1, 5, replace = TRUE)  # random Use values for each practice

  for (practice in practices) {
    ID_list <- c(ID_list, id)
    Size_list <- c(Size_list, size)
    Practice_list <- c(Practice_list, practice)
    Use_list <- c(Use_list, use[practice == practices])

# Create data frame
synth_data <- data.frame(
  ID = ID_list,
  Size = Size_list,
  Practice = factor(Practice_list, levels = c("A", "B", "C", "D", "E")),
  Use = Use_list

synth_data$Practice <- as.factor(synth_data$Practice) # convert Practice to a factor
1. GLMM and predicted probabilities
use_mod <- glmer(Use ~ Practice * log(Size) + (1|ID), data = synth_data, family = binomial(link = "logit"), glmerControl(optimizer = "bobyqa"))

I logged Size to deal with the convergence/rescaling warnings, which is also a problem with the actual data. I will also note that I want to use the interaction because it best represents the hypotheses I am testing (i.e., whether or not it's significant is somewhat irrelevant here).

Estimated marginal means and visualization of predicted probabilities
# Specify the range of values for the continuous variable (Size)
size_values <- seq(min(synth_data$Size), max(synth_data$Size), by = 1) 

# Calculate estimated marginal means and save as data frame
use_emm <- emmeans(use_mod, ~ Practice * log(Size), at = list(Size = size_values), type = "response")
use_emm_df <- as.data.frame(use_emm)

# Visualize predicted probabilities with asymptotic 95% CIs
ggplot(use_emm_df, aes(x = log(Size), y = prob, fill = Practice)) +
  geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL), alpha = 0.2) +
  geom_line(linewidth = 1) +
  facet_grid(~ Practice) +
  guides(fill = FALSE)
2. Odds ratios
# Include key Sizes for the custom contrasts
use_emm_v2 <- emmeans(use_mod, ~ Practice * log(Size), at = list(Size = c(50, 100, 250)))

# Extract odds ratios (by Practice) and asymptotic 95% CIs and save as data frame
use_emm_or <- contrast(use_emm_v2, "pairwise", by = "Practice", type = "response") 
use_emm_or_ci <- confint(use_emm_or, adjust = "tukey")
use_emm_or_ci_df <- as.data.frame(use_emm_or_ci)

# Visualize ORs with asymptotic 95% CIs
ggplot(use_emm_or_ci_df, aes(x = contrast, y = odds.ratio, colour = contrast)) +
  geom_point(position = position_dodge(width = 0.2), size = 3) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey30") +
  facet_grid(~ Practice) +
  guides(colour = FALSE)

(By changing the contrast to by = "Size" I can also obtain comparisons of ORs by Size, rather than by Practice.)

  • $\begingroup$ Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. $\endgroup$
    – Community Bot
    Oct 17, 2023 at 0:45

1 Answer 1


Some quick partial answers...

  1. The statistics, estimates, SEs, etc. produced by emmeans are based on the model. If the model appropriately accounts for ID effects, then so do the emmeans() results. That said, with a mixed model, it might be advisable to do a bias adjustment in the back-transformed results. See the vignettes for examples/discussion; for instance, this example.

  2. It will use the Size values you specified. It does this by creating a model matrix for model predictions with the predictor combinations you specify, based on the model formula. The square-rooting is done when that model matrix is constructed.

  3. Yes, these are Wald tests/intervals. One way to do a kind of parametric bootstrap would be to construct

bemm <- regrid(emm, transform = "response", N.sim = 2000)
boddr <- regrid(pairs(emm, by = "Practice"), transform = "response", N.sim = 2000)

This creates a sample from the estimated sampling distribution of the estimated probabilities or odds ratios. (again, you might want to add sigma and bias.adj = TRUE). Anything else you do subsequently with emmeans functions treats these like MCMC results; for example, summary(bemm) gives you HPD intervals. Or you can use coda::as.mcmc() to obtain a matrix, each column being the sample you obtained.

  • $\begingroup$ Thanks a lot for your helpful response. Relatedly, would adding an adjustment for simultaneous CIs be a reasonable way to adjust the Wald CIs? For example, emmeans(use_mod, ~ Practice * log(Size), at = list(Size = size_values), type = "response", adjust = "mvt"). As I set a seq from min to max for the continuous log(Size) predictor, I get the message "mvt method for 465 estimates." I'm not actually assessing contrasts at all of those log(Size) points, but simply trying to show non-Wald CIs for my figures. Is this reasonable? The geom_ribbons are larger with this approach, as expected. $\endgroup$
    – user398696
    Nov 28, 2023 at 18:32
  • $\begingroup$ @user398696 I think this is OK. The mvt method takes into account all the correlations among the estimates, and amounts to the Studentized maximum modulus. One check on this is if you halve or double the number of means requested, the ribbons shouldn't change very much since you still have a pretty dense sequence of estimates. This approach does require a lot of computation due to the randomized algorithm used to obtain critical values. One note: the specs in emmeans only pay attention to the variables there, so you'll get the same results with Size there in place of log(Size) $\endgroup$
    – Russ Lenth
    Dec 6, 2023 at 17:26

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.