Can I return credible intervals from `predict`?

Interfaces brms
1

Please also provide the following information in addition to your question:

  • Operating System: Ubuntu 19.04
  • brms Version: 2.9.0 (Github)

I’m working with a categorical model of the form:
My model is:
mod <- brm(Reason ~ V1 + V2 +V3 + (1|ID),
data = Data ,
family = categorical(link = “logit”),
prior = priors)

I am using predict.brmsfit to return predicted probabilities of each ‘Reason’ with:
predict(mod, newdata = NULL, re_formula = ~ (1|ID),
summary = T)

and using cbind to attach the returned predicted probabilities to the original data frame for visualization.

Is it possible to also return 95% credible intervals surrounding the predicted probabilities?

2

The quantiles in the output are already boundaries of a credible interval by definition.

3

How do I return the quantiles? I’m only getting a single value for each response category (row)

4

ah I see. for categorical models use fitted().

5

Thanks!

fitted(mod, newdata = NULL, re_formula = ~ (1|ID), summary = T)

returned the complete matrix I was looking for

1 Like
6

Hello!

I’m doing something similar and found this answer. I work on a multinomial multilevel model of electoral preference with brm (I use the MRP method). With the model I generate a probability prediction on a database with 480 rows (let’s call it data_est).

In a first attempt, I made the prediction with the function

“predict(model, ndraws=2500,newdata=data_est,summary=T)”

and I obtain a table with the probabilities for each row of each category (P(Y=partyA) , P(Y=partyB),…). With them, I took the mean as a point estimator and the quantiles (.05 and .975) as the “credibility interval”, but they are very wide.

Seeing this post I made a prediction now with the “fitted” function in the “data_est” database and the summary gives me the estimate (“Estimate.P(Y=partyA)”, “Estimate.P(Y=partyB)”, . …), the error (“Est.Error.P(Y=partyA)”, “Est.Error.P(Y=partyB)”, …) and the quantiles corresponding to each category (“Q2.5. P(Y=partyA)” and “Q97.5.P(Y=partyA)”, …) per row. I have a question, to obtain my “general” estimator, is it correct to only obtain the (weighted) average? Something like

partyA=mean(“Estimate.P(Y=partyA)”)

or is there something else to do? And for the quantiles? (low_partyA=mean(“Q2.5.P(Y=partyA)”))?
or am i doing something wrong?

Thank you so much.

7

It sounds like you might want to use (…summary = F) and then do post-processing yourself (recommend tidybayes) to obtain quantities of interest.

8

hey! thanks for answering!

with “summary=T” , i get probabilities for every category (column) for each “individual” (row), estimates and quantiles, something like…

          Estimate.P(Y=pA)     Est.Error.P(Y=pA)      Q2.5.P(Y=pA)    Q97.5.P(Y=pA) ...
                   <dbl>                  <dbl>             <dbl>       <dbl>
 1                0.0742                 0.0337            0.0177       0.153
 2                0.0712                 0.0320            0.0161       0.149
 3                0.0632                 0.0284            0.0133       0.128
 4                0.120                  0.0504            0.0288       0.231
 5                0.0851                 0.0367            0.0196       0.170

…

with “summary=F” i obtain all ndraws with the posible result (1, 2, etc corresponding to 1=partyA, 2=partyB…)

      draw1         draw2        draw3        draw4...
1    7               9             2              7
2    1               7             9              9   
3    7               9             9              7
4    2               1             7              1
5    7               4             9              2

…

i think it’s ok to take the first one because i get the summary data and the estimated prediction to every party, but my doubt is about the quantiles, with “predict” function only obtain the estimated by row (not quantiles or “error”), and with this column i get mean and quantiles (in my first attempt like paul.buerkner said answers above), i get someting like this:

party    low (.025)           mean         up(.975)   
A           .04                .08            .16
B            .26               .40            .58

…
but this values are to wide.

For “fitted” function, if i summarise Estimates.P(Y=pA) obtain a similar value like my mean above, but my question is if it´s ok make the same (mean value) for Q2.5.P(Y=pA) and Q97.5.P(Y=pA)? because if i do it i get this:

party    low (.025)          mean         up(.975)
A            .06               .08            .11
B            .34               .40            .46

…
and this is more closed.

Tks!