We present the modeling framework implemented by the package. Section 1.1 outlines the bare-bones version of the model, which is elaborated on in Sections 1.2, 1.3 and 1.4. Section 1.5 extends the model and introduces multilevel modeling, treating infections as parameters, and accounting for population effects.
We now formulate the basic version of the model for one homogeneous population. The same model can be used for multiple regions or groups jointly. Suppose we observe a non-negative time series of count data \(Y = (Y_1, \ldots Y_n)\) for a single population. This could for example be daily death or case incidence. \(Y_t\) is modeled as deriving from past new infections \(i_s\), \(s < t\), and some parameter \(\alpha_t > 0\), a multiplier, which in most contexts represents an instantaneous ascertainment rate. The general model can be expressed as \[\begin{align} Y_t & \sim p(y_t , \phi), \tag{1.1} \\\\ y_t & = \alpha_t \sum_{s < t} i_s \pi_{t-s} \tag{1.2}, \end{align}\] where \(y_t\) is the expected value of the data distribution and \(\phi\) is an auxiliary parameter. \(\pi_{k}\) is typically the time distribution from an infection to an observation, which we refer to as the infection to observation distribution. More generally, however, \(\pi_k\) can be used to obtain any linear combination of past infections. New infections \(i_t\) at times \(t>0\) are modeled through a renewal equation, and are tempered by a non-negative parameter \(R_t\) which represents the reproduction number at time \(t\). Formally \[\begin{align} i_t &= R_t \sum_{s < t} i_s g_{t-s}, \tag{1.3} \end{align}\] where \(g_k\) is a probability mass function for the time between infections. The recursion is initialized with seeded infections \(i_{v:0}\), \(v < 0\), which are treated as unknown parameters. All parameters are assigned priors, i.e. \[\begin{equation} i_{v:0}, R, \phi, \alpha \sim p(\cdot), \end{equation}\] where \(R = (R_1, \ldots, R_n)\) and \(\alpha = (\alpha_1, \ldots, \alpha_n)\). The posterior distribution is then proportional to prior and likelihood, i.e. \[\begin{equation} p(i_{v:0}, R, \phi, \alpha \mid Y) \propto p(i_{v:0})p(R)p(\phi)p(\alpha) \prod_{t} p(Y_t \mid y_t, \phi). \end{equation}\] This posterior distribution is represented in a Stan program, and an adaptive Hamiltonian Monte Carlo sampler (Hoffman and Gelman 2014) is used to approximately draw samples from it. These samples allow for inference on the parameters, in addition to simulating data from the posterior predictive distribution.
Reproduction numbers \(R\) and multipliers \(\alpha\) can be modeled flexibly with Bayesian regression models, and by sharing parameters, are the means by which multiple regions or groups are tied together through multilevel models. One can, for example, model \(R\) as depending on a binary covariate for a control measure, say full lockdown. The coefficient for this can be partially pooled between multiple populations. The effect is to share information between groups, while still permitting between group variation.
As mentioned, \(Y_t\) is usually a count of some event type occurring at time \(t\). These events are precipitated by past infections. Prototypical examples include daily cases or deaths. \(\alpha_t\) is a multiplier, and when modeling count data, it typically is interpreted as an ascertainment rate, i.e. the proportion of events at time \(t\) that are recorded in the data. For case or death data this would be the infection ascertainment rate (IAR) or the infection fatality rate (IFR) respectively.
The multiplier \(\alpha\) plays a similar role for observations as \(R\) does for infections; tempering expected observations for time-specific considerations. As such, epidemia treats \(\alpha\) in a similar manner to reproductions number, and allows the user to specify a regression model for it. Section 1.4 discusses this in detail in the context of reproduction numbers, and this discussion is not repeated here. The model schematic](model-schematic.html) details the model for \(\alpha\), as well as for observational models in general.
The sampling distribution \(p(y_t, \phi)\) (Equation (1.1)) should generally be informed by parts of the data generating mechanism not captured by the mean \(y_t\): i.e. any mechanisms which may induce additional variation around \(y_t\). Options for \(p(y_t, \phi)\) include the Poisson, quasi-Poisson and negative-binomial families. The Poisson family has no auxiliary parameter \(\phi\), while for the latter two families this represents a non-negative dispersion parameter which is assigned a prior.
epidemia allows simultaneous modeling of multiple observation vectors. In this case, we simply superscript \(Y_t^{(l)}\), \(\alpha_{t}^{(l)}\) and \(\pi^{(l)}\), and assign independent sampling distributions for each type. Separate regression models are then specified for each multiplier \(\alpha_{t}^{(l)}\). Leveraging multiple observation types can often enhance a model. For example, high quality death data existed during the first wave of the Covid-19 pandemic in Europe. Case data gradually increased in reliability over time, and has the advantage of picking up changes in transmission dynamics much quicker than death data.
Infections \(i_t\) propagate over time through the discrete renewal equation (1.3). This is self-exciting: past infections give rise to new infections. The theoretical motivation for this lies in counting processes and is explained in more detail in Bhatt et al. (2020). The equation is connected to Hawkes processes and the Bellman Harris branching process (Bellman and Harris 1948, 1952; Mishra et al. 2020). Such processes have been used in numerous previous studies (Fraser 2007; Cori et al. 2013; Nouvellet et al. 2018; Cauchemez et al. 2008), and are also connected to compartmental models such as the SEIR model (Champredon, Dushoff, and Earn 2018).
Equation (1.3) implies that infections \(i_t\), \(t > 0\) are deterministic given \(R\) and seeded infections \(i_{v:0}\). epidemia sets a prior on \(i_{v:0}\) by first assuming that daily seeds are constant over the seeding period. Formally, \(i_{k} = i\) for each \(k \in \{v,\ldots 0\}\). The parameter \(i\) can be assigned a range of prior distributions. One option is to model it hierarchically; for example as \[\begin{align} i & \sim \text{Exp}(\tau^{-1}), \tag{1.4}\\\\ \tau & \sim \text{Exp}(\lambda_0) \tag{1.5}, \end{align}\] where \(\lambda_0 > 0\) is a rate hyperparameter. This prior is uninformative, and allows seeds to be largely determined by initial transmission rates and the chosen start date of the epidemic.
Several extensions to the infection model are possible in epidemia, including extending (1.3) to better capture dynamics such as super-spreading events, and also adjusting the process for the size of the remaining susceptible population. These extensions are discussed in Section 1.5.2 and 1.5.3 respectively. The basic infection model is shown in the model schematic.
Reproduction numbers are modeled flexibly. One can form a linear predictor consisting of fixed effects, random effects and autocorrelation terms, which is then transformed via a suitable link function. Formally \[\begin{equation} R = g^{-1}(\eta), \end{equation}\] where \(g\) is a link function and \(\eta\) is a linear predictor. In full generality, \(\eta\) can be expressed as \[\begin{equation} \eta = \beta_0 + X \beta + Z b + Q \gamma, \tag{1.6} \end{equation}\] where \(X\) is an \(n \times p\) model matrix, \(Z\) is an \(n \times q\) model matrix for the \(q\)-vector of group-specific parameters \(b\). \(Q\) is an \(n \times r\) model matrix for the \(r\)-vector of autocorrelation terms. The columns of \(X\) are predictors explaining changes in transmission. These could, for example, be binary vectors encoding non-pharmaceutical interventions, as in Flaxman et al. (2020). A number of families can be used for the prior on \(\beta\), including normal, cauchy, and hierarchical shrinkage families. The parameters \(b\) are modeled hierarchically as \[\begin{equation} b \sim N(0, \Sigma), \end{equation}\] where \(\Sigma\) is a covariance matrix that is itself assigned a prior. The particular form for \(\Sigma\), as well as its prior is discussed in more detail in this article. These partially-pooled parameters are particularly useful when multiple regions are being modeled simultaneously. In this case, they allow information on transmission rates to be shared between groups.
\(Q\) is a binary matrix specifying which of the autocorrelation terms in \(\gamma\) to include for each period \(t\). Currently, epidemia supports only random walk processes. However multiple such processes can be included, and can have increments that occur at a different time scale to \(R\); for example weekly increments can be used.
Choosing an appropriate link function \(g\) is difficult. \(R_t\) is non-negative, but is clearly not able to grow exponentially: regardless of the value of the linear predictor \(\eta_t\), one expects \(R_t\) to be bounded by some maximum value \(K\). In other words, \(R_t\) has some carrying capacity. One of the simplest options for \(g\) is the log-link. This satisfies non-negativity, and also allows for easily interpretable effect sizes; a one unit change in a predictor scales \(R_t\) by a constant factor. Nonetheless, it does not respect the carry capacity \(K\), often placing too much prior mass on large values of \(R_t\). With this in mind, epidemia offers an alternative link function satisfying \[\begin{equation} g^{-1}(x) = \frac{K}{1 + e^{-x}} \tag{1.7}. \end{equation}\] This is a generalization of the logit-link, and we refer to it as the scaled-logit.
Various extensions to the basic model just presented are possible, including multilevel modeling, adding variation to the infection process, and explicitly accounting for population effects. These are discussed in turn.
Consider modeling the evolution of an epidemic across multiple regions or populations. Of course, separate models can be specified for each group. This approach is fast as each model can be fit in parallel. Nonetheless, often there is little high quality data for some groups, particularly in the early stages of an epidemic. A joint model can benefit from improved parameter estimation by sharing signal across groups. This can be done by partially or fully pooling effects underlying reproduction numbers \(R\).
We give an example for concreteness. Suppose the task is to infer the effect of a series of \(p\) control measures on transmission rates. Letting \(R^{(m)}\) be the vector of reproduction numbers for the \(m\) group, one could write \[\begin{equation} R^{(m)} = g^{-1}\left( \beta_0 + b_0^{(m)} + X^{(m)} (\beta + b^{(m)}) \right), \end{equation}\] where \(X^{(m)}\) is an \(n \times p\) matrix whose rows are binary vectors indicating which of the \(p\) measures have been implemented in the \(m\) group at that point in time. The parameters \(b_{0}^{(m)}\) allow each region to have its own initial reproduction number \(R_0\), while \(b^{(m)}\) allow for region-specific policy effects. These parameters can be partially pooled by letting \[\begin{equation} (b_0^{(m)}, b^{(m)}) \sim N(0,\tilde{\Sigma}), \end{equation}\] for each \(m\), and assigning a hyperprior to the covariance matrix \(\tilde{\Sigma}\).
In addition to hierarchical modeling of parameters making up \(R\), seeded infections can also be modeled hierarchically. Equations (1.4) and (1.5) are replaced with \[\begin{align} i^{(m)} &\sim \text{Exp}(\tau^{-1}), \\\\ \tau & \sim \text{Exp}(\lambda_0), \end{align}\] where \(i^{(m)}\) is the daily seeded infections for the \(m\) group.
Recall the renewal equation (Equation (1.3)) which describes how infections propagate in the basic model. Infections \(i_t\) for \(t > 0\) are a deterministic function of seeds \(i_{v:0}\) and reproduction numbers \(R\). If infections counts are large, then this process may be realistic enough. However, when infection counts are low, there could variation in day-to-day infections caused by a heavy tailed offspring distribution and super-spreader events. This may cause actual infections to deviate from those implied by the renewal equation. Although the expected number of offspring of any given infection is driven by \(R\), in practice the actual number of offspring can exhibit considerable variation around this. To capture this randomness, replace Equation (1.3) with \[\begin{align} i_t &\sim p(i_t', d), \tag{1.8} \\\\ i_{t}' &= R_t \sum_{s < t} i_s g_{t-s}. \end{align}\] This treats \(i_t\) as latent parameters which must be sampled. Instead, the is described by the renewal equation. \(p(i_t', d)\) is parameterised by the mean and the coefficient of dispersion \(d\), which is assigned a prior. This extension can be motivated formally through counting processes. Please see Bhatt et al. (2020) for more details.
Nothing in Equation (1.3) prevents cumulative infections from exceeding the total population size \(P\). In particular if \(R_t > 1\) then infections can grow exponentially over time. This does not always present a problem for modeling. Indeed the posterior distribution usually constrains past infections to reasonable values. Nonetheless, forecasting in the basic model will be unrealistic if projected infections grow too large. As the susceptible population diminishes, the transmission rate is expected to fall.
epidemia can apply a simple transformation to ensure that cumulative infections remain bounded by \(P\), and that transmission rates are adjusted for changes in the susceptible population. Let \(S_t \in [0,P]\) be the number of susceptible individuals in the population at time \(t\). Just like infections, this is treated as a continuous quantity. \(S_t\) consists of those who have not been infected by time \(t\), and have not been removed from the susceptible class by other means; i.e. vaccination.
Let \(i'_t\) denote from the model. This is given by (1.3) in the basic model or by (1.8) if the extension of Section 1.5.2 is applied. These are interpreted as the number of infections if the entire population were susceptible. These are adjusted with
\[\begin{equation} i_t = S_{t-1}\left( 1 - \exp\left(-\frac{i'_t}{P}\right)\right). \tag{1.9} \end{equation}\]
The motivation for this is provided in Bhatt et al. (2020). Equation (1.9) satisfies intuitive properties: if \(i'_t = 0\) then \(i_t = 0\), and as \(i'_t \to \infty\) we have that \(i_t \to S_{t-1}\). All infections at time \(t\) are then removed from the susceptible population, so that \[\begin{equation} S_t = S_{t-1} - i_t \tag{1.10} \end{equation}\] We are left to define \(S_{v-1}\), the susceptible population the day before modeling begins. If this is the start of an epidemic, it is natural to take \(S_{v-1} = P\). Nonetheless, it is often of interest to begin modeling later, when a degree of immunity already exists exists within the population. In this case, epidemia allows the user to assign a prior distribution to \(S_{v-1} / P\). This must lie between \(0\) and \(1\).
Previous infection is one avenue through which individuals are removed from the susceptible population. Immunity can also be incurred through vaccination. epidemia provides a basic way to incorporate such effects.
Let \(v_t\) be the proportion of the susceptible population at time \(t\) who are removed through some means other than infection. These are individuals who have never been infected but may have been previously vaccinated, and their immunity is assumed to have developed at time \(t\).
epidemia requires \(v_t\) to be supplied by the user. Then (1.10) is replaced with \[\begin{equation} S_t = \left(S_{t-1} - i_t\right) \left(1 - v_t \right). \end{equation}\] Of course, \(v_t\) is a difficult quantity to estimate. It requires the user to estimate the time-lag for a jab to become effective, and to also adjust for potentially different efficacies of jabs and doses. Recognizing this, we allow the update \[\begin{equation} S_t = \left(S_{t-1} - i_t\right) \left(1 - v_t \xi \right), \end{equation}\] where \(\xi\) is a noise term that is assigned a prior distribution. \(\xi\) helps to account for potentially systematic biases in calculating vaccine efficacy.
Bellman, R., and T. E. Harris. 1948. “On the Theory of Age-Dependent Stochastic Branching Processes.” Proceedings of the National Academy of Sciences. https://doi.org/10.1073/pnas.34.12.601.
Bellman, Richard, and Theodore Harris. 1952. “On Age-Dependent Binary Branching Processes.” The Annals of Mathematics. https://doi.org/10.2307/1969779.
Bhatt, Samir, Neil Ferguson, Seth Flaxman, Axel Gandy, Swapnil Mishra, and James A Scott. 2020. “Semi-Mechanistic Bayesian Modeling of COVID-19 with Renewal Processes.” arXiv Preprint arXiv:2012.00394. https://arxiv.org/abs/2012.00394.
Cauchemez, Simon, Alain Jacques Valleron, Pierre Yves Boëlle, Antoine Flahault, and Neil M. Ferguson. 2008. “Estimating the impact of school closure on influenza transmission from Sentinel data.” Nature. https://doi.org/10.1038/nature06732.
Champredon, David, Jonathan Dushoff, and David J. D. Earn. 2018. “Equivalence of the Erlang-distributed SEIR epidemic model and the renewal equation.” SIAM Journal on Applied Mathematics. https://doi.org/10.1137/18M1186411.
Cori, Anne, Neil M. Ferguson, Christophe Fraser, and Simon Cauchemez. 2013. “A new framework and software to estimate time-varying reproduction numbers during epidemics.” American Journal of Epidemiology. https://doi.org/10.1093/aje/kwt133.
Flaxman, Seth, Swapnil Mishra, Axel Gandy, H Juliette T Unwin, Thomas A Mellan, Helen Coupland, Charles Whittaker, et al. 2020. “Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe.” Nature. https://doi.org/10.1038/s41586-020-2405-7.
Fraser, Christophe. 2007. “Estimating individual and household reproduction numbers in an emerging epidemic.” PLoS ONE. https://doi.org/10.1371/journal.pone.0000758.
Hoffman, Matthew D., and Andrew Gelman. 2014. “The no-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo.” Journal of Machine Learning Research 15.
Mishra, Swapnil, Tresnia Berah, Thomas A Mellan, H Juliette T Unwin, Michaela A Vollmer, Kris V Parag, Axel Gandy, Seth Flaxman, and Samir Bhatt. 2020. “On the derivation of the renewal equation from an age-dependent branching process: an epidemic modelling perspective.” arXiv Preprint arXiv:2006.16487.
Nouvellet, Pierre, Anne Cori, Tini Garske, Isobel M. Blake, Ilaria Dorigatti, Wes Hinsley, Thibaut Jombart, et al. 2018. “A simple approach to measure transmissibility and forecast incidence.” Epidemics. https://doi.org/10.1016/j.epidem.2017.02.012.