Custom brms families for Bayesian modelling of pairs of counts that are positively correlated. Applications include counts recorded by two sources that measure the same quantity, such as two observers, two instruments, or two reporting channels, where the question is how much, and how systematically, the two sources disagree.
The pair bipois() or binegbin(), respectively.
Correlation between the two counts is induced via an unobserved shared latent component,
with an additional component each that only one source recorded.
The native parameters
— means and (for negative binomial) dispersions for each of the three
latent components — can be reparameterised into four interpretable
coordinates: the midpoint of the two expected counts, the congruence
between the two sources, the bias toward one of them, and (for negative
binomial) the overdispersion of each component.
The package also allows the model to be fit when one source is only
partially observed (the *_partialobs variants); posterior predictions
can then be made for the missing values of that source
based on the observed values of the other.
| both counts recorded | first count sometimes unrecorded | |
|---|---|---|
| equidispersed (Poisson latent counts) | bipois() |
bipois_partialobs() |
| overdispersed (negative-binomial latent counts) | binegbin() |
binegbin_partialobs() |
Each row of the table above is one family (Poisson or negative binomial).
The two constructors for each family return the same
likelihood and the same post-processing methods, and differ only in whether the
y1_obs) enters the likelihood as data or as a constant.
The indicator is not supplied in cases where both counts are fully observed.
For a row where both counts are observed, the likelihood term is the
full joint distribution of mu, shapes), the rate and
(for negative binomial) dispersion specific to each source (lambdaone,
shapexone, lambdatwo, shapextwo), and, through
lambdaone and lambdatwo, the bias between the two sources.
For a row where only mu,
shapes, lambdatwo, shapextwo) and not the parameters that pertain
only to lambdaone, shapexone).
# install.packages("pak")
pak::pak("anhsmith/bicountbrms")Stan and a C++ toolchain are required. On Windows, install
Rtools45.
Either rstan or cmdstanr works as the brms backend.
The second count is provided as supplementary integer data via the vint()
function, and the Stan function for the family is specified via the stanvars
argument. Both are required.
library(brms)
library(bicountbrms)
fit <- brm(
bf(y1 | vint(y2) ~ 1,
lambdaone ~ 1, lambdatwo ~ 1,
shapes ~ 1, shapexone ~ 1, shapextwo ~ 1),
family = binegbin(),
stanvars = binegbin_stanvars(),
data = dat
)Where the first count was not recorded on some rows, each row supplies a second
integer through vint(), in the order vint(y2, y1_obs), and the family and its
stanvars are the _partialobs constructors:
fit_po <- brm(
bf(y1 | vint(y2, y1_obs) ~ 1,
lambdaone ~ 1, lambdatwo ~ 1,
shapes ~ 1, shapexone ~ 1, shapextwo ~ 1),
family = binegbin_partialobs(),
stanvars = binegbin_partialobs_stanvars(),
data = dat
)y1_obs is 1 where the first count was recorded and 0 where it was not. On a
row with y1_obs == 0, y1 may hold any non-negative integer — 0 is the
conventional placeholder — because the likelihood term for that row does not
depend on it. Do not use NA: brms drops such rows before fitting, taking their
observed posterior_predict() and posterior_epred() then
return a draw and
brms leaves every distributional parameter of a custom family except mu
without a prior. Set them; see Choosing priors.
| Code | Math | Meaning |
|---|---|---|
y1 (the response) |
count from source 1 | |
y2 (via vint()) |
count from source 2 | |
| — | latent count recorded by both sources | |
| — |
|
latent counts only one source recorded |
mu |
rate of the shared component (not a response mean) | |
lambdaone, lambdatwo
|
|
rates of the two source-exclusive components |
shapes |
NB2 dispersion of the shared component | |
shapexone, shapextwo
|
|
NB2 dispersions of the two exclusive components |
y1_obs (via vint()) |
— | the _partialobs() indicator: 1 where 0 otherwise |
brms::custom_family() requires one distributional parameter to be named mu,
and rejects any name ending in a digit. Here mu is the rate of the shared
component rather than the mean of either count, and the two exclusive
rates are spelled out (one or two).
- Getting started — one fit, end to end: simulate, fit, check, predict.
- The families and their parameters — the construction, the moments, and a fitted demonstration that the dispersions of the two exclusive components can be told apart when they differ.
- Choosing priors — what brms leaves improper, and which of the
class,dparandnlparslots each prior belongs in. - A worked partially observed fit — imputing a count that was never recorded, scored against values withheld from the model.
- Mapping native parameters to interpretable coordinates — the midpoint of the two expected counts, congruence and source bias, with shrinkage priors.
- Migration and errata — what each release changed, and two corrected results.
Modelling the differences directly. Where only the disagreement
resp_trunc(), the companion package
skellambrms supplies Skellam,
discrete Laplace and discrete normal families that model pairedcountbrms; see
Migration and errata.
citation("bicountbrms")The package is archived on Zenodo under the concept DOI 10.5281/zenodo.22239120, which resolves to the current version. To cite the exact version an analysis used, take the version-specific DOI from the Zenodo record for that release.
Holgate (1964) states the trivariate-reduction representation and gives maximum-likelihood estimation for it; the regression form is due to Karlis and Ntzoufras (2003). Kirkpatrick and Neale (2016) and Kirkpatrick (2022) use the same construction with negative-binomial components.
This package was developed by Sea Through Science for a fisheries electronic-monitoring project led by Johanna Pierre of JPEC Consulting, with support from The Nature Conservancy, the Belize Blue Bonds for Ocean Conservation and the Belize Fund for a Sustainable Future.
Holgate, P. (1964). Estimation for the bivariate Poisson distribution. Biometrika 51(1–2), 241–245. https://doi.org/10.2307/2334210
Karlis, D. and Ntzoufras, I. (2003). Analysis of sports data by using bivariate Poisson models. Journal of the Royal Statistical Society: Series D (The Statistician) 52(3), 381–393. https://doi.org/10.1111/1467-9884.00366
Kirkpatrick, R. M. (2022). RMKdiscrete: Sundry discrete probability distributions. R package version 0.2. https://CRAN.R-project.org/package=RMKdiscrete
Kirkpatrick, R. M. and Neale, M. C. (2016). Applying multivariate discrete distributions to genetically informative count data. Behavior Genetics 46(2), 252–268. https://doi.org/10.1007/s10519-015-9757-z
Each article listed under Documentation ends with its own reference list.