Fitting a Negative Binomial distribution as a Gamma Mixture of Poisson distributions with brms

I have some data which I know should be approximately negative binomial distributed. These are counts, with corresponding ids. At the moment, I have a single observation for each id and no additional covariates but I am expecting to ultimately have repeated observations and want to account for covariates which differ by id.

For my current, simple version of the problem I can run a negative binomial regression and get a good fit to the overall data e.g.


dat <- rnbinom(100,size=1.66,prob=0.25)

dat %>%
  as_tibble() %>%
  mutate(id=row_number()) %>%
  rename(count=value) -> dat

#fit nbd
brms_nbd_model <- brm(count ~ 1,

This works fine with both my real, and simulated data.

The NBD distribution can be thought of as mixture of poisson distributions, where the poisson rate parameters are gamma distributed.

For my use case, I’m interested in how the rates of each id vary. So I am looking to fit the equivalent to the negative binomial model as mixture using brms.

I am struggling to do this, I can’t find an equivalent example of how to mix different types of distributions. Is this possible using the brms interface?

I have tried the following apporaches which don’t quite work:

  1. I can fit a poisson likelihood to my counts, and allow the rates to vary by id:
brms_mixed_model <- brm(count ~ (1|id),

This runs, but isn’t quite what I want as the rates by id are normally distributed about the average so it’s not equivalent to an NBD. Visually the fit is worse than the topline NBD based on pp_check.

  1. I have tried extending this approach to make use of the non-linear model syntax:
brms_mixed_model2 <- brm(bf(count ~ lambda,
                           lambda ~ (1|id),

I’m not 100% clear on how everything is interpreted here (does the family part mean lambda is treated as poisson?), but the results look equivalent to the previous attempt.

  1. I think what I need to do is to be able to specify that count is poisson, and that lambda is gamma distributed but it doesn’t look like I can do this using “mixture()”:
brms_mixed_model3 <- brm(bf(count ~ lambda,
                            lambda ~ (1|id),
#Error: Cannot mix families with real and integer support.
  1. I’ve also tried to do the same thing by using brmsformula without success:
complex_form <- bf(count ~ lambda,nl=TRUE) + poisson()+ lf(lambda ~ (1|id))+Gamma(link="log")

brms_mixed_model4 <- brm(complex_form,
#Error: Family 'gamma' requires response greater than 0.

In summary - is it possible to do what I’m trying to do directly in brms? If so, would anyone be able to point me to an example or modify my syntax above to move in the right direction?

It looks like there’s an example of how to do this directly in stan here but I would prefer to use brms if it is possible.

Thanks in advance for any help.

You may want to set a Gamma prior for lambda (=rate) and fit a Poisson. That would give you the negbin.

Maybe I am overlooking something, but why don’t you just use the negative binomial distribution as provides by brms?

Search for negbinomal here