Returns a brms custom family for the joint distribution of a pair
of counts, (y1, y2), constructed via trivariate reduction: y1 = N_s + N1, y2 = N_s + N2, with N_s ~ Poisson(mu), N1 ~ Poisson(lambdaone), N2 ~ Poisson(lambdatwo)
mutually independent given their rates. All three rates are link = "log".
Modelling the pair jointly avoids regressing the difference on one of its
own components (a d = y1 - y2 ~ y2 design), which induces regression to
the mean.
Use this when both counts were recorded on every row. If the first count is
missing on some rows, use bipois_partialobs(), which is the same family
with an observation flag – same name, same three dpars, same likelihood,
same post-processing.
y1 is the response of the family; y2 is passed in as supplementary
integer data via the brms vint() addition term, since the brms
custom_family() machinery is built around a single declared response
column – see Details.
Use in a brm() call as: brm( bf(y1 | vint(y2) ~ ...), family = bipois(), stanvars = bipois_stanvars(), data = dat )
Usage
bipois()
bipois_stanvars()
log_lik_bipois(i, prep)
posterior_predict_bipois(i, prep, ...)
posterior_epred_bipois(prep)Value
bipois() returns a brms custom_family object. bipois_stanvars()
returns a stanvars object holding the Stan code for the corresponding
_lpmf. log_lik_bipois() returns a numeric vector of log-densities,
one per posterior draw, for observation i.
posterior_predict_bipois() returns an integer vector of simulated
y1 values, one per posterior draw, for observation i.
posterior_epred_bipois() returns a numeric matrix of E[y1 | y2], with
posterior draws in rows and observations in columns.
Details
Naming note. brms::custom_family() requires a dpar literally named
"mu" (stop2("All families must have a 'mu' parameter."),
unconditional), whatever the family actually calls that quantity. Here it is
bound to lambda_shared, the rate of the
component shared between y1 and y2 – not a mean of either
response individually. lambdaone (rate of the source-1-exclusive component)
and lambdatwo (rate of the source-2-exclusive component) are the other two
dpars, plainly named (no forced reinterpretation needed for those two).
The spelling lambdaone rather than lambda1. custom_family() rejects dpar names
ending in a digit (stop2("'dpars' should not end with a number.")), as
well as dots and underscores. The documentation therefore writes these
rates as \(\lambda_1\) and \(\lambda_2\) while the
code must spell them lambdaone/lambdatwo. See the notation table in
the package README.
y2 as supplementary data rather than a second response. The brms
custom_family() API supports exactly one declared response column
(Y) plus optional supplementary integer/real data (vint()/vreal()
addition terms) – the same mechanism used for, e.g., binomial trial
counts in the brms custom-families vignette. There is no
undeclared-response concept for a genuinely joint two-count
likelihood; vint(y2) is the correct fit for that gap, not a
workaround. This does mean y2 is not itself treated as
brms-modelled response data (no missing-value handling, no
resp_*() addition terms apply to it) – it is fixed, observed
per-row data, consistent with the fact that every row this constructor is
for comes from the paired (both-observed) subset.
Order of dpars matters for the generated Stan call. brms generates
target += bipois_lpmf(Y[n] | mu[n], lambdaone[n], lambdatwo[n], vint1[n], 1) – dpars in the order declared here, then the two vars
entries. bipois_stan_funs (stanfunctions via bipois_stanvars())
declares bipois_lpmf with exactly this argument order; changing the order
here without changing the Stan signature (or vice versa) silently swaps
which rate governs which count. The trailing 1 is the observation flag,
supplied as a literal because every row of a fully paired design has both
counts.
See also
bipois_partialobs() for partially observed pairs; binegbin()
for the overdispersed counterpart.