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.

nbbp takes two approaches to mitigating this problem. First, we use Stan’s built-in negative binomial log-density function neg_binomial_2_lpmf to compute the final size PMF. Let $\mathrm{NegativeBinomial}(x; R, k)$ be the negative binomial PMF evaluated at $x$ with mean $R$ and concentration $k$. The PMF of the final size distribution can be written

$$\mathrm{Pr}(c \mid R, k) = -\log(c) + \mathrm{NegativeBinomial}(c - 1; R c, k c)$$

Compared to alternative representations of the likelihood (for example, directly translating equation 9 of Blumberg & Lloyd-Smith (2013) to Stan code), using neg_binomial_2_lpmf appeared to lead to the lowest frequency of CDF values exceeding 1.

Second, when an invalid CDF value (i.e., greater than 1) is encountered, nbbp dynamically rescales the CDF in Stan. Let $c_\mathrm{max}$ be the maximum outbreak size for which we need to evaluate $\mathrm{Pr}(c \mid R, k)$. If $\sum_{x=1}^{c_\mathrm{max}} \mathrm{Pr}(x \mid R, k) > 1$, then we instead use:

$$\widetilde{\mathrm{Pr}}(c \mid R, k) = \frac{1 - \epsilon}{\sum_{x = 1}^{c_\mathrm{max}} \mathrm{Pr}(x \mid R, k)} \cdot \mathrm{Pr}(c \mid R, k)$$

As such, the rescaled CDF is safely less than one. In practice, we use Stan’s machine_precision() for $\epsilon$. Note that due to the implementation of the likelihood in Stan, $\widetilde{\mathrm{Pr}}$ is only used for likelihoods of (completely observed) censored and size-conditioned outbreaks, not partially-observed outbreaks.

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.