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.