# Computing a 95% prediction interval from brms model

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


Whatâ€™s the issue with posterior_predict() precisely?

Hi Angelo,

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 ]

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