Manually generate posterior predictions from negative binomial model fit with brms

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.

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 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).