Plot random effects for intercept and slope with tidybayes

Modeling
1

I want to plot the random effects of a model with varying intercepts and slopes. This post was very helpful in terms of the appropriate tidybayes syntax. However, adding varying slopes to the model results in unexpected behaviour from gather_draws; specifying r_id[id,Intercept] or r_id[id,time] makes no difference and the resulting plot does not correspond to the random effects for either the intercept or the slope.

I am using @JohsEnevoldsen’s reprex with the only modification the addition of random slopes for time:

library(tidyverse)
library(brms)
library(patchwork)

set.seed(1)
d <- 
  expand_grid(
    id = factor(1:10),
    time = 1:4
  ) %>% 
  mutate(
    z = rep(rnorm(10, sd = 0.5), each = 4),
    Y = 1.5 + 0.5 * time + z + rnorm(40, sd = 0.2),
  )

M1 <- brm(Y ~ 1 + time + (1 + time| id),
          data = d, family = gaussian(),
          iter = 1000,
          backend = "cmdstanr")


ranef(M1)
$id
, , Intercept

      Estimate Est.Error        Q2.5       Q97.5
1  -0.12182594 0.2141703 -0.54930118  0.29942240
2   0.06140011 0.1663396 -0.28008318  0.38041707
3  -0.33533598 0.1801811 -0.72191107 -0.01113552
4   0.44677237 0.1893477  0.07975836  0.86135812
5  -0.11140675 0.1807069 -0.47692082  0.23308177
6  -0.29537786 0.1831022 -0.67463705  0.04450496
7  -0.08908667 0.1802029 -0.46727428  0.23957137
8   0.37610969 0.1886941  0.02324225  0.77461247
9   0.25854348 0.1904829 -0.10047515  0.65007647
10 -0.17810860 0.1772599 -0.54263985  0.15150215

, , time

       Estimate  Est.Error        Q2.5       Q97.5
1  -0.127122340 0.07349263 -0.28767610 0.002235723
2   0.016880230 0.05182165 -0.08242850 0.124009700
3   0.002815841 0.05538314 -0.10253180 0.119227900
4   0.068828232 0.06146326 -0.04963366 0.194749350
5   0.043993295 0.05598179 -0.05237502 0.161750675
6  -0.042610859 0.05758365 -0.15697500 0.065933037
7   0.055915037 0.05640605 -0.03752997 0.176907225
8  -0.016174280 0.05552775 -0.13380868 0.085501485
9  -0.030336484 0.05775368 -0.16932777 0.071325792
10  0.016324616 0.05190020 -0.08350246 0.124140150

get_variables(M1)
plot_posterior_draws <- function(M) {
  
  rand.int <- gather_draws(M, r_id[id,Intercept]) %>% 
    mutate(id = factor(id)) %>% 
    ggplot(aes(x = .value, y=id)) +
    stat_halfeye() +
    labs(title="Random Intercepts")
  
  rand.slo <- gather_draws(M, r_id[id,time]) %>% 
    mutate(id = factor(id)) %>% 
    ggplot(aes(x = .value, y=id)) +
    stat_halfeye() +
    labs(title="Random Slopes")
  
  rand.int / rand.slo
}

plot_posterior_draws(M1)

rando
rando759×797 12.7 KB

@mjskay am I doing this wrong?

2

gather_draws will by default return both the random terms. For example,

r$> M1 |> gather_draws(r_id[id, term]) |> count(id, term)
# A tibble: 20 × 4
# Groups:   id, term, .variable [20]
      id term      .variable     n
   <int> <chr>     <chr>     <int>
 1     1 Intercept r_id       2000
 2     1 time      r_id       2000
 3     2 Intercept r_id       2000
 4     2 time      r_id       2000
 5     3 Intercept r_id       2000
 6     3 time      r_id       2000
 7     4 Intercept r_id       2000
 8     4 time      r_id       2000
 9     5 Intercept r_id       2000
10     5 time      r_id       2000
11     6 Intercept r_id       2000
12     6 time      r_id       2000
13     7 Intercept r_id       2000
14     7 time      r_id       2000
15     8 Intercept r_id       2000
16     8 time      r_id       2000
17     9 Intercept r_id       2000
18     9 time      r_id       2000
19    10 Intercept r_id       2000
20    10 time      r_id       2000

Your plotting function is then combining draws for both terms.
You could alter the function to something like

plot_posterior_draws <- function(M) {
  gather_draws(M, r_id[id, term]) %>%
    mutate(id = factor(id)) %>%
    ggplot(aes(x = .value, y = id)) +
    facet_wrap(~term) +
    stat_halfeye() +
    labs(title = "Random terms")
}

Or do explicit filtering of the gathered draws.

1 Like
3

I see! So in gather_draws(M, r_id[id, anything]), the anything bit does not already do the filtering, it can be literally anything (but not empty; that drops the column) and filtering has to be done separately downstream.

Thanks!