Brms conditional_effects - how to plot multiple credible interval probabilities and proportional raw points?

Modeling
1

Hi everyone,

I’m currently trying to generate nice plots from my brms model using conditional_effects() and I run into two issues that I would like you to ask for advice for.

  1. I can pick a credible interval probability with the command “prob=0.95”, but I would like to plot several credible interval probabilities such as here: https://www.andrewheiss.com/blog/2021/11/10/ame-bayes-re-guide/index_files/figure-html/civlib-effect-plot-1.png How should I go about this?

  2. My dependent variable is proportional and thus looks like this: successes|trials(number of trials). The number of trials is variable. Now, when I would like to plot raw points by setting points=TRUE, I get plotted the number of successes, which are values way outside of the range of predicted probabilities (which range only from 0 to 1). How can I plot the raw probabilities or raw posterior points rather than the raw data of successes using the points=TRUE argument?

Thanks for your help!
Joe

2

Welcome to the forum @JoeL1!

To make it easier for us to answer your questions, can you please provide a reproducible example (sample data plus the code you ran)?

3 
library(tidyverse)
library(brms)

inv_logit <- function(x) {
    1 / (1 + exp(-x))
}

# Fake data simulation
n <- 100

data <- tibble(
    f = factor(rbinom(n, 1, 0.5)),
    x = rnorm(n),
    b = rnorm(n)
)

model_matrix <- model.matrix(~ 0 + f + f:x + b, data)

pars <- c(
    "f0" = 1.0,
    "f1" = 2.0,
    "b" = 5.0,
    "f0:x" = 1.0,
    "f1:x" = -1.0
)

model_data <- cbind(data, prob = inv_logit(model_matrix %*% pars)) %>%
    mutate(trials = rpois(n, 10)) %>%
    mutate(y = rbinom(n, size = trials, prob = prob))

fit <- brm(
    y | trials(trials) ~ 0 + f + f:x + b,
    model_data,
    family = "binomial"
)

cond_effect <- conditional_effects(
    fit,
    effects = c("x"),
    conditions = tibble(f = factor(c(0, 1)))
    prob=0.95
)

plot(cond_effect, plot = FALSE, points=TRUE)[[1]] +
    ylab("Probability of Success") +
    xlab("x") +
    theme_classic()
4

Here’s some code that I hope will answer your questions. First some set-up, based on your reproducible example:

library(tidyverse)
library(brms)
library(ggdist)
theme_set(theme_bw())

# FYI: brms has built-in inv_logit_scaled() function
inv_logit <- function(x) {
    1 / (1 + exp(-x))
}

# Fake data simulation
n <- 100

set.seed(5847)
data <- tibble(
    f = factor(rbinom(n, 1, 0.5)),
    x = rnorm(n),
    b = rnorm(n)
)

model_matrix <- model.matrix(~ 0 + f + f:x + b, data)

pars <- c(
    "f0" = 1.0,
    "f1" = 2.0,
    "b" = 5.0,
    "f0:x" = 1.0,
    "f1:x" = -1.0
)

model_data <- cbind(data, prob = inv_logit(model_matrix %*% pars)) %>%
    mutate(trials = rpois(n, 10)) %>%
    mutate(y = rbinom(n, size = trials, prob = prob))

fit <- brm(
    y | trials(trials) ~ 0 + f + f:x + b,
    model_data,
    family = "binomial",
    file="fit"
)

# Note that I've used the make_conditions() function to add informative 
#  labels to the facets for f
cond_effect <- conditional_effects(
    fit,
    effects = c("x"),
    conditions = make_conditions(fit, "f"),
    prob=0.95
)
#> Setting all 'trials' variables to 1 by default if not specified otherwise.

plot(cond_effect, plot = FALSE, points=TRUE)[[1]] +
    ylab("Probability of Success") +
    xlab("x")

conditional_effects() returns a list of data frames for plotting. If the focal variable is continous (like x) it returns a data frame with 100 rows for each condition, in this case f=0 and f=1. It sets other variables to their mean if continuous or reference level if factor. conditional_effects() also sets the number of trials to 1 by default (as reported in the message returned to the console).

glimpse(cond_effect$x)
#> Rows: 200
#> Columns: 11
#> $ x          <dbl> -1.9648819, -1.9202264, -1.8755708, -1.8309152, -1.7862596,…
#> $ y          <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,…
#> $ trials     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ f          <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
#> $ b          <dbl> 0.0317888, 0.0317888, 0.0317888, 0.0317888, 0.0317888, 0.03…
#> $ cond__     <fct> f = 0, f = 0, f = 0, f = 0, f = 0, f = 0, f = 0, f = 0, f =…
#> $ effect1__  <dbl> -1.9648819, -1.9202264, -1.8755708, -1.8309152, -1.7862596,…
#> $ estimate__ <dbl> 0.2989807, 0.3089284, 0.3188315, 0.3290534, 0.3393385, 0.35…
#> $ se__       <dbl> 0.10252110, 0.10229261, 0.10216625, 0.10196309, 0.10156782,…
#> $ lower__    <dbl> 0.1371737, 0.1450188, 0.1533625, 0.1623946, 0.1717314, 0.18…
#> $ upper__    <dbl> 0.5369587, 0.5427293, 0.5483094, 0.5544026, 0.5606309, 0.56…

We can create customized plots by creating a data frame similar to conditional effects. In this case, I’ll create a data frame of prediction values for a few different numbers of trials (you can explore other combinations of values by updating this data frame as appropriate):

pred.dat = crossing(x=seq(min(model_data$x), max(model_data$x), length=100),
                    f = c(0,1),
                    trials = c(1,5,10), 
                    b = mean(model_data$b))

Then we extract the posterior draws for each combination of values (each row) in pred.dat:

post.epred = epred_draws(fit, newdata=pred.dat, ndraws=500) 
post.epred
#> # A tibble: 300,000 × 9
#> # Groups:   x, f, trials, b, .row [600]
#>        x     f trials      b  .row .chain .iteration .draw .epred
#>    <dbl> <dbl>  <dbl>  <dbl> <int>  <int>      <int> <int>  <dbl>
#>  1 -1.96     0      1 0.0318     1     NA         NA     1  0.221
#>  2 -1.96     0      1 0.0318     1     NA         NA     2  0.409
#>  3 -1.96     0      1 0.0318     1     NA         NA     3  0.262
#>  4 -1.96     0      1 0.0318     1     NA         NA     4  0.252
#>  5 -1.96     0      1 0.0318     1     NA         NA     5  0.212
#>  6 -1.96     0      1 0.0318     1     NA         NA     6  0.169
#>  7 -1.96     0      1 0.0318     1     NA         NA     7  0.203
#>  8 -1.96     0      1 0.0318     1     NA         NA     8  0.261
#>  9 -1.96     0      1 0.0318     1     NA         NA     9  0.303
#> 10 -1.96     0      1 0.0318     1     NA         NA    10  0.515
#> # ℹ 299,990 more rows

Now we can create a plot similar to the plot you linked to from Andrew Heiss’s blog:

ggplot(post.epred, aes(x = x, y = .epred)) +
  stat_lineribbon(.width=c(0.5, 0.8, 0.95), linewidth=0.3) + 
  geom_point(data=model_data %>% 
               filter(trials %in% unique(pred.dat$trials)), 
             aes(y=y), colour="blue", size=1, alpha=0.5) +
  scale_fill_brewer(palette = "Reds") +
  theme(legend.position = "bottom") +
  facet_grid(rows=vars(trials), cols=vars(f), scales="free_y", 
             labeller=label_both) +
  scale_y_continuous(breaks=1:20) +
  labs(y="Number of successes")

Note that your model is predicting the number of successes, rather than the probability of success. But we can plot the probability as well. Also, the data used to fit the model doesn’t have any observations for cases where the number of trials equals 1, so there are no points plotted in the trials=1 facets.

ggplot(post.epred %>% mutate(.epred = .epred/trials), 
       aes(x = x, y = .epred)) +
  stat_lineribbon(.width=c(0.5, 0.8, 0.95), linewidth=0.3) + 
  geom_point(data=model_data %>% 
               filter(trials %in% unique(pred.dat$trials)) %>% 
               mutate(y = y/trials), 
             aes(y=y), colour="blue", size=1, alpha=0.5) +
  scale_fill_brewer(palette = "Reds") +
  theme(legend.position = "bottom") +
  facet_grid(rows=vars(trials), cols=vars(f), labeller=label_both) +
  labs(y="Probability of success")

Created on 2025-06-14 with reprex v2.1.1

1 Like