# Manually generate posterior predictions from negative binomial model fit with brms

**URL:** <https://discourse.mc-stan.org/t/manually-generate-posterior-predictions-from-negative-binomial-model-fit-with-brms/29563>\
**Category:** brms\
**Tags:** posterior-predictive, brms\
**Created:** [November 22, 2022, 6:07pm UTC](https://discourse.mc-stan.org/t/manually-generate-posterior-predictions-from-negative-binomial-model-fit-with-brms/29563 "2022-11-22T18:07:13Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![qdread](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/qdread/32/4665_2.png) [@qdread](https://discourse.mc-stan.org/u/qdread)\
**Post date:** [November 22, 2022, 6:07pm UTC](https://discourse.mc-stan.org/t/manually-generate-posterior-predictions-from-negative-binomial-model-fit-with-brms/29563/1 "2022-11-22T18:07:13Z")

</div>

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.

```nohighlight
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.

```

---

<div class="post-metadata">

**Author:** ![scholz](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/scholz/32/14298_2.png) [@scholz](https://discourse.mc-stan.org/u/scholz)\
**Post date:** [November 24, 2022, 9:54am UTC](https://discourse.mc-stan.org/t/manually-generate-posterior-predictions-from-negative-binomial-model-fit-with-brms/29563/2 "2022-11-24T09:54:53Z")

</div>

`predict.brmsfit()` wraps `posterior_predict`, which gives you draws from the posterior predictive distribution. It sounds like what you are doing is closer to `posterior_epred`, which gives you “expected values of the posterior predictive distribution" (see [Expected Value of Posterior vs. Posterior of Expected Value with epred - #3 by JimBob](https://discourse.mc-stan.org/t/expected-value-of-posterior-vs-posterior-of-expected-value-with-epred/28502/3) for the difference).  
In addition, remember that in bayesian modeling, we are usually modeling an expected value or median for each individual observation, so there is some aleatoric uncertainty in your predictions unless you just go with the posterior mean (which is again the predict vs epred thing).

---

<div class="post-metadata">

**Author:** ![qdread](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/qdread/32/4665_2.png) [@qdread](https://discourse.mc-stan.org/u/qdread)\
**Post date:** [December 14, 2022, 5:05pm UTC](https://discourse.mc-stan.org/t/manually-generate-posterior-predictions-from-negative-binomial-model-fit-with-brms/29563/3 "2022-12-14T17:05:47Z")

</div>

Thanks very much for this helpful clarification. Do you have any idea how I can reproduce the behavior of `posterior_predict()`?

---

<div class="post-metadata">

**Author:** ![scholz](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/scholz/32/14298_2.png) [@scholz](https://discourse.mc-stan.org/u/scholz)\
**Post date:** [December 14, 2022, 5:30pm UTC](https://discourse.mc-stan.org/t/manually-generate-posterior-predictions-from-negative-binomial-model-fit-with-brms/29563/4 "2022-12-14T17:30:10Z")

</div>

You can simply run each sample through the rng of your likelihood. Negbinom in your case. Just manually build you linear predictor and apply the inverse link to it before putting it into the rng. You might have to transform between the parameterization of the rng and the one brms uses as well.

---

<div class="post-metadata">

**Author:** ![qdread](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/qdread/32/4665_2.png) [@qdread](https://discourse.mc-stan.org/u/qdread)\
**Post date:** [December 14, 2022, 7:35pm UTC](https://discourse.mc-stan.org/t/manually-generate-posterior-predictions-from-negative-binomial-model-fit-with-brms/29563/5 "2022-12-14T19:35:51Z")

</div>

Thanks a lot. I believe your approach is the correct solution and I’ll mark it as such, but just to confirm, this is the new code I used to manually generate the draws from the posterior distributon and get the quantiles.

```plaintext
my_pred_fn <- function(pace, plac) {
  linpred <- as.matrix(post[, 1:3]) %*% c(1, pace, plac)
  draws <- rnbinom(n = nrow(post), mu = exp(linpred), size = post$shape)
  quantile(draws, probs = percentiles)
}

```

When I set the same random seed before `predict.brmsfit` and before my own version, they generate identical results:

```plaintext
set.seed(234)
brms_predictions <- predict(exmod, newdata = newdat, probs = percentiles, robust = TRUE)

set.seed(234)
my_predictions <- t(mapply(my_pred_fn, newdat$pace, newdat$plac))

all.equal(unname(brms_predictions[, 3:5]), unname(my_predictions)) #TRUE

```

---

<div class="post-metadata">

**Author:** ![scholz](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/scholz/32/14298_2.png) [@scholz](https://discourse.mc-stan.org/u/scholz)\
**Post date:** [December 14, 2022, 7:44pm UTC](https://discourse.mc-stan.org/t/manually-generate-posterior-predictions-from-negative-binomial-model-fit-with-brms/29563/6 "2022-12-14T19:44:27Z")

</div>

On mobile but looks good to me and if you get the same as predict that’s a good sign :)