Predict risk given success/trials data


I’m fitting a logistic risk model on hospital data, given Region and Class predictors, with cases and denominators for each observation, with the aim of making simulations.

the model is:
brm(cases | trials(num) ~ (1 | Region / Class), family = binomial(), data = DF, chains = 8, cores = 8, iter = 8000)

I would like to predict the per hospital risk, but predict.brmsfit() produces only the predicted count of events, and does not take a type = 'response' argument as predict.glm() to predict risk. I could use posterior_linpred() and then invlogit() but it doesn’t include the residual error in the prediction, which I need.

What am I missing?

You want to use the fitted method.

Uhm, fitted() returns to me a non-integer (?) number of successes, not the probability of success:

fitted(model, nsamples = 1, summary = F) %>% as.vector() %>% summary()
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  0.4842   3.7819   9.2655  15.9890  17.0964 142.8172

Because you didnt set trials to 1. There are multiple discourse threads or github issues that explain how you do it.

I didn’t know about it. I try searching a lot on google before asking here but nothing came up.


Hmm strange. Essentially you do conditions = data.frame(ntrials = 1) if ntrials is the name of you trials variable in your data set.

The uncertainty will increase proportionally.

1 Like