# Semi-Mechanistic Bayesian modeling of COVID-19 with Renewal Processes

Samir Bhatt  MRC Centre for Global Infectious Disease Analysis, Imperial
College, London, UK Neil Ferguson  MRC Centre for Global Infectious Disease
Analysis, Imperial College, London, UK Seth Flaxman  Department of
Mathematics, Imperial College, London, UK  Axel Gandy  Department of
Mathematics, Imperial College, London, UK  Swapnil Mishra  MRC Centre for
Global Infectious Disease Analysis, Imperial College, London, UK James A.
Scott  Department of Mathematics, Imperial College, London, UK

###### Abstract

We propose a general Bayesian approach to modeling epidemics such as COVID-19.
The approach grew out of specific analyses conducted during the pandemic, in
particular an analysis concerning the effects of non-pharmaceutical
interventions (NPIs) in reducing COVID-19 transmission in 11 European
countries. The model parameterizes the time varying reproduction number
$R_{t}$ through a regression framework in which covariates can e.g. be
governmental interventions or changes in mobility patterns. This allows a
joint fit across regions and partial pooling to share strength. This
innovation was critical to our timely estimates of the impact of lockdown and
other NPIs in the European epidemics, whose validity was borne out by the
subsequent course of the epidemic. Our framework provides a fully generative
model for latent infections and observations deriving from them, including
deaths, cases, hospitalizations, ICU admissions and seroprevalence surveys.
One issue surrounding our model’s use during the COVID-19 pandemic is the
confounded nature of NPIs and mobility. We use our framework to explore this
issue. We have open sourced an R package epidemia implementing our approach in
Stan. Versions of the model are used by New York State, Tennessee and Scotland
to estimate the current situation and make policy decisions.

##  1 Introduction

This article presents a general framework for semi-mechanistic Bayesian
modeling of infectious diseases using renewal processes. The term semi-
mechanistic relates to statistical estimation within some constrained
mechanism. Variants of this general model have been used in specific analyses
of Covid-19 [12, 13, 44, 26, 42, 31, 32, 39, 28], and continue to be used in
ongoing work to make policy decisions. The present article motivates and
discusses the key statistical and epidemiological features of this framework,
starting from a counting process setup. Various extensions of the basic model
are considered, including a latent infection process. We discuss limitations
and applications of the modelling framework to stimulate further research.

The model uses a flexible regression-based framework for parameterising
transmission and ascertainment rates. This allows the fitting of multilevel
models [16, 20, 23] for several regions simultaneously. Such partial pooling
of parameters has specific advantages in the context of infectious diseases.
Suppose we wish to estimate the effect of NPIs [10, 13] or mobility [1, 27] on
transmission rates. Estimating these effects separately in different regions
could lead to noisy estimates for at least two reasons. There is typically
little high quality data at the early stages of an epidemic. Such data is
generally correlated, reducing the information content that can be used to
infer such an effect. In addition, NPIs often occur in quick succession and
their effects are confounded [19, 46]. This is exacerbated by the random times
between infections (the generation distribution) and between infections and
observations, which smooths the observed data, making it more difficult to
attribute changes in transmission rates to a particular NPI. Alternatively,
one could pool the effect across all groups. This ignores group-level
variations and can lead to poor predictive performance, in particular
underestimating variance for previously unmodeled regions. One could augment
such a model with group-level indicators, but this results in a large number
of parameters, which are difficult to estimate and leads to overfitting with
classical estimation techniques. Partial pooling provides a natural solution
to this.

Sometimes the inferential goal is not to assess the effect of a covariate on
outcomes, but rather to infer transmission rates over time. Previous studies
have focused on estimating reproduction numbers from case data [11, 36, 5, 15,
21, 9], sometimes directly substituting observed case counts for the unknown
number of infected individuals [45]. However, the emergence of SARS-CoV-2 has
highlighted shortcomings of methods that rely on just case data. Limited
testing capacity at the early stages of the pandemic led to only a small
proportion of infections being detected and reported [24]. Those tested were
typically more likely to have been hospitalised or were at higher risk of
infection or death. In particular this proportion, referred to as the
infection ascertainment rate (IAR) is country-specific and likely to have
changed over time due to changes in testing policies and capacity. If
unaccounted for, it will lead to biases in the inferred transmission rates.

This highlights the need for more flexible observational models, whereby more
varied types of data can be incorporated, and their idiosyncrasies accounted
for. Daily death data has been used in [13] to recover reproduction numbers in
the early stages of the SARS-CoV-2 pandemic, and has been seen as more
reliable than case data. However, there have been clear variations in
definitions and reporting across time and countries. It is therefore important
to appropriately model noise within the the observational models. Our
framework allows for multiple types of data including deaths, cases,
hospitalizations, ICU admissions and the results of seroprevalence surveys.
This improves robustness of inferred parameters to biases in any one type of
data.

The model uses discrete renewal processes to propagate infections within
modeled populations. These have been used in a number of previous studies [14,
9, 30, 7], and are linked to other popular approaches to infectious disease
modeling. [8] show that the renewal equation leads to identical dynamics as
Erlang-Distributed Susceptible-Exposed-Infected-Recovered (SEIR) compartmental
models, when a particular form is used for the generation distribution. A
special case of this is the standard Susceptible-Infected-Recovered (SIR)
model [22]. The approach is also connected to counting processes including the
Hawkes process and the Bellman-Harris process [3, 4]. [3, 29] derive the
renewal equation as the expectation of an age-dependent branching process.
Age-dependence allows for more realistic dynamics than age-insensitive
processes, like the Galton-Watson process [2, 17]. More complex branching
processes such as the Crump-Mode-Jagers branching process could also be
considered. Hawkes processes are also related to renewal processes, with the
expectation of the Hawkes intensity function resulting in the renewal equation
[37].

We describe the general model in detail, and start by considering the bare-
bones version in Section 2. The motivation for the model lies in continuous-
time counting processes, and this connection is discussed in Section 3.
Sections 4 and 5 present the infection and observation processes in more
detail, and consider important extensions of the basic model. Section 6
considers how to use the framework for multilevel modeling. Section 7 compares
our approach to standard time series models, and outlines the key challenges
involved in modeling with our framework. Section 8 considers the specific
aspect of confounding and causality when estimating the effects of variables
on transmission rates. Section 9 has a brief discussion.

##  2 Model Overview

We now formulate a basic version of the model for one homogeneous population.
The same model can be used for multiple regions or groups jointly. Let
$R_{t}>0$ be the reproduction number at time $t>0$. This determines the rate
at which infections grow. Infections $i_{v},\ldots,i_{0}$ for some $v\leq 0$
are given a prior distribution. For $t>0$, we let new infections $i_{t}$ be
defined by

$\displaystyle i_{t}$ $\displaystyle=R_{t}\sum_{s<t}i_{s}g_{t-s},$ (1)

where the generation time, the lag between infections, is given through a
probability mass function $g$, i.e. $g_{t}\geq 0$ and
$\sum_{t=1}^{\infty}g_{t}=1$.

Observations occur at certain times $t>0$ . In general, there may be multiple
types; case and death counts for example. Each such type is driven by its own
time-varying ascertainment rate $\alpha_{t}>0$. The mean of the observations
at time $t$ is linked to past infections by

$\displaystyle y_{t}$ $\displaystyle=\alpha_{t}\sum_{s\leq t}i_{s}\pi_{t-s},$
(2)

where $\pi$ is a distribution for the lag between an infection and when it
gives rise to an observation. The sampling distribution of the observations
with these means is typically nonegative and discrete, and may depend on
auxiliary parameters. When multiple types are observed, we can superscript the
quantities as $y^{(l)}_{t},\alpha^{(l)}_{t}$ and $\pi^{(l)}$ and assign
independent sampling distributions for each type.

Transmission rates $R_{t}$ and ascertainment rates $\alpha_{t}$ can be modeled
flexibly using Bayesian regression models, and through sharing of parameters,
are the means through which we tie together multiple regions or groups using
multilevel modeling. One can, for example, model transmission rates as
depending on a binary covariate for an NPI, say full lockdown. The coefficient
for this can be partially pooled between these groups. The effect is to share
information between groups, while still permitting between group variation.

##  3 Motivation from continuous time

Our model can be motivated from a continuous time perspective as follows.
Infections give rise to additional infections in the future, referred to as
offspring. Letting $N^{I}(t)$ denote the number of infections occurring up to
time $t$, defined by its intensity

$\lambda(t)=R(t)\int_{s<t}g(t-s)N^{I}(ds),\quad t>0,$ (3)

where $g$ is the density of a probability distribution on $\mathbb{R}^{+}$
defining the time between infections, and where $\\{R(t):t>0\\}$ is a non-
negative stochastic process. The process can be initialised by assuming values
for $N^{I}(t)$ for $t$ in the seeding period $[v,0]$.

Equation (3) is similar to the Hawkes intensity, however the memory kernel $g$
is scaled by a time-specific factor $R(t)$. The integrand $g$ allows the
intensity to increase due to previous infection events, while $R(t)$ tempers
the intensity for other time-specific considerations. If $R(t^{\prime})=R(t)$
for all $t^{\prime}$ then Equation (3) reduces to a Hawkes process. Under this
assumption, since $g$ integrates to unity, the expected number of offspring is
simply $R(t)$, and so this is the instantaneous reproduction number or
alternatively the branching factor of the Hawkes process. The generation time,
defined as the time from an infection to a secondary infection, is distributed
according to $g$ and so $g$ is the generation distribution.

Observations are precipitated by past infections; a given infection may lead
to observation events in the future. Letting $N^{Y}(t)$ be the count of some
observation type over time defined by the intensity

$\lambda_{y}(t)=\alpha(t)\int_{s<t}\pi(t-s)N^{I}(ds),$ (4)

for $t>0$, where $\pi:\mathbb{R}^{+}\to\mathbb{R}^{+}$ is a function and
$\\{\alpha(t):t\geq 0\\}$ is a non-negative stochastic process. This is
similar to Equation (3), however the intensity increases due to past
infections, rather than past observations.

Consider the special case where $\pi$ is a probability density and where
$\alpha(t^{\prime})=\alpha(t)$ for all $t^{\prime}$. The average number of
observation events attributable to a single infection is then $\alpha(t)$, and
so this is an instantaneous ascertainment rate. $\pi$ is then interpreted as
the distribution for the time from an infection to an observation, and
therefore we call it the infection to observation distribution.

##  4 Infection Process

Starting from the continuous model, we now describe a discrete model, which
results in the formulation of Section 2. This discrete model is more amenable
to inference. Let $I_{t}$ be the number of new infections at time $t$; this is
the equivalent of $N^{I}_{t}-N^{I}_{t-1}$ in the continuous model. As basic
modelling block we use the following discrete version of (3):

$\operatorname{E}[I_{t}|R_{1:t},I_{v:t-1}]=R_{t}L_{t},$ (5)

where $L_{t}:=\sum_{s<t}I_{s}g_{t-s}$ is the case load or total infectiousness
by time $t>0$. Moreover, letting
$i_{t}:=\operatorname{E}[I_{t}|R_{1:t},I_{v:0}]$ and taking the conditional
expectation given reproduction numbers $R_{1:t}$ and seeded infections
$I_{v:0}$ on both sides of (5) gives

$i_{t}=R_{t}\operatorname{E}[L_{t}|R_{1:t},I_{v:0}]=R_{t}\sum_{s<t}\operatorname{E}[I_{s}|R_{1:s},I_{v:0}]g_{t-s}=R_{t}\sum_{s<t}i_{s}g_{t-s},$

which is Equation (1). This is a discrete renewal equation, which can
alternatively be interpreted as an $\text{AR}(t)$-process with known
coefficients $g_{k}$. From this point of view, the basic model in Section 2
uses $i_{t}$ as synonymous with actual infections. Since infections are simply
a deterministic function of other parameters, there is no need to treat them
as unknown latent parameters to sample. This can lead to lower sampling times
and faster convergence.

###  4.1 Modeling Latent Infections

The model of Section 2 can be extended by replacing each $i_{t}$ with the
actual infections from the counting process $I_{t}$, and then assigning a
prior to $I_{t}$. Although sampling can be slower, this has certain
advantages. When past infection counts are low, significant variance in the
offspring distribution can imply that new infections $I_{t}$ has high
variance. This is not explicitly accounted for in the basic model. In
addition, this approach cleanly separates infections and observations; the
latter being modeled conditional on actual infections. The sampling
distribution can then focus on idiosyncrasies relating to the observation
process.

We assign a prior to $I_{t}$ conditional on previous infections and current
transmission $R_{t}$. The expected value for this is given by Equation (5).
Appendix 10.1 shows that assuming the variance of the prior to be a constant
proportion $d$ of this mean is equivalent to letting $d$ be the coefficient of
dispersion for the offspring distribution. $d>1$ implies overdispersion, and
can be used to account for super-spreading events, which has been shown to be
an important aspect for modeling Covid-19. The parameter $d$ can be assigned a
prior.

Any two parameter family can be used to match these first two moments. Letting
this be continuous rather than discrete allows inference to proceed using
Hamiltonian Monte Carlo, whereby new values for $I_{t}$ are proposed
simultaneously with all other parameters. Possible candidates include log-
normal, gamma and the Weibull distributions. If an explicit distribution for
the offspring distribution is desired, one can show that assuming a Gamma
distribution with rate $\lambda$ for this results in a Gamma distribution for
$I_{t}$ with rate $\lambda$. The coefficient of dispersion is then simply
$D=\lambda^{-1}$.

###  4.2 Population Adjustments

If $R_{t}$ remains above unity over time, infections grow exponentially
without limit. In practice, infections should be bounded from above by
$S_{0}$, the initial susceptible population. All else being equal,
transmission rates are expected to fall as the susceptible population is
diminished.

Consider first the model using $I_{t}$, which was described in Section 4.1.
Equation 5 can be replaced with

$\operatorname{E}[I_{t}|R_{1:t},I_{v:t-1}]=(S_{0}-I_{t-1})\left(1-\exp\left(-\frac{R_{u,t}L_{t}}{S_{0}}\right)\right),$
(6)

where $R_{u,t}$ is an unadjusted reproduction number, which does not account
for the susceptible population. This satisfies intuitive properties. As the
unadjusted expected infections $R_{u,t}L_{t}$ approaches infinity, the
adjusted expected value approaches the remaining susceptible population. The
motivation for and derivation of Equation (6) is provided in Appendix 10.2. In
short, this is the solution to a continuous time model whose intensity is a
simplification of Equation (3). We must also ensure that the distribution of
$I_{t}$ cannot put positive mass above $S_{0}-I_{t-1}$. A simple solution is
to use truncated distributions. Of course, this adjusts the mean value from
Equation (6), however this is unlikely to be significant unless the
susceptible population is close to depletion.

In the basic model, one can apply the adjustment to $i_{t}$ by replacing
$L_{t}$ in Equation (6) with

$\operatorname{E}(L_{t}|R_{1:t},I_{v:0})=\sum_{s<t}i_{s}g_{t-s}.$ (7)

##  5 Observations

Observations are modeled in discrete time, analogous to how we treated
infections in Section 4. Letting $\pi:\mathbb{N}\to\mathbb{R}^{+}$ and
$Y_{t}:=N^{Y}_{t}-N^{Y}_{t-1}$, the discrete analogue to Equation (4) is

$\operatorname{E}[Y_{t}|\alpha_{t},I_{v:t}]=\alpha_{t}\sum_{s\leq
t}I_{s}\pi_{t-s}.$ (8)

Taking the expected value of the above given seeded infections, transmission
rates and the current ascertainment rate gives

$\operatorname{E}[Y_{t}|\alpha_{t},R_{1:t},I_{v:0}]=\alpha_{t}\sum_{s\leq
t}i_{s}\pi_{t-s},$ (9)

which is recognisable as Equation (2). Thus we have two possible expressions
for the mean of $Y_{t}$, one given actual infections, and the other given
expected infections $i_{t}$. The basic model of Section 2 uses the latter,
while the extension in Section 4.1 uses the former.

We assume that $Y_{t}\sim\mathcal{F}\left(y_{t},\phi\right)$, where
$\mathcal{F}$ is a non-negative discrete family parameterised by its mean
$y_{t}$ and potentially an auxiliary parameter $\phi$. This could be a Poisson
distribution, where there is no auxiliary parameter. Using a quasi-Poisson or
negative binomial instead allows for overdispersion. This can be useful to
capture, for example, day-to-day variation in ascertainment rates when
infection counts are low. The mean $y_{t}$ can be taken to be either (8) or
(9), the latter being used in the basic version of the model. Hidden in this
formulation is the assumption that each $Y_{t}$ is conditionally independent
given $y_{t}$. Using multiple observation series $Y^{(l)}_{t}$ can help to
improve the model inferences and identifiability of certain parameters. We
simply assume that each such series is conditionally independent given the
underlying infection process.

##  6 Multilevel Models

Transmission rates can be modeled quite generally within the framework. If the
aim is simply to estimate transmission in a single region over time, one
approach could be to let $R_{t}=g^{-1}(\gamma_{t})$, where $g$ is a link
function and $\gamma_{t}$ is some autocorrelation process, for example a
random walk. Suppose, however, that transmission is modeled in $M$ regions and
the goal is to estimate the effect of a series of NPIs on transmission.
Letting $R^{(m)}_{t}$ denote transmission in region $m$ at time $t$, we could
let

$R^{(m)}_{t}=g^{-1}\left(\beta^{(m)}_{0}+\sum_{l=1}^{p}x^{(m)}_{t}\beta^{(m)}_{k}\right),$
(10)

where $x^{(m)}_{t}$ are binary encodings of NPIs, and $\beta^{(m)}_{0}$ and
$\beta_{k}^{(m)}$ are region-specific intercepts and effects respectively. The
intercepts are used to allow regions to have their own baseline transmission
rates. Collecting these group specific parameters into $\beta^{(m)}$, we can
partially pool them by letting
$\beta^{(m)}\sim\mathcal{N}\left(0,\Sigma\right)$, for each group $m$, and the
assigning a prior to the covariance matrix $\Sigma$. This could be an inverse-
Wishart prior, or alternatively, $\Sigma$ can be decomposed into variances and
a correlation matrix, which are each given separate priors [41].

One possible option for $g$ is is the log-link. This provides easily
interpretable effect sizes; a one unit change in a covariate multiplies
transmission by a constant factor. However, this can lead to prior mass on
unreasonably high transmission rates. With this in mind, an alternative is to
use a generalisation of the logit link for which

$g^{-1}(x)=\frac{K}{1+e^{-x}},$ (11)

and where $K$ is the maximum possible value for transmission rates. This
serves a similar purpose to the carrying capacity in a logistic growth model.

The ascertainment rate $\alpha_{t}$ can also be modeled with similar
considerations to the above. This flexibility is useful, particularly because
these quantities are likely to change as an epidemic progresses. This has been
clearly seen during the Covid-19 epidemic, where the infection ascertainment
rate may have increased over time due to increased testing capacity and
improved track and trace systems. Multilevel models can in theory also be
specified through $\alpha_{t}$.

##  7 Forecasting, epidemiological constants, and seeding

A key benefit of using a semi-mechanistic approach is that forecasts are
constrained by plausible epidemiological mechanisms. For example, in the
absence of any further interventions or behavioural changes, and looking at a
medium term forecast of just incidence (daily new cases/infections), a
traditional time series forecasting approach may predict a constant function
based on observing broadly constant incidence, but semi-mechanistic approach
would expect a monotonic decrease based on a constant rate of transmission and
the effect of herd immunity. The performance of epidemiologically constrained
models is generally good [6]; this is perhaps not surprising as examining the
discrete renewal equation shows that these models correspond to
autoregressive(n) filters with a convex combination of coefficients specified
by the generation interval. However, similar to financial forecasting, the
predictive capability of epidemic models are likely to be better interpreted
as scenarios rather than actual predictions due to the rapidly adaptive
landscape of policy.

A second benefit of epidemic models is to provide a plausible mechanism to
explain (non causally) the changes observed in noisy data. For example, in
estimating the effect of an intervention on observed death data, we need
consider what that intervention affects i.e the rate of transmission or
$R_{t}$. As we have described above, we can connect the rate of transmission
to latent infections to deaths such that we have an epidemiologically
motivated mechanism. While we can statistically estimate parameters for how
the intervention affects $R_{t}$, certain important parameters will be
entirely unidentifiable and need to be fixed as constants or with very tight
priors. For example, to reliably estimate the number of infections, an
infection fatality rate needs to be chosen. A failure to choose an appropriate
infection fatality rate can result in a bimodal posterior where changes can
either be attributed to herd immunity or to interventions. From a statistical
perspective, is it difficult to disentanlge which mode of the posterior best
represents reality, and hence it is sensible to first estimate a plausible
infection fatality rate and then use this within the semi-mechanistic model. A
second example is the onset of symptoms-to-death distribution. Given the lag
between transmission, infection and deaths, the effect of an intervention is
dependent on the onset of symptoms-to-death distribution.

Infection seeding is a fundamentally challenging aspect of epidemic modelling.
Estimating the initial effect of seeding is crucial to understanding a
baseline rate of growth ($R_{0}$) from which behaviours and interventions can
modify. This seeding is heavily confounded by importation and under
ascertainment. Both these factors can influence estimates of the initial
growth rates, and this in turn can affect the impact of changes in
transmission as time progresses. We have proposed heuristic approaches to
mitigate issues with early seeing, but principled statistical approaches need
to be developed. In particular, Bayesian pair plots show strong correlation
between seeding parameters and $R_{0}$, which can potentially lead again to a
bimodal posterior where initial growth dynamics can be explained through
$R_{0}$ or via initial seed infections.

##  8 Confounding and Causality: Estimating the Effect of Interventions

Section 6 showed that changes in transmission rates over time can be explained
by parameterising these rates in terms of predictor variables, such as NPIs
and mobility. Clarifying the effects of interventions on transmission is
important, if only because of their clear economic and human consequences.
This is however a significant challenge because the effects are potentially
confounded with unobserved behavioural changes, and they are hard to identify.
Identification is difficult because interventions are highly correlated when
they occur in quick succession. Moreover, the random time between an infection
and its recording as a case or death leads to observations being less
informative about the effect of any particular intervention.

[13] estimated the effectiveness of NPIs across 11 European countries, and
used partial pooling of effect sizes to address the identification problem. At
that time, little data existed other than information on deaths and the timing
of interventions. NPIs, which were coded as a binary set of mandatory
government measures (e.g. school closures, ban on public events, lockdown),
could not fully explain the patterns seen in some countries (e.g. Sweden), and
especially at the subnational level. Mobility data became available in April
and was used to model the epidemic in Italy, Brazil and the USA [44, 26, 42].
Such data is useful as it may help account for behavioral changes that
confound the effects of NPIs. However since mobility affects transmission, is
linked to the introduction of NPIs and potentially also to voluntary
behavioural measures, we expect it to be a confounder.

Section 8.2 extends the model in [13] to investigate further this issue of
confounding, and models both NPIs and mobility jointly. This is in keeping
with standard practice in regression/ANOVA: expanding a model to take into
account more explanatory variables. Nonetheless, NPIs may partially affect
transmission via a path through mobility. A joint model of mobility and NPIs
does not account for this. Therefore, in Section 8.3 we take a first step in
assessing causal considerations through a simple mediation analysis. We begin
however by exploring the relationship between interventions and mobility.

###  8.1 Interventions and Mobility

Regressing average mobility on NPIs in a Bayesian linear model (no intercept
or partial pooling) we find a correlation of over $85\%$ with a mean absolute
error of $0.1\%$. Given mobility generally ranges from -1 to 1, this is a good
overall fit. Figure 1a shows that these fits visually correspond well with
changes in average mobility. One could conjecture that mobility and NPIs are
lagged, but lagging NPI dates either forwards or backwards in time does not
result in a better fitting model (see Figure 1b). Indeed Figure 1b does
support a hypothesis that the timing of NPIs and changes in mobility are
strongly linked. The coefficient sizes from this regression are entirely
consistent with [13] finding that the NPI with the largest effect size is
lockdown (see Figure 1c). This simple analysis does not model transmission,
but does provide strong evidence that mobility and NPIs do not provide
conflicting narratives. We note, to perform this regression as fairly as
possible we used a hierarchical shrinkage prior [35] that performs both
shrinkage and variable selection.

(a)

(b)

(c)

Figure 1: Simple regression of mobility against NPIs. Figure (a) shows the fit
for the United Kingdom with mobility in black and the fit for NPIs in red. (b)
Shows the effect on mean absolute error from lagging the NPIs. (c) Shows the
coefficient effect sizes from the regression. Y axis are NPIs and X axis the
regression effect sizes

###  8.2 Controlling for Mobility

Section 8.1 found a large correlation between interventions and mobility,
demonstrating that mobility is a potential confounder. Here we control for
this by jointly modeling NPIs and mobility. This is done using the same 11
European countries, sets of NPIs and death data as used in [13].

A two-stage approach [18] is used, whereby $R_{t}$ is first estimated using a
daily random walk. In the second stage this is regressed on NPIs and mobility.
The random walk can in theory select any arbitrary function of $R_{t}$ that
best describes the data without any prior information about which
interventions happened when or how well they worked. Given these estimates of
$R_{t}$ for all 11 European countries, we run a simple partial pooling model
to see if interventions and/or mobility can reproduce the trends in $R_{t}$.
The model used is a linear regression with country level intercepts (to
account for variation in $R_{0}$), and both joint and country specific effect
sizes for interventions/mobility. As with the earlier analysis we use a
hierarchical shrinkage prior on the coefficients [35].

Three variations of the model are considered: NPIs only (NPI_only), mobility
only (Mobility_only), and both NPIs and mobility (NPI+Mobility). MCMC
convergence diagnostics in all cases did not indicate problems. We found the
best fitting model to be NPI+Mobility. Relative to this model the expected log
posterior difference ($\pm$ standard error) in WAIC of the model with only
NPIs is $-5.2\pm 4$, ad $-565.6\pm 49.2$ with only mobility. Therefore, in
fits to the estimated $R_{t}$, the model with mobility alone is substantially
worse than the models with NPIs. Controlling for mobility does not appear to
significantly change the estimated effects of NPIs. As in [13], the largest
effect size is attributed to lockdown, as seen in Figure 2. This is true with
and without the inclusion of the mobility variable. This analysis could be
improved using Bayesian leave-one-out cross-validation [43] to account for the
time series nature of the data.

(a)

(b)

(c)

Figure 2: Regression of NPIs and/or mobility against nonparametric $R_{t}$.
Figure (a) shows the fit for NPIs only. (b) for mobility only, (c) NPIs and
mobility. Mobility only was not _significantly_ preferred by WAIC. Y axis are
covariates and X axis the regression effect sizes

An advantage of the two-stage approach is that it is scalable to a large
number of regions. $R_{t}$ can be estimated in each region in parallel using
separate models. Partial pooling can still be leveraged to estimate effects in
the second stage. Once $R_{t}$ has been estimated, any number of interesting
statistical analyses can be conducted. Nonetheless the estimated $R_{t}$ is
not entirely non-parametric; it is clearly influenced by its random walk form
in the first stage. This analysis could be extended by considering a range of
alternative priors for $R_{t}$. More importantly however, this approach has
not considered causal relationships between NPIs and mobility. This is the
focus of the next example.

###  8.3 Causal Mediation

The effect of interventions in the previous analysis holds mobility constant.
However, we intuitively expect that part of these effects occur indirectly
through changeing mobility. We can hypothesise that changes in mobility are
both an effect of NPIs and a cause of reductions in transmission. Causal
mediation analysis provides a means to disentangle the total effect of a
variable into a direct and indirect effect. The indirect effect occurs via
some mediator, which in this case is hypothesised to be mobility. Further
information about causal mediation can be found here [33, 34].

Only lockdown is considered here because performing causal mediation with all
NPIs is challenging and lockdown is consistently the NPI with the largest
effect size in Section 8.1, Section 8.2 and in [13]. Briefly, to perform
causal mediation we consider two transmission models

$\displaystyle R_{t}^{(m)}$
$\displaystyle=\tilde{R}^{1}_{m}\exp\left(\left(\beta^{1}_{1}+\beta^{1}_{1,m}\right)L_{t,m}+\epsilon^{1}_{t,m}\right),$
(12) $\displaystyle R_{t}^{(m)}$
$\displaystyle=\tilde{R}^{2}_{m}\exp\left(\left(\beta^{2}_{1}+\beta^{2}_{1,m}\right)L_{t,m}+\left(\beta^{2}_{2}+\beta^{2}_{2,m}\right)M_{t,m}+\epsilon^{2}_{t,m}\right),$
(13)

where $L_{t,m}$ is a binary indicator for lockdown and $M_{t,m}$ is mobility
in country $m$ respectively. $\tilde{R}^{i}_{m}$ and $\epsilon^{i}_{t,m}$ are
country specific parameters modeling baseline transmission and a weekly random
walk respectively. All other aspects of both models are the same as in [13].
Model (12) includes effects for lockdown, while (13) additionally considers
mobility. $\beta^{1}_{1}$ is the total effect for lockdown, while
$\beta^{2}_{1}$ is the partial effect when controlling for mobility. The
mediated effect is therefore $\beta^{1}_{1}-\beta^{2}_{1}$. This quantifies
the effect of lockdown via the path through mobility. We find this mediated
effect reduces $R_{t}$ by 18.3% with a 95% credible interval of
$[12.2\%,44.4\%]$. The posterior probability of the effect being greater than
$0$ is 89.6% . Individual coefficients are shown in Figure 3.

(a)

(b)

Figure 3: Mediation analysis (a) Lockdown, (b) Lockdown and Mobility

These mediation results suggest a causal link between lockdown and mobility
that eventually leads to reduced transmission rates. They also suggest that
the mediated effect is far less than the total effect of lockdown, suggesting
lockdown will have other causal pathways. Of course, mobility is also mediated
through other pathways, and a principled causal analysis is out of the scope
of this article. The exclusion of other NPIs may introduce omitted variable
bias. Nonetheless, this simple analysis with lockdown adds support to causal
considerations.

##  9 Discussion

This article has discussed a class of statistical models for epidemics such as
Covid-19 which are able to capture key epidemiological mechanisms. The model
has appeared in various forms for specific analyses during the Covid-19 crisis
and, at the time of writing, continues to be used to inform public policy. By
presenting it in a general form and discussing key modelling difficulties we
hope to stimulate discussion around it. One key difficulty within the
framework is dealing with confounded variables, particularly those used to
explain changes in transmission during the early stages of an epidemic. The
analyses in Section 8 make a first step in dealing with these, and support the
central finding of [13]: that lockdown and other NPIs together served to
control the first wave of the epidemic in 11 European countries.

A number of model enhancements have not been discussed here. These include
explicitly accounting for importations, and allowing for uncertainty in the
generation and infection to observation distributions. The model can be fit
using Stan [38], but the adaptive Hamiltonian Monte Carlo used often has
difficulty converging when latent infections are modeled directly, or when
multiple regions are jointly modeled. We conjecture that convergence may be
improved by carefully choosing initial parameters for the sampler. Future
research could explore whether alternative samplers can be developed to fit
these models more pragmatically.

## References

  * [1] Hamada S. Badr et al. “Association between mobility patterns and COVID-19 transmission in the USA: a mathematical modelling study” In _The Lancet Infectious Diseases_ , 2020 DOI: 10.1016/S1473-3099(20)30553-3
  * [2] Robert Bartoszynski “Branching Processes and the Theory of Epidemics” In _Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability_ , 1967
  * [3] R. Bellman and T. E. Harris “On the Theory of Age-Dependent Stochastic Branching Processes” In _Proceedings of the National Academy of Sciences_ , 1948 DOI: 10.1073/pnas.34.12.601
  * [4] Richard Bellman and Theodore Harris “On Age-Dependent Binary Branching Processes” In _The Annals of Mathematics_ , 1952 DOI: 10.2307/1969779
  * [5] Luís M.A. Bettencourt and Ruy M. Ribeiro “Real time bayesian estimation of the epidemic potential of emerging infectious diseases” In _PLoS ONE_ , 2008 DOI: 10.1371/journal.pone.0002185
  * [6] Cristina Carias et al. “Forecasting the 2014 West African Ebola Outbreak” In _Epidemiologic Reviews_ 41.1, 2019, pp. 34–50 DOI: 10.1093/epirev/mxz013
  * [7] Simon Cauchemez et al. “Estimating the impact of school closure on influenza transmission from Sentinel data” In _Nature_ , 2008 DOI: 10.1038/nature06732
  * [8] David Champredon, Jonathan Dushoff and David J.D. Earn “Equivalence of the Erlang-distributed SEIR epidemic model and the renewal equation” In _SIAM Journal on Applied Mathematics_ , 2018 DOI: 10.1137/18M1186411
  * [9] Anne Cori, Neil M. Ferguson, Christophe Fraser and Simon Cauchemez “A new framework and software to estimate time-varying reproduction numbers during epidemics” In _American Journal of Epidemiology_ , 2013 DOI: 10.1093/aje/kwt133
  * [10] Benjamin J. Cowling et al. “Impact assessment of non-pharmaceutical interventions against coronavirus disease 2019 and influenza in Hong Kong: an observational study” In _The Lancet Public Health_ , 2020 DOI: 10.1016/S2468-2667(20)30090-6
  * [11] N. M. Ferguson, C. A. Donnelly and R. M. Anderson “Transmission intensity and impact of control policies on the foot and mouth epidemic in Great Britain” In _Nature_ , 2001 DOI: 10.1038/35097116
  * [12] S Flaxman et al. “Report 13: Estimating the number of infections and the impact of non-pharmaceutical interventions on COVID-19 in 11 European countries” In _Spiral_ Imperial College London], 2020 URL: https://doi.org/10.25561/77731
  * [13] Seth Flaxman et al. “Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe” In _Nature_ , 2020 DOI: 10.1038/s41586-020-2405-7
  * [14] Christophe Fraser “Estimating individual and household reproduction numbers in an emerging epidemic” In _PLoS ONE_ , 2007 DOI: 10.1371/journal.pone.0000758
  * [15] Christophe Fraser et al. “Pandemic potential of a strain of influenza A (H1N1): Early findings” In _Science_ , 2009 DOI: 10.1126/science.1176062
  * [16] Andrew Gelman and Jennifer Hill “Data Analysis Using Regression and Multilevel/Hierarchical Models” In _Data Analysis Using Regression and Multilevel/Hierarchical Models_ , 2006 DOI: 10.1017/cbo9780511790942
  * [17] Wayne Getz and James Lloyd-Smith “Basic methods for modeling the invasion and spread of contagious diseases”, 2006 DOI: 10.1090/dimacs/071/05
  * [18] Nils Haug et al. “Ranking the effectiveness of worldwide COVID-19 government interventions” In _Nat. Hum. Behav._ 4.12 Nature Publishing Group, 2020, pp. 1303–1312 DOI: 10.1038/s41562-020-01009-0
  * [19] SARS Expert Committee HKSAR “Chronology of the SARS epidemic in Hong Kong.”, 2003
  * [20] Joop J Hox, Mirjam Moerbeek and Rens Schoot “Multilevel analysis: Techniques and applications” Routledge, 2010
  * [21] Heath A. Kelly et al. “Pandemic (H1N1) 2009 influenza community transmission was established in one Australian state when the virus was first identified in North America” In _PLoS ONE_ , 2010 DOI: 10.1371/journal.pone.0011341
  * [22] A. G. Kermack, William Ogilvy and McKendrick “A contribution to the mathematical theory of epidemics” In _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ , 1927 DOI: 10.1098/rspa.1927.0118
  * [23] Ita Kreft and Jan Leeuw “Introducing Multilevel Modeling” In _Introducing Multilevel Modeling_ , 2011 DOI: 10.4135/9781849209366
  * [24] Ruiyun Li et al. “Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)” In _Science_ , 2020 DOI: 10.1126/science.abb3221
  * [25] Robert M. May “Simple mathematical models with very complicated dynamics” In _Nature_ , 1976 DOI: 10.1038/261459a0
  * [26] Thomas A Mellan et al. “Subnational analysis of the COVID-19 epidemic in Brazil” In _medRxiv_ Cold Spring Harbor Laboratory Press, 2020 DOI: 10.1101/2020.05.09.20096701
  * [27] Andrew C Miller et al. “Mobility trends provide a leading indicator of changes in SARS-CoV-2 transmission” In _medRxiv_ Cold Spring Harbor Laboratory Press, 2020 DOI: 10.1101/2020.05.07.20094441
  * [28] Swapnil Mishra et al. “A COVID-19 Model for Local Authorities of the United Kingdom” In _medRxiv_ , 2020 DOI: 10.1101/2020.11.24.20236661
  * [29] Swapnil Mishra et al. “On the derivation of the renewal equation from an age-dependent branching process: an epidemic modelling perspective”, 2020 arXiv:2006.16487 [q-bio.PE]
  * [30] Pierre Nouvellet et al. “A simple approach to measure transmissibility and forecast incidence” In _Epidemics_ , 2018 DOI: 10.1016/j.epidem.2017.02.012
  * [31] NYS Press Office “Amid Ongoing COVID-19 Pandemic, Governor Cuomo Announces State is Bringing in International Experts to Help Advise the State’s Reopening Plan”, 2020 URL: https://www.governor.ny.gov/news/amid-ongoing-covid-19-pandemic-governor-cuomo-announces-state-bringing-international-experts
  * [32] Andrew M Olney et al. “Estimating the Effect of Social Distancing Interventions on COVID-19 in the United States” In _medRxiv_ Cold Spring Harbor Laboratory Press, 2020
  * [33] Judea Pearl “Causal inference in statistics: An overview” In _Stat. Surv._ 3 The author, under a Creative Commons Attribution License, 2009, pp. 96–146
  * [34] Judea Pearl “The causal mediation formula—a guide to the assessment of pathways and mechanisms” In _Prev. Sci._ 13.4 Springer, 2012, pp. 426–436
  * [35] Juho Piironen and Aki Vehtari “Sparsity information and regularization in the horseshoe and other shrinkage priors” In _Electron. J. Stat._ 11.2 The Institute of Mathematical Statisticsthe Bernoulli Society, 2017, pp. 5018–5051
  * [36] Steven Riley et al. “Transmission dynamics of the etiological agent of SARS in Hong Kong: Impact of public health interventions” In _Science_ , 2003 DOI: 10.1126/science.1086478
  * [37] Marian-Andrei Rizoiu et al. “Expecting to Be HIP: Hawkes Intensity Processes for Social Media Popularity” In _Proceedings of the 26th International Conference on World Wide Web_ , WWW ’17 Perth, Australia: International World Wide Web Conferences Steering Committee, 2017, pp. 735–744 DOI: 10.1145/3038912.3052650
  * [38] Stan Development Team “The {Stan} Core Library”, 2018 URL: http://mc-stan.org/
  * [39] The Scottish Government “Coronavirus (COVID-19): modelling the epidemic”, 2020 URL: https://www.gov.scot/collections/coronavirus-covid-19-modelling-the-epidemic/
  * [40] W. A. Thompson, Howard M. Taylor and Samuel Karlin “An Introduction to Stochastic Modeling.” In _Journal of the American Statistical Association_ , 1985 DOI: 10.2307/2287941
  * [41] Tomoki Tokuda et al. “Visualizing distributions of covariance matrices” In _Columbia Univ., New York, USA, Tech. Rep_ Citeseer, 2011, pp. 18–18
  * [42] H. Juliette T. Unwin et al. “State-level tracking of COVID-19 in the United States” In _Nature Communications (in Press)_ , 2020
  * [43] Aki Vehtari, Andrew Gelman and Jonah Gabry “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC” In _Stat. Comput._ 27.5, 2017, pp. 1413–1432
  * [44] Michaela A. C. Vollmer et al. “A sub-national analysis of the rate of transmission of COVID-19 in Italy” In _medRxiv_ Cold Spring Harbor Laboratory Press, 2020 DOI: 10.1101/2020.05.05.20089359
  * [45] Jacco Wallinga and Peter Teunis “Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures” In _American Journal of Epidemiology_ , 2004 DOI: 10.1093/aje/kwh255
  * [46] WHO “SARS: chronology of a serial killer.”, 2003

##  10 Appendix

###  10.1 Offspring Dispersion

Define the offspring distribution of any given infection to be the
distribution of the random number of offspring attributable to that infection.
We show that assuming the variance of these distributions are a constant
proportion of the mean implies, under suitable independence assumptions, the
same result for new infections $I_{t}$ for all time points.

Assume some ordering over infections at each period, and let $O_{t}^{(i)}$
denote the number of offspring of the $i$th infection at time $t$. This can be
decomposed as

$O_{t}^{(i)}=\sum_{s=t+1}^{\infty}O_{ts}^{(i)},$ (14)

where $O_{ts}^{(i)}$ are the number of offspring of $i$ birthed at time $s$.
The branching process behind Equation (5) implies that $O_{ts}^{(i)}$ has mean
$R_{s}g_{s-t}$. Assume that $\\{O^{(i)}_{ts}:s\geq t\\}$ are mutually
independent and have variance which is a fixed proportion $d$ of the mean. By
Equation (14), this implies the same variance relationship for $O_{t}^{(i)}$.
In particular, if $R_{s}=R_{t}$ for $s>t$ then $O_{t}^{(i)}$ has mean $R_{t}$
and variance $dR_{t}$. New infections at time $t$ can be expressed as

$I_{t}=\sum_{s=1}^{t-1}\sum_{i=1}^{I_{s}}O_{st}^{(i)}.$ (15)

Assume that all $O_{st}^{(i)}$ appearing in Equation (15) are mutually
independent conditional on everything occurring up to time $t-1$, the result
clearly follows by taking the variance of both sides of Equation 15 given
$R_{t}$ and $I_{v:t-1}$.

###  10.2 Population Adjustment

Here we motivate Equation (6), which is used to adjust transmission rates for
the size of the infectable population. The most obvious starting point for
such an adjustment would be to let

$\operatorname{E}[I_{t}|R_{t},I_{v:t-1}]=\left(\frac{S_{0}-I_{t-1}}{S_{0}}\right)R_{u,t}L_{t},$
(16)

where $R_{u_{t}}$ is defined as in Section 4.2. This is similar in form to a
discrete logistic growth model. Such models are well known as examples of
simple models that exhibit chaotic dynamics [25]. In particular, it is
possible that the expected value on the left hand side exceeds the remaining
susceptible population. Intuitively, this issue occurs because multiple
infections can occur simultaneously in the discrete model. We therefore
propose solving this by using a population adjustment motivated by the
solution to a continuous time model whose intensity is a simplification of
Equation (3).

Suppose we observe $I_{v:t-1}$ and current transmission $R_{t}$. We evolve
infections from time $t-1$ to $t$ continuously, and hence avoid overshooting.
Define a continuous time counting $\tilde{I}(s)$ process starting at time
$t-1$ by the intensity

$\tilde{\lambda}(s)=\left(\frac{S_{0}-\tilde{I}(s)}{S_{0}}\right)R_{u,t}L_{t},$
(17)

for $s\geq t-1$, and with initial condition $\tilde{I}(t-1)=I_{t-1}$.
Supplementary 11.1 shows that

$\operatorname{E}[\tilde{I}(t)]=I_{t-1}+(S_{0}-I_{t-1})\left(1-\exp\left(-\frac{R_{u,t}L_{t}}{S_{0}}\right)\right),$
(18)

which is the motivation for Equation (6).

##  11 Online Supplement

###  11.1 Proof of Equation (18)

Without loss of generality, we prove the result for time $t=1$. The argument
remains the same for all $t>1$.

From [40, Lemma 5.5], we have

$\operatorname{E}[\tilde{I}(s)]=\tilde{I}(0)+\int^{s}_{0}\operatorname{E}[\tilde{\lambda}(l)]dl\text{
for }s\geq 0.$

The following lemma derives an expression for the the expected intensity on
the right hand side.

######  Lemma 11.1.

The expected intensity takes the form

$E[\tilde{\lambda}(s)]=\tilde{\lambda}(0)\exp\left(-\frac{R_{u,1}L_{1}}{S_{0}}s\right),$

for all $s\geq 0$.

###### Proof of Lemma 11.1.

Fix $s\geq 0$, some small $\Delta>0$ and let
$h(s):=\operatorname{E}[\tilde{\lambda}(s)]$. We have from Equation (17) that

$h(s+\Delta)=\left(\frac{S_{0}-\operatorname{E}[\tilde{I}(s+\Delta)]}{S_{0}}\right)R_{u,1}L_{1}.$
(19)

We can write

$\operatorname{E}[\tilde{I}(s+\Delta)|\tilde{\lambda}(s)]=\operatorname{E}[\tilde{I}(s)|\tilde{\lambda}(s)]+\tilde{\lambda}(s)\Delta+\mathcal{O}(\Delta),$

and taking expectations on both sides,

$\operatorname{E}[\tilde{I}(s+\Delta)]=\operatorname{E}[\tilde{I}(s)]+h(s)\Delta+\mathcal{O}(\Delta).$

Substituting this into (19) and rearranging gives

$\begin{split}h(s+\Delta)&=\left(\frac{S_{0}-\operatorname{E}[\tilde{I}(s+\Delta)]}{S_{0}}\right)R_{u,1}L_{1}-\frac{R_{u,1}L_{1}}{S_{0}}\left(h(s)\Delta+\mathcal{O}(\Delta)\right),\\\
&=h(s)-\frac{R_{u,1}L_{1}}{S_{0}}\left(h(s)\Delta+\mathcal{O}(\Delta)\right).\end{split}$

Rearranging gives

$\frac{h(s+\Delta)-h(s)}{\Delta}=-\frac{R_{u,1}L_{1}}{S_{0}}\left(h(s)+\frac{\mathcal{O}(\Delta)}{\Delta}\right).$

Taking the limit as $\Delta\to 0$ and rearranging gives the differential
equation

$\frac{h^{\prime}(s)}{h(s)}=-\frac{R_{u,1}L_{1}}{S_{0}}.$

Integrating both sides gives

$\log(h(s))=-\frac{R_{u,1}L_{1}}{S_{0}}s+C.$

Using that $h(0)=\tilde{\lambda}(0)$ gives the constant
$C=\log(\tilde{\lambda}(0))$. Plugging in yields the required result. ∎

Hence,

$\begin{split}\operatorname{E}[\tilde{I}(s)]&=I_{0}+\tilde{\lambda}(0)\int_{0}^{s}\exp\left(-\frac{R_{u,1}L_{1}}{S_{0}}l\right)dl\\\
&=I_{0}+\tilde{\lambda}(0)\frac{S_{0}}{R_{u,1}L_{1}}\left(1-\exp\left(-\frac{R_{u,1}L_{1}}{S_{0}}s\right)\right)\\\
&=I_{0}+(S_{0}-\tilde{I}(s))\left(1-\exp\left(-\frac{R_{u,1}L_{1}}{S_{0}}s\right)\right).\end{split}$

Letting $s=1$ gives the required result.