Implementation details

Conditioning and multimodality of the likelihood surface

As can be seen in the vignette “nbbp,” the likelihood surface sometimes has a mode at 0. It can also be multimodal when conditioning is applied.

Conditioning on observation sizes

In particular, multimodality can be induced when conditioning on observing chains of at least some size. An example of this can be seen using the pneumonic plague data considered in “Advanced data” where the likelihood must be conditioned on observing at least two cases to match the data-collection process.

library(nbbp)
library(ggplot2)
data(pneumonic_plague)

pneumonic_plague[4] <- Inf

lnl_surface <- compute_likelihood_surface(
  pneumonic_plague,
  condition_geq = rep(2, 19),
  r_grid = seq(1.0, 1.35, length.out = 100),
  k_grid = exp(seq(log(0.01), log(9999), length.out = 100))
)

lnl_surface |>
  dplyr::mutate(
    log_dens = ifelse(log_dens >= max(log_dens) - 1, log_dens, NA)
  ) |>
  ggplot(ggplot2::aes(x = r, y = k, fill = log_dens)) +
  theme_minimal() +
  geom_tile() +
  scale_y_continuous(transform = "log") +
  scale_fill_viridis_c(option = "magma")

Inferring an unconstrained $R$

Some applications of the negative binomial branching process model have conditioned the likelihood on extinction to study the $R > 1$ case. In this way, it can be guaranteed that all chain sizes are finite, even when infinite chains are possible. However, Waxman and Nouvellet (2019) show that when conditioning there are two indistinguishable modes, one with $R > 1$ and one with $R < 1$. For this reason, we treat non-extinct chains separately in both stan code for inference and R code for simulation. In inference, a non-finite chain provides guarantees that $R > 1$, this is the necessary condition for a negative binomial branching process to never go extinct. Therefore, the presence or absence of such chains is itself informative (see the vignette “Advanced data” for an application). In simulation, this means that when $R > 1$, some samples may be Inf, indicating a non-extinct chain.

Mathematically, we partition the model outcomes into finite and infinite chains. That is, a chain either goes extinct (has finite size) or doesn’t, and this happens according to the extinction probability $\text{Pr}(\mathcal{E} \mid R, k)$. Formally, the possible chain sizes are $c \in \{0, 1, \dots\} \cup \infty$. We can write the likelihood for one chain size $c$ as $$\begin{equation} \text{Pr}(c \mid R, k) = \mathbb{I}(c \in \{0, 1, \dots\}) \left[ \text{Pr}(c \mid R, k, \mathcal{E}) \text{Pr}(\mathcal{E} \mid R, k) \right] + \mathbb{I}(c \in \{\infty\}) \left[ 1.0 - \text{Pr}(\mathcal{E} \mid R, k) \right] \end{equation}$$ where the chain size which can take finite values $c = 0, 1, \dots$ as well as the value “$c = \infty$”, $\mathbb{I}$ is the indicator function, and $\mathcal{E}$ ($\mathcal{E}^{\text{C}}$) is the event that the chain goes (does not go) extinct. For conditionally finite chains, we can write $$\begin{equation} \text{Pr}(c \mid R, k, \mathcal{E}) = \frac{\text{Pr}(c, \mathcal{E} \mid R, k)}{\text{Pr}(\mathcal{E} \mid R, k)} = \frac{\text{Pr}(c \mid R, k)}{\text{Pr}(\mathcal{E} \mid R, k)}. \end{equation}$$ Using this expression reduces the likelihood of chain size $c$ to $$\begin{equation} \text{Pr}(c \mid R, k) = \mathbb{I}(c \in \{0, 1, \dots\} ) \left[\text{Pr}(c \mid R, k) \right] + \mathbb{I}(c \in \{\infty\}) \left[ 1.0 - \text{Pr}(\mathcal{E} \mid R, k) \right]. \end{equation}$$

Numerical approximations

Probability of extinction

Except for very specific values of $k$, the extinction probability is not available in closed form. In stan, we use Newton’s method to obtain a numerical solution when it is needed.

CDF

While the PMF for the chain-size distribution is available in closed form, there is no closed form solution for the CDF.

Compromises for efficiency

For efficiency, when drawing random chain sizes the package uses base::sample() on a pre-computed, truncated, set of values. This enables rapidly sampling many realizations for a pair of $R, k$ values. Truncation is hard-coded at 1,000,000 which in practice appears to cover at least 99.999% of the probability mass. The most mass is unaccounted for when $R = 1$ and $k \to 0$.

Censoring

When brute-force computing the CDF (by summing the PMF), sometimes it evaluates to numbers slightly greater than 1. In practice, this appears to happen exclusively for $R < 1$ and the maximum seems to be less than $1 + 10^{-12}$. This is noticeable when censoring observations, where for chain size $c$ we need $\text{Pr}(C >= c) = 1 - \text{CDF}(c - 1)$. When the CDF exceeds one, this is negative and MCMC cannot appropriately explore the $R < 1$ region. To avoid these issues, the stan model code checks that the CDF does not exceed 1 over the range it needs to be evaluated. If it ever exceeds 1, we rescale so that the maximum CDF is 1.0 - machine_precision().

Incompletely observed chains

Incompletely-observed chains require either inferring the unobserved true size via data augmentation or marginalizing the likelihood over the true size. We adopt the latter approach, which is implemented for a binomial sampling model of observed chain sizes as discussed in the “Advanced data” vignette. We must therefore approximate the infinite sum over the unobserved true size $u$ given the observed size $c$ and a per-case observation probability $p$:

$$\mathrm{Pr}(c \mid R, k, p) = \sum_{u = c}^\infty \mathrm{Pr}(c \mid u, p) \mathrm{Pr}(u \mid R, k)$$ In practice, we must evalutate the finite sum up to some $u_\mathrm{max}$, leaving an error, $$\begin{aligned} \mathrm{Pr}(c \mid R, k, p) &= \sum_{u = c}^{u_\mathrm{max}} \mathrm{Pr}(c \mid u, p) \mathrm{Pr}(u \mid R, k) + \sum_{u = u_\mathrm{max} + 1}^{\infty} \mathrm{Pr}(c \mid u, p) \mathrm{Pr}(u \mid R, k)\\ &= \sum_{u = c}^{u_\mathrm{max}} \mathrm{Pr}(c \mid u, p) \mathrm{Pr}(u \mid R, k) + \epsilon \end{aligned}$$

Under the hood, nbbp uses a dynamic approximation to choose $u_\mathrm{max}$ such that the absolue error is no larger than a pre-specified $\epsilon$. As we show below, by choosing $u_\mathrm{max}$ in terms of the quantile of a negative binomial distribution, we can achieve this bound.

The bound can most readily be understood in two parts. First, and introducing a later-convenient factor of $p$, $$\begin{aligned} p \epsilon &= p \sum_{u = u_m + 1}^{\infty} \mathrm{Pr}(c \mid u, p) \mathrm{Pr}(u \mid R, k)\\ &\leq p\sum_{u = u_m + 1}^{\infty} \mathrm{Pr}(c \mid u, p) \end{aligned}$$

Rewriting the RHS, $$\begin{aligned} p \sum_{u = u_m + 1}^{\infty} \mathrm{Pr}(c \mid u, p) &= p \sum_{u = u_m + 1}^{\infty} {u \choose c} p^{c} (1 - p)^{u - c}\\ &= \sum_{u = u_m + 1}^{\infty} {u \choose c} p^{c + 1} (1 - p)^{u - c}\\ &= \sum_{x = x_m + 1}^{\infty} {x + c \choose c} p^{c + 1} (1 - p)^x\\ &= \sum_{x = x_m + 1}^{\infty} {x + n - 1 \choose c} p^{n} (1 - p)^x\\ &= \sum_{x = x_m + 1}^{\infty}\frac{(x + n - 1)!}{x! \, \Gamma(n)} p^{n} (1 - p)^x \end{aligned}$$

This can be recognized as the upper tail probability of a negative binomial distribution, in R’s parameterization, with prob parameter $p$ and size parameter $n = c + 1$, evaluated at $x_m + 1 = u_m - c + 1$. Thus, and using the more common lower tail, we can write $$p \epsilon \leq 1 - \mathrm{CDF}_\mathrm{NB}(u_m - c; c + 1, p)$$ Ergo, if we choose to set $x_m = u_m - c$ to the $(1 - p \epsilon)$ quantile of said negative binomial, we will keep the error no larger than $\epsilon$. That is, $$u_m - c := \mathrm{CDF}^{-1}_\mathrm{NB}(1 - p \epsilon; c + 1, p)$$

In practice, this bound appears quite lenient. Numerical experiments suggest that, at default package settings for selecting $\epsilon = 10^{-5}$, the error is perhaps roughly $9 \times 10^{-10}$ on average and typically (much) smaller.

Approaches for confidence intervals

While profiling of the likelihood is a common approach for generating confidence intervals for the negative binomial branching process model, we found that univariate profiling produced CIs with poor coverage properties for much of parameter space. Except near $R = 1$, the parametric bootstrap (repeatedly simulating data from the fit model and fitting to the simulated data) provides notably better coverage. Thus the default method for CI generation is a hybrid approach, using the parametric bootstrap unless the profile-based interval would be superior. This is implemented by checking whether the profiling confidence interval crosses 1, in which case it is used, or not, in which case parametric bootstrapping is used.