Getting started with bicountbrms

library(brms)
library(bicountbrms)

This vignette demonstrates how to use one family, binegbin(), from simulation to prediction: counts are simulated from known parameters, fitted, and the credible intervals checked against the values that generated them. The other joint families (bipois(), and the partially observed binegbin_partialobs() and bipois_partialobs()) follow the same pattern; the README compares them.

Everything here is a brms custom family (Bürkner 2017), so the fitting, prediction and model-comparison interfaces are those of brms. The non-linear formula syntax used further down is documented in Bürkner (2018).

The generative model

Each observed count is the sum of a latent count that both sources recorded and one that only that source recorded. This construction is trivariate reduction: Holgate (1964) and Karlis and Ntzoufras (2003) use it with Poisson components, Kirkpatrick and Neale (2016) and Kirkpatrick (2022) with negative-binomial components. binegbin() takes all three latent counts to be independent and negative binomial:

$\begin{aligned} N_s &\sim \mathrm{NB2}(\mu,\ \phi_s) \\ N_1 &\sim \mathrm{NB2}(\lambda_1,\ \phi_{x1}) \qquad N_2 \sim \mathrm{NB2}(\lambda_2,\ \phi_{x2}) \\[4pt] y_1 &= N_s + N_1 \qquad y_2 = N_s + N_2 \end{aligned}$

$N_s$ is the latent count recorded by both sources; $N_1$ and $N_2$ are the two source-exclusive components, the counts recorded by one source only. $N_s$ induces the correlation between the pair. It is never observed, and is marginalised out analytically in the likelihood.

NB2(m, phi) is neg_binomial_2 in Stan and dnbinom(size = phi, mu = m) in R: mean m, variance m + m^2/phi. Larger phi means less overdispersion.

Six dpars, all with a log link:

dpar	role
`mu`	rate of the shared component, $\mu$
`lambdaone`	rate of the exclusive component of the first source, $\lambda_1$
`lambdatwo`	rate of the exclusive component of the second source, $\lambda_2$
`shapes`	dispersion $\phi_s$ of the shared component
`shapexone`	dispersion $\phi_{x1}$ of the exclusive component of the first source
`shapextwo`	dispersion $\phi_{x2}$ of the exclusive component of the second source

brms::custom_family() rejects any family whose dpars do not include one named literally mu. The requirement is on the name, not on the position: mu may appear anywhere in the vector. Here it is bound to the rate of the shared component, so it is not the mean of either response, nor the mean of their difference. E[y1] = mu + lambdaone.

The two rates are spelled lambdaone/lambdatwo rather than lambda1/lambda2 because custom_family() rejects dpar names ending in a digit. The maths below writes them $\lambda_1$ and $\lambda_2$ .

All six use a log link, the conventional log-linear rate parameterisation for this construction (Karlis and Ntzoufras 2003).

Simulate from known parameters

set.seed(20260724)

n <- 400
truth <- list(
  mu        = 8,
  lambdaone = 3,
  lambdatwo = 2,
  shapes    = 4,
  shapexone = 6,
  shapextwo = 2
)

n_shared <- rnbinom(n, size = truth$shapes, mu = truth$mu)
n1       <- rnbinom(n, size = truth$shapexone, mu = truth$lambdaone)
n2       <- rnbinom(n, size = truth$shapextwo, mu = truth$lambdatwo)

dat <- data.frame(
  y1 = n_shared + n1,
  y2 = n_shared + n2
)

str(dat)
#> 'data.frame':    400 obs. of  2 variables:
#>  $ y1: num  11 8 7 15 5 13 9 15 13 8 ...
#>  $ y2: num  11 7 5 16 6 11 9 16 9 7 ...
cor(dat$y1, dat$y2)
#> [1] 0.8518386

The correlation is positive and substantial because both counts include the same n_shared.

Fit

Two things are specific to this package.

The second count is supplied through vint(). brms::custom_family() declares a single response column, so only y1 can be the response. y2 is passed alongside as supplementary integer data, and the Stan signature of the family reads it from there.

stanvars is not optional. binegbin_stanvars() injects the binegbin_lpmf Stan function into the functions block of the model. Without it, the generated model will not compile.

fit <- brm(
  bf(
    y1 | vint(y2) ~ 1,
    lambdaone     ~ 1,
    lambdatwo     ~ 1,
    shapes        ~ 1,
    shapexone     ~ 1,
    shapextwo     ~ 1
  ),
  family   = binegbin(),
  stanvars = binegbin_stanvars(),
  data     = dat,
  prior    = c(
    prior(normal(2, 1), class = "Intercept"),
    prior(normal(2, 1), class = "Intercept", dpar = "lambdaone"),
    prior(normal(2, 1), class = "Intercept", dpar = "lambdatwo"),
    prior(normal(2, 1), class = "Intercept", dpar = "shapes"),
    prior(normal(2, 1), class = "Intercept", dpar = "shapexone"),
    prior(normal(2, 1), class = "Intercept", dpar = "shapextwo")
  ),
  chains  = 2,
  iter    = 1500,
  refresh = 0,
  backend = backend,
  seed    = 1
)

One weakly-informative prior serves all six, since every dpar here is a log-linked positive quantity of similar magnitude. On the natural scale normal(2, 1) is lognormal, with 95% of its mass between 1 and 52.

The six truths span 0.69 to 2.08 on the log scale, so no single prior sits away from all of them: the true shared rate falls near the prior mean, the true second exclusive rate 1.3 prior SDs below it. Recovery here is therefore not on its own evidence that the prior is uninfluential. The article Mapping native parameters to interpretable coordinates tests that directly, shifting a prior median sevenfold and moving the corresponding posterior median by 0.005.

To adapt these priors, shift the mean to the scale of the counts rather than increasing the SD. All six dpars are log-linked, so a normal prior is lognormal on the natural scale, and increasing its SD moves mass to implausible values rather than making the prior neutral (Smith et al. 2020, supplement 3, Fig. S3-1).

Recovery of the generating values

All six dpars are log-linked, so each posterior intercept exponentiates back onto the natural scale.

draws <- as_draws_df(fit)

pars <- c(
  mu        = "b_Intercept",
  lambdaone = "b_lambdaone_Intercept",
  lambdatwo = "b_lambdatwo_Intercept",
  shapes    = "b_shapes_Intercept",
  shapexone = "b_shapexone_Intercept",
  shapextwo = "b_shapextwo_Intercept"
)

recovery <- do.call(rbind, lapply(names(pars), function(p) {
  x <- exp(draws[[pars[[p]]]])
  data.frame(
    dpar  = p,
    truth = truth[[p]],
    est   = median(x),
    lower = unname(quantile(x, 0.025)),
    upper = unname(quantile(x, 0.975))
  )
}))
recovery$covered <- with(recovery, truth >= lower & truth <= upper)

knitr::kable(recovery, digits = 2)

dpar	truth	est	lower	upper	covered
mu	8	7.91	7.14	8.61	TRUE
lambdaone	3	3.33	2.81	3.97	TRUE
lambdatwo	2	2.29	1.77	2.94	TRUE
shapes	4	3.58	2.68	4.73	TRUE
shapexone	6	7.16	3.31	21.90	TRUE
shapextwo	2	3.13	1.44	10.05	TRUE

The three rates should land tightly on their true values. The three dispersions are estimated from an aggregate mean–variance mismatch rather than from any directly observed quantity, so their intervals are wider; the dispersions of the two exclusive components especially, since they are identified only through the part of the spread in the pair that the shared component cannot explain. Wide but covering is the expected result here, not a warning sign.

Predict and score

posterior_predict() draws new y1 conditional on the observed y2 in each row. It samples the discrete conditional distribution of $N_s$ given y2, then adds a fresh draw of the exclusive component — the exact conditional, not an approximation.

yrep <- posterior_predict(fit, ndraws = 200)
dim(yrep)
#> [1] 200 400

data.frame(
  quantity = c("mean", "sd", "max"),
  observed = c(mean(dat$y1), sd(dat$y1), max(dat$y1)),
  predicted = c(mean(yrep), mean(apply(yrep, 1, sd)), mean(apply(yrep, 1, max)))
) |>
  knitr::kable(digits = 2)

quantity	observed	predicted
mean	11.25	11.24
sd	5.59	5.43
max	37.00	33.23

log_lik() gives the pointwise log-likelihood of the joint pair, evaluated by an independent R implementation of the same marginalisation sum that the Stan function computes. It feeds loo() and waic() in the usual way (Vehtari et al. 2017). The Pareto $k$ diagnostic reported alongside the estimate flags observations whose importance weights are unreliable (Vehtari et al. 2024).

ll <- log_lik(fit)
dim(ll)
#> [1] 1500  400

loo(fit)
#> 
#> Computed from 1500 by 400 log-likelihood matrix.
#> 
#>          Estimate   SE
#> elpd_loo  -2202.8 22.6
#> p_loo         5.8  0.6
#> looic      4405.7 45.2
#> ------
#> MCSE of elpd_loo is 0.1.
#> MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.5]).
#> 
#> All Pareto k estimates are good (k < 0.69).
#> See help('pareto-k-diagnostic') for details.

Quantities returned by `posterior_epred()`

Both families return the same quantity:

$\mathrm{E}[y_1 \mid y_2] = \mathrm{E}[N_s \mid y_2] + \lambda_1,$

the expected first count conditional on the observed second count, not the marginal expectation $\mathrm{E}[y_1]$ .

For bipois(), $\mathrm{E}[N_s \mid y_2]$ is closed form, because conditioning a sum of independent Poissons on its total gives a Binomial split. For binegbin() there is no such shortcut, so the conditional is evaluated over its support $k = 0, \ldots, y_2$ and summed. Both are exact, and each agrees with the mean of its own posterior_predict() draws to within Monte Carlo error.

The conditional expectation is returned for every row, including — under binegbin_partialobs() and bipois_partialobs() — those where $y_1$ was never observed. Those two constructors exist so that the fitted model can impute $y_1$ on those rows.

posterior_epred(fit) needs no special handling: it dispatches to the family method, as do fitted() and conditional_effects().

bipois() — the non-overdispersed Poisson sibling. Same construction, three dpars instead of six. Use it when the margins are not overdispersed; compare the two with loo().
The symmetric model. The fit above let the dispersions of the two exclusive components differ. To impose $\phi_{x1} = \phi_{x2}$ instead, route both through one non-linear parameter — nlf(shapexone ~ shapexx), nlf(shapextwo ~ shapexx), shapexx ~ 1, with nl = TRUE — and set its prior with nlpar = "shapexx". That is the model this package fitted under a single shapex dpar before 0.10.0, term for term.
binegbin_partialobs() — for data where the first count is missing on some rows. It is the same family: same name, same six dpars, same likelihood and the same post-processing. What it adds is a second vint() integer, a 0/1 flag saying whether $y_1$ was recorded. Rows where it was not are scored by the integrated-out marginal of this same joint rather than dropped, so they still inform $\mu$ , $\phi_s$ , $\lambda_2$ , $\phi_{x2}$ and any group-level effects; afterwards the fit imputes the missing count conditional on the observed one. Note that $\lambda_1$ and $\phi_{x1}$ are then identified by the paired rows alone. Despite its former name, this is unrelated to the brms cens() addition term, which means a value known to lie in a set.
The $(M, f, \delta)$ reparameterisation — the midpoint $M$ of the two expected counts, congruence $f$ , and source bias $\delta$ , with the inverse dispersions converted to overdispersions $\kappa = 1/\sqrt{\phi}$ , where $\kappa = 0$ is the Poisson limit. It is fitted through nlf() rather than a separate family, and binegbin_mfd_to_dpars() converts between the two. Because $\delta = 0$ and $\kappa = 0$ are finite, interpretable nulls, these coordinates admit shrinkage priors towards the simpler model unless the data support otherwise (Simpson et al. 2017). The rate parameterisation has no such nulls. Worked end to end in the article Mapping native parameters to interpretable coordinates.
bipois_partialobs() — the same partial observation, equidispersed: the Poisson special case of binegbin_partialobs(). Its integrated-out marginal is closed form, $y_2 \sim \mathrm{Poisson}(\mu + \lambda_2)$ , since a sum of independent Poissons is Poisson. Use it when the margins are not overdispersed, rather than fitting binegbin_partialobs() with its dispersions pressed against the Poisson boundary.
Modelling the difference instead. Where only the disagreement $d = y_1 - y_2$ is of interest, or where truncation via resp_trunc() is required, the Skellam and discrete Laplace/normal families in the companion package skellambrms (Smith 2026) model $d$ directly. The difference of two independent Poisson counts is Skellam-distributed (Skellam 1946); for the Bayesian treatment of count differences generally, see Karlis and Ntzoufras (2006). Note what that reduction discards: the shared component cancels out of $d$ , so the midpoint and the congruence are no longer estimable.

References

Bürkner, Paul-Christian. 2017. “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.

Bürkner, Paul-Christian. 2018. “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal 10 (1): 395–411. https://doi.org/10.32614/RJ-2018-017.

Holgate, P. 1964. “Estimation for the Bivariate Poisson Distribution.” Biometrika 51 (1–2): 241–45. https://doi.org/10.2307/2334210.

Karlis, Dimitris, and Ioannis Ntzoufras. 2003. “Analysis of Sports Data by Using Bivariate Poisson Models.” Journal of the Royal Statistical Society: Series D (The Statistician) 52 (3): 381–93. https://doi.org/10.1111/1467-9884.00366.

Karlis, Dimitris, and Ioannis Ntzoufras. 2006. “Bayesian Analysis of the Differences of Count Data.” Statistics in Medicine 25 (11): 1885–905. https://doi.org/10.1002/sim.2382.

Kirkpatrick, Robert M. 2022. RMKdiscrete: Sundry Discrete Probability Distributions. https://CRAN.R-project.org/package=RMKdiscrete.

Kirkpatrick, Robert M., and Michael C. Neale. 2016. “Applying Multivariate Discrete Distributions to Genetically Informative Count Data.” Behavior Genetics 46 (2): 252–68. https://doi.org/10.1007/s10519-015-9757-z.

Simpson, Daniel, Håvard Rue, Andrea Riebler, Thiago G. Martins, and Sigrunn H. Sørbye. 2017. “Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors.” Statistical Science 32 (1): 1–28. https://doi.org/10.1214/16-STS576.

Skellam, J. G. 1946. “The Frequency Distribution of the Difference Between Two Poisson Variates Belonging to Different Populations.” Journal of the Royal Statistical Society 109 (3): 296. https://doi.org/10.2307/2981372.

Smith, A. N. H., D. Acuña-Marrero, P. Salinas-de-León, E. S. Harvey, M. D. M. Pawley, and M. J. Anderson. 2020. “Instantaneous Vs. Non-Instantaneous Diver-Operated Stereo-Video (DOV) Surveys of Highly Mobile Sharks in the Galápagos Marine Reserve.” Marine Ecology Progress Series 649: 111–23. https://doi.org/10.3354/meps13447.

Smith, Adam N. H. 2026. skellambrms: Skellam and Discrete-Difference brms Custom Families for Count Differences. https://doi.org/10.5281/zenodo.22231870.

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.

Vehtari, Aki, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry. 2024. “Pareto Smoothed Importance Sampling.” Journal of Machine Learning Research 25 (72): 1–58.