I exported the parameters from a model fit with brms so that I can calculate posterior predictions from the model “by hand.” This is because I eventually want to generate predictions in a GUI app coded in a different language where I will be able to import the parameters but won’t have access to the `predict.brmsfit()`

method. However I believe that I am generating the predictions incorrectly. Can someone please take a look at my code and tell me what is incorrect?

My approach is to extract the posterior draws, multiply the beta coefficients times the x values to get posterior mu values, use the function `qnbinom()`

with those mu values and the posterior shape parameter values to get the appropriate negative binomial quantile, then take the median of the resulting vector to get the point estimate. I get very close to the correct answers returned by `predict.brmsfit()`

but there are a few discrepancies. Any help would be appreciated.

```
exdat <- structure(list(pace = c(18.43, 18.43, 18.43, 18.43, 18.43, 18.43,
41.64, 41.64, 41.64, 41.64, 41.64, 41.64, 0, 0, 0, 0, 0, 0, 8.09,
8.09, 8.09, 8.09, 8.09, 8.09, 21.42, 21.42, 21.42, 21.42, 21.42,
21.42, 44.72, 44.72, 44.72, 44.72, 44.72, 44.72, 0, 0, 0, 0,
0, 0, 9.18, 9.18, 9.18, 9.18, 9.18, 9.18, 23.19, 23.19, 23.19,
23.19, 23.19, 23.19, 46.82, 46.82, 46.82, 46.82, 46.82, 46.82,
0, 0, 0, 0, 0, 0, 9.57, 9.57, 9.57, 9.57, 9.57, 9.57, 23.97,
23.97, 23.97, 23.97, 23.97, 23.97, 48.08, 48.08, 48.08, 48.08,
48.08, 48.08),
plac = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.79,
1.79, 1.79, 1.79, 1.79, 1.79, 2.52, 2.52, 2.52, 2.52, 2.52, 2.52,
3.22, 3.22, 3.22, 3.22, 3.22, 3.22, 4.02, 4.02, 4.02, 4.02, 4.02,
4.02, 11.18, 11.18, 11.18, 11.18, 11.18, 11.18, 11.75, 11.75,
11.75, 11.75, 11.75, 11.75, 12.61, 12.61, 12.61, 12.61, 12.61,
12.61, 13.47, 13.47, 13.47, 13.47, 13.47, 13.47, 31.67, 31.67,
31.67, 31.67, 31.67, 31.67, 31.94, 31.94, 31.94, 31.94, 31.94,
31.94, 32.47, 32.47, 32.47, 32.47, 32.47, 32.47, 33.25, 33.25,
33.25, 33.25, 33.25, 33.25),
FLR = c(10L, 10L, 3L, 4L, 2L, 2L,
7L, 7L, 2L, 2L, 1L, 1L, 10L, 10L, 6L, 6L, 4L, 4L, 9L, 8L, 2L,
3L, 2L, 2L, 6L, 6L, 2L, 2L, 1L, 1L, 4L, 4L, 2L, 2L, 1L, 1L, 4L,
4L, 2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 1L, 1L, 1L,
1L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L, 1L, 1L, 2L, 2L, 1L,
1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L)),
row.names = c(NA, -84L), class = "data.frame")
options(mc.cores = 2, brms.backend = 'cmdstanr')
exmod <- brm(FLR ~ pace + plac,
data = dat, family = negbinomial(link = 'log', link_shape = 'log'),
chains = 2, iter = 2000, warmup = 1000, seed = 1110,
prior = prior(normal(0, 2), class = b))
percentiles <- c(0.5, 0.1, 0.9)
newdat <- expand.grid(pace = c(0, 5, 10, 20), plac = c(0, 5, 10, 20))
# brms predictions
brms_predictions <- predict(exmod, newdata = newdat, probs = percentiles, robust = TRUE)
# my interpretation of how to achieve the predictions
post <- as_draws_df(exmod)
my_pred_fn <- function(pace, plac) {
mus <- as.matrix(post[, 1:3]) %*% c(1, pace, plac)
sapply(percentiles, function(p) median(qnbinom(p, mu = exp(mus), size = post$shape)))
}
my_predictions <- t(mapply(my_pred_fn, newdat$pace, newdat$plac))
brms_predictions[, 3:5] - my_predictions # Close but not identical.
```