Computing a 95% prediction interval from brms model

Interfaces brms
1

Hi,
I am new to brms. I am doing a meta-analysis of proportions. I am able to get 95% intervals from the posterior distribution (see code below). But I also want a 95% prediction interval (example from bayesmeta package in attached image). I think it should work the same as I did below, but just with sampling from the predictive distribution. However I have trouble accessing the predictive distribution.
predict() or posterior_predict() do not do what I want or I am using them wrong.

Can anybody help me out here?
Thanks!!

data <- tribble(
  ~study,	~n,	~event,	
  "Buttermann 2004", 77, 2,	
  "Osterman 2006", 37, 0,	
  "Weinstein 2006, 2008", 253, 10,	
  "Peul 2007, 2008", 187, 2,	
  "McMorland 2010", 28,	0,	
  "Bailey 2020, 2021", 80, 2,	
  "Wilby 2021", 96,	0,	
  "Gadjradj 2022", 249,	8,	
  "Arts 2009", 166, 14,	
  "Teli 2009", 70, 2	
)

model <- brm(event | trials (n) ~ 1 + (1|study),
             data = data,
             family = binomial(link = "logit"))

as_draws_df(model, variable = c("b_Intercept")) %>% 
  rename(mu = b_Intercept) %>% 
  summarise(estimate = median(mu), 
            ci_lb = quantile(mu, 0.025), 
            ci_ub = quantile (mu, 0.975))

Screenshot 2023-11-16 at 21.48.38
Screenshot 2023-11-16 at 21.48.381378Ă—574 32.7 KB

2

What’s the issue with posterior_predict() precisely?

3

Hi Angelo,
thanks for your reply.

This is what I get for my specific use case without specifying anything.
Seem to predict the events for each study???

What I need however would be a distribution of the true proportion of events in a future study. I am imagining an output similar to what I get with: as_draws_df(model, variable = c(“b_Intercept”))
The mean/ median of the distribution would be the same, but it should be wider.

posterior_predict(model)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 2 1 3 3 0 3 0 9 4 0
[2,] 1 0 7 2 0 3 1 9 11 0
[3,] 0 0 21 0 0 2 0 7 20 4
[4,] 0 0 12 0 0 4 0 14 16 4
[5,] 0 0 12 2 0 0 3 11 15 0
[6,] 4 0 9 4 0 10 0 12 8 1
[7,] 4 0 5 3 0 1 0 10 17 2
[8,] 3 1 8 0 3 1 0 7 13 0
[9,] 2 0 15 4 0 2 0 10 10 3
[10,] 0 1 8 3 0 6 0 9 9 1
[11,] 0 0 6 3 0 2 2 6 14 4
[12,] 4 0 5 0 2 0 0 8 21 1
[13,] 0 2 12 2 1 8 0 11 17 3
[14,] 4 0 14 4 0 1 0 6 13 1
[15,] 1 0 10 2 0 2 2 6 10 1
[16,] 2 1 5 0 0 2 2 12 12 2
[17,] 2 1 7 4 1 1 3 10 14 4
[18,] 2 0 8 4 0 0 1 10 35 6
[19,] 1 0 7 2 0 5 1 5 9 0
[20,] 3 0 9 1 1 4 2 5 19 3
[21,] 2 0 7 1 1 1 0 8 12 0
[22,] 2 0 7 3 2 2 1 7 10 3
[23,] 1 1 15 2 0 1 1 7 13 0
[24,] 3 1 5 5 0 0 1 11 13 0
[25,] 5 0 4 3 0 3 4 7 7 1
[26,] 6 1 13 4 3 1 2 4 8 0
[27,] 3 3 9 5 0 2 2 15 17 10
[28,] 2 1 11 6 3 4 2 15 24 2
[29,] 1 1 8 2 0 5 0 7 10 0
[30,] 2 0 8 2 0 0 2 7 13 1
[31,] 0 0 2 1 2 0 1 7 12 2
[32,] 2 1 8 3 1 0 1 5 14 1
[33,] 3 1 10 6 4 3 1 3 10 2
[34,] 4 0 6 1 0 2 1 1 8 2
[35,] 0 0 6 4 0 2 1 11 7 0
[36,] 1 1 9 1 1 1 1 5 8 2
[37,] 1 1 4 7 0 1 2 8 7 0
[38,] 8 1 5 2 1 3 0 7 11 1
[39,] 0 1 15 3 0 1 1 5 15 8
[40,] 2 0 3 2 0 4 0 4 4 1
[41,] 3 1 7 6 2 2 0 3 11 3
[42,] 0 0 11 0 1 2 3 8 8 2
[43,] 2 0 7 4 1 2 0 4 4 0
[44,] 0 3 10 1 1 3 11 9 18 1
[45,] 2 2 12 4 0 0 0 3 19 6
[46,] 2 1 10 7 0 2 0 8 11 7
[47,] 0 0 10 7 1 1 5 7 11 0
[48,] 0 0 9 1 0 3 0 9 13 2
[49,] 1 1 9 5 0 5 3 8 11 1
[50,] 2 0 11 4 0 4 2 23 11 2
[51,] 1 3 15 1 1 2 2 5 18 1
[52,] 3 0 10 1 1 3 1 6 13 2
[53,] 1 1 5 1 0 1 2 7 17 3
[54,] 2 4 14 0 0 2 0 9 10 3
[55,] 3 0 15 13 1 1 1 5 6 6
[56,] 2 0 9 1 2 3 1 7 14 6
[57,] 5 1 8 3 0 7 2 7 13 3
[58,] 2 0 4 4 1 3 0 2 14 3
[59,] 1 2 9 5 2 3 0 6 15 1
[60,] 4 0 8 3 2 5 1 5 11 1
[61,] 2 0 5 2 1 1 0 8 10 0
[62,] 1 2 14 9 0 2 0 14 20 3
[63,] 2 0 7 4 3 4 2 11 16 2
[64,] 1 0 16 2 3 1 0 6 8 4
[65,] 3 1 8 4 0 4 0 9 9 3
[66,] 6 0 8 3 0 2 2 8 10 2
[67,] 1 0 8 2 0 1 0 8 13 5
[68,] 2 0 18 4 0 3 0 12 14 1
[69,] 3 0 9 1 0 1 0 11 9 7
[70,] 0 1 7 2 0 3 0 10 17 2
[71,] 1 1 11 4 2 4 3 10 6 0
[72,] 0 0 3 4 0 0 0 5 19 0
[73,] 5 0 9 1 1 2 3 11 15 2
[74,] 4 1 3 3 2 5 0 14 17 2
[75,] 3 0 13 4 1 6 1 17 13 2
[76,] 4 1 5 4 0 4 5 5 20 1
[77,] 2 2 10 5 2 5 1 7 2 4
[78,] 2 0 6 1 0 2 2 13 5 0
[79,] 2 1 10 5 3 2 5 8 2 5
[80,] 2 0 12 4 2 5 3 6 17 4
[81,] 3 5 8 2 1 3 0 8 3 1
[82,] 1 0 1 4 1 0 1 6 10 0
[83,] 2 0 9 6 0 6 1 6 6 4
[84,] 3 2 3 3 0 0 0 10 16 1
[85,] 1 0 20 3 0 1 0 7 7 0
[86,] 1 0 6 4 1 0 2 6 13 1
[87,] 2 0 8 7 1 2 4 10 14 2
[88,] 1 1 5 0 0 5 4 11 11 2
[89,] 2 0 7 6 0 1 0 6 7 0
[90,] 1 1 11 0 4 1 0 12 18 5
[91,] 3 1 6 0 1 0 0 5 16 5
[92,] 2 1 12 2 0 1 2 6 27 0
[93,] 2 0 14 8 2 3 0 17 16 2
[94,] 2 0 8 5 0 4 1 14 5 1
[95,] 1 1 8 2 0 3 1 10 11 0
[96,] 3 0 7 4 0 0 0 11 16 3
[97,] 5 0 9 4 1 3 0 6 7 0
[98,] 2 0 11 3 1 0 0 5 6 1
[99,] 2 0 8 3 0 5 1 14 19 1
[100,] 0 4 15 4 0 0 3 2 12 2
[ reached getOption(“max.print”) – omitted 3900 rows ]

bayesmeta
bayesmeta1150Ă—1002 21.8 KB

4

Hi Florian,

That’s exactly what posterior_predict() should do!
The model does model events based on an event rate which variates by study according to a binomial distribution:

E_s|N_s \sim \mathcal{Binomial}(p_s, N)\\ log(\frac{p_s}{1-p_s}) = \beta_0 + \mathcal{N}(0,\sigma^2_s)

with E_s and N_s being the number of events and the denominator of each study s, \beta_0 the grand mean rate, \sigma^2_s the variance of the random deviates of each study from the mean rate, and finally, p_s the expected study event rate. From your answer, I believe you want to get this last (p_s), not the predicted number of events E_s.
If that’s the case, you want to use posterior_epred() instead, which gives the expected value distribution by observation, not the number of events.

1 Like