Skip to contents

Joint bivariate custom family with negative-binomial components for brms

Source: R/binegbin.R
binegbin.Rd

Overdispersed sibling of bipois(). Returns a brms custom family for the joint distribution of a pair of counts (y1, y2) via trivariate reduction with Negative-Binomial (rather than Poisson) latent components: y1 = N_s + N1, y2 = N_s + N2, with N_s ~ NB2(mu, shapes), N1 ~ NB2(lambdaone, shapexone), N2 ~ NB2(lambdatwo, shapextwo) mutually independent given their rates. NB2(m, phi) has mean m and variance m + m^2/phi (Stan neg_binomial_2; R dnbinom(size = phi, mu = m)).

Use this when both counts were recorded on every row. If the first count is missing on some rows, use binegbin_partialobs(), which is the same family with an observation flag – same name, same six dpars, same likelihood, same post-processing.

Six dpars: the three rates (mu = shared rate, lambdaone/lambdatwo = the rates of the two source-exclusive components) plus three dispersions – shapes for the shared component and shapexone/shapextwo for the two exclusive components. All six use link = "log". Supply the exclusive rates through a non-linear formula without an explicit exp() (the log link applies it): nlf(lambdaone ~ lamx) gives lambdaone = exp(lamx).

See the binegbin.R file header for why Negative-Binomial components are used instead of an observation-level random effect on bipois() – briefly, the random-effect version fails synthetic recovery, because with one observed pair but three latent deviates per unit the exclusive-component deviates act as residual absorbers.

Use in a brm() call as: brm( bf(y1 | vint(y2) ~ 1, mu ~ 1 + (1 | vessel), nlf(lambdaone ~ lamx), nlf(lambdatwo ~ lamx), lamx ~ 1, shapes ~ 1, shapexone ~ 1, shapextwo ~ 1, nl = TRUE), family = binegbin(), stanvars = binegbin_stanvars(), data = dat )

Usage

binegbin()

binegbin_stanvars()

log_lik_binegbin(i, prep)

posterior_predict_binegbin(i, prep, ...)

posterior_epred_binegbin(prep)

Value

binegbin() returns a brms custom_family object. binegbin_stanvars() returns a stanvars object holding the Stan code for the corresponding _lpmf. log_lik_binegbin() returns a numeric vector of log-densities, one per posterior draw, for observation i. posterior_predict_binegbin() returns an integer vector of simulated y1 values, one per posterior draw, for observation i. posterior_epred_binegbin() returns a numeric matrix of E[y1 | y2], with posterior draws in rows and observations in columns.

Details

Tying the dispersions of the two exclusive components. Up to 0.9.1, binegbin() had a single shapex, shared by both exclusive components, which imposed shapexone == shapextwo. Equal dispersions are a modelling choice rather than a property of the construction. 0.10.0 removes the constraint. To recover the constraint, route shapexone and shapextwo through one non-linear parameter:

bf(y1 | vint(y2) ~ 1,
   mu ~ 1 + (1 | vessel),
   nlf(shapexone ~ shapexx),
   nlf(shapextwo ~ shapexx),
   shapexx ~ 1, ..., nl = TRUE)

The resulting likelihood is term-for-term the pre-0.10.0 five-dpar one, and a package test pins that. Stored five-dpar fits keep working unchanged: their single shapex resolves to both shapexone and shapextwo when post-processed. Set a prior with nlpar = "shapexx" rather than dpar = "shapex".

Names forced by custom_family(). Identical conventions to bipois() – mu is the mandatory brms dpar name, here bound to the rate of the shared component (lambda_shared), not a mean of either response; y2 is supplied as supplementary integer data through vint() because custom_family() declares a single response column. See bipois() for the full explanation, including why the rates are spelled lambdaone/ lambdatwo in code but written \(\lambda_1\)/ \(\lambda_2\) in the documentation.

Order of dpars matters for the generated Stan call. brms generates target += binegbin_lpmf(Y[n] | mu[n], lambdaone[n], lambdatwo[n], shapes[n], shapexone[n], shapextwo[n], vint1[n], 1) – dpars in the order declared here, then the two vars entries. binegbin_stan_funs declares binegbin_lpmf with exactly this signature; reordering one without the other silently swaps which rate or dispersion governs which component. The trailing 1 is the observation flag, supplied as a literal because every row of a fully paired design has both counts.

See also

binegbin_partialobs() for partially observed pairs; bipois() for the equidispersed Poisson counterpart.