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"
)

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:

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 extact the posterior draws for combination of values 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 also plot the probability as well.

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