# Bayesian inference for high-dimensional discrete-time epidemic models:
spatial dynamics of the UK COVID-19 outbreak

Alison C Hale  Department of Mathematics and Statistics, Lancaster University,
Lancaster, LA1 4YW, U.K<EMAIL_ADDRESS>Barry S. Rowlingson
Lancaster Medical School, Lancaster University, Lancaster, LA1 4YW, UK
<EMAIL_ADDRESS>Christopher Suter  Google Research
<EMAIL_ADDRESS>Jonathan M. Read  Lancaster Medical School, Lancaster
University, Lancaster, LA1 4YW, UK<EMAIL_ADDRESS>Gareth
Roberts  Department of Statistics, University of Warwick, Coventry, CV4 7AL,
U.K<EMAIL_ADDRESS>Chris P Jewell  Department of Mathematics
and Statistics, Lancaster University, Lancaster, LA1 4YW. U.K.
<EMAIL_ADDRESS>
###### Abstract

In the event of a disease outbreak emergency, such as COVID-19, the ability to
construct detailed stochastic models of infection spread is key to determining
crucial policy-relevant metrics such as the reproduction number, true
prevalence of infection, and the contribution of population characteristics to
transmission. In particular, the interaction between space and human mobility
is key to prioritising outbreak control resources to appropriate areas of the
country. Model-based epidemiological intelligence must therefore be provided
in a timely fashion so that resources can be adapted to a changing disease
landscape quickly. The utility of these models is reliant on fast and accurate
parameter inference, with the ability to account for large amount of censored
data to ensure estimation is unbiased. Yet methods to fit detailed spatial
epidemic models to national-level population sizes currently do not exist due
to the difficulty of marginalising over the censored data. In this paper we
develop a Bayesian data-augmentation method which operates on a stochastic
spatial metapopulation SEIR state-transition model, using model-constrained
Metropolis-Hastings samplers to improve the efficiency of an MCMC algorithm.
Coupling this method with state-of-the-art GPU acceleration enabled us to
provide nightly analyses of the UK COVID-19 outbreak, with timely information
made available for disease nowcasting and forecasting purposes.

###### keywords:

Stochastic Compartmental Models, COVID-19, Bayesian statistics, MCMC

##  1 Introduction

During the COVID-19 pandemic, accurate situational awareness and the ability
to project case numbers forward in time was an important aspect of adaptive
infection control policy. In the UK, as elsewhere, the epidemic was
characterised by dramatically fluctuating case numbers over space and time,
with population behaviour interacting with pathogen transmission to create a
complex and highly-variable disease landscape. Understanding this variability
in case numbers became key to controlling the epidemic, and quantitative
modelling of the drivers of infection transmission were key to determining
quantities such as the time-varying reproduction number, intrinsic growth
rate, and the impact of intervention strategies.

Epidemics dynamics are largely driven by the characteristics of the population
at risk: the behavioural determinants of how individuals interact,
socioeconomics and spatially-varying population density and environment in
which people live (Grenfell et al., 2002). In particular, the spatial
distribution and mobility of a population has a marked impact on the pattern
of disease transmission, and understanding how epidemic behaviour responds to
individuals’ mobility is a key aspect of outbreak control (Keeling et al.,
2001). At both a local and national level, disease transmission may be assumed
to be facilitated by individual-to-individual contact, and understanding how
contact behaviour spreads infection is key to designing successful control
interventions whilst minimising public disturbance (Keeling and Eames, 2005).

Stochastic individual-level state-transition models of epidemics have been
extensively used to study the effects of spatial population structure on
epidemic spread (Keeling et al., 2001). They are a special case of the more
general class of state-transition processes, which consider individuals as
transitioning between a number of discrete, mutually-exclusive states
reflecting the expected natural history of infection according to a stochastic
process (Bartlett, 1964). For example, the classic SEIR model (Figure 1)
assumes individuals begin as _susceptible_ to infection, and transition
sequentially to _exposed_ (i.e. infected but not yet infectious), _infectious_
, and finally _removed_ (i.e. read or recovered with immunity from further
infection). Of particular interest is the infection rate, or hazard rate of
transitioning from _susceptible_ to _exposed_ , which given a suitable model
can provide valuable insights into how interactions between individuals
modulate the propensity for disease transmission as a function of spatial
separation, as well as how individual-level covariates affect disease
susceptibility and infectivity (e.g. Jewell et al., 2009; Smieszek et al.,
2011). These models are particularly popular due to their interpretability and
compatibility with decision-making questions (e.g. Tildesley et al., 2006;
Palomo-Briones et al., 2022).

Though powerful, individual-level models are computationally intensive, with
fast forward simulation and inference methods requiring model-specific
optimisations, and inference methods requiring application-specific
approximations to allow fast computation of likelihood functions (Deardon et
al., 2010; Brand et al., 2015; Sellman et al., 2018; Probert et al., 2018).
Additionally, whilst such methods have been used successfully for population
sizes up to approximately 200,000 (commonly suiting livestock disease
applications as in Probert et al. (2018)), individual-level covariate data is
rarely available to support modelling at the individual level for national and
international human populations. In such situations, the population may be
stratified by group-level covariates of interest. Stratification by space
leads to the metapopulation model, in which individuals are aggregated into a
number of discrete spatial units, with individuals within a spatial unit
sharing a set of common covariates (Levins, 1969). The exchangability of
individuals within each metapopulation hence offers performance gains in both
computatation time and storage.

For models with modest numbers of metapopulations (low-dimensional models), a
popular approach is to represent the state transition model as a system of
ordinary differential equations, using a continous mean-field approximation to
the discrete space of numbers of individuals transiting through the state-
transition graph (e.g. Lipshtat et al., 2021). However, increasing numbers of
metapopulations (high-dimensional models) divides the population into smaller
units, such that the continuous-space approximation of the ODE systems fail.
In this case, a stochastic implementation is required that captures the random
nature of infection transmission among small numbers of individuals within the
metapopulation. Whilst metapopulation models can be implemented as a
continuous time-inhomogeneous Poisson process (e.g. Minh et al., 2011), a more
computationally convenient setup is the discrete-time “chain-binomial” model
(Becker, 1981). This has the additional computational performance advantage
that complexity may be controlled by packaging transition events into time-
quanta, rather than having to compute for each individual transition event in
continuous time (Diekmann et al., ).

In a real-time outbreak response context, the utility of epidemic models is
maximised through principled parameter inference, providing important policy
information through improved situational awareness and forecasting of the
ongoing epidemic in space and time (Jewell et al., 2009; Birrell et al.,
2021). Given complete epidemic information, i.e. the epidemiological states
and transition times for each individual over the course of the epidemic, the
likelihood function of a given state-transition model is tractable and
amenable to conventional inference methods. However, in all real-world
situations inference is complicated by the inability to directly observe
individuals’ epidemiological states, or certain transition events. For
example, whereas the times of individual symptom onset or recovery might be
well-recorded, infection events are typically unobserved (Chapman et al.,
2018). The demands on inference methods for real-time epidemic model fitting
therefore present a considerable statistical challenge: not only must the
method be rapid in its own right, but must _also_ be flexible enough to
marginalise over potentially many censored transition events and
epidemiological states (Swallow et al., 2022).

In this paper we shall adopt a Bayesian approach to inference for stochastic
metapopulation models, as demonstrated by our COVID-19 application. This has
many advantages: it provides a principled and coherent calculus for the
measurement of uncertainty both in parameters and predictions of future
epidemic scenarios; it gives a natural way of taking into account missing and
censored data using data augmentation; it provides a natural framework for
experts to input knowledge to inform inference (via prior elicitation); and
finally its implementation is facilitated through powerful MCMC techniques.
That being said, MCMC methods are often seriously challenged by spatial
applications at scale due to the strong correlations inherently present within
such models. For individual-level models, data-augmentation methodology has
become state-of-the-art for unbiased parameter estimation, allowing efficient
marginalisation over the censored data (Neal and Roberts, 2004; Jewell et al.,
2009; Probert et al., 2018; Chapman et al., 2018). However, these methods
scale poorly to large population sizes, and are therefore unsuitable for the
COVID-19 at hand. Whilst particle filter methods are beginning to show promise
for fitting larger-scale epidemic models, they are inherently constrained by
poor performance as the dimensionality of the population stratification
increases (Rimella et al., 2023). Inspired by existing data-augmentation
methodology, we therefore develop a novel state-of-the art MCMC method suited
to metapopulation models, and which are amenable to rapid implementation
through efficient multi-core computing. Our innovations allow us to analyse
national-scale outbreaks in a timely manner, to provide estimates of
quantities such as the degree of interaction between population strata and
reproduction number, and predict the ongoing disease trajectory, we are
interested in capturing spatiotemporal dependence between disease prevalence
and incidence.

Between May 2020 and March 2022, the modelling approach to the UK’s COVID-19
outbreak we describe in this paper was used to provide information to the UK
government, via the Scientific Pandemic Influenza Group on Modelling,
Operational sub-group (SPI-M-O) of the Scientific Advisory Group on
Emergencies (SAGE).

pt$\mathrm{S}$$\mathrm{E}$$\mathrm{I}$$\mathrm{R}$$\lambda^{\scaleto{\mathrm{SE}}{0.75ex}}_{i}(t)$$\lambda^{\scaleto{\mathrm{EI}}{0.75ex}}_{i}(t)$$\lambda^{\scaleto{\mathrm{IR}}{0.75ex}}_{i}(t)$
Figure 1: SEIR model showing states (boxes), transitions (arrows), and
associated transition rates.

###  1.1 Spatial COVID-19 cases in the UK

Data on the daily number of new cases of COVID-19 in each of 382 Local
Authority Districts (LADs) in the UK is available from the UK Government
Coronavirus website (Public Health England, 2020, 2021). This dataset contains
the number of people $y_{it}$ who submitted a positive test sample for
COVID-19 on day $t$ in LAD $i$, where a positive test is defined as positive
by PCR or by LFD followed by confirmatory PCR. In practice we further
spatially aggregated Cornwall and the Isles of Scilly due to the latter’s
small population size, and likewise City of London and City of Westminster for
a total of 380 discrete spatial units.

From May 2020 until March 2022, we ran our analysis daily on an 84 day sliding
window. In practice, for an analysis on day $t$, we used data in the interval
$[t-88,t-4)$, discarding the latest 4 days’ worth of data due to significant
recording lag. For the purposes of exposition in this paper, however, we
retrict our results to the 84 day window between 7th June 2021 and 31th August
2021 inclusive. This covered the emergence of the SARS-CoV-2 Delta (B.1.617.2)
strain, and was after social distancing regulations had been relaxed; Figure
2.

Connectivity between LADs was measured using freely available Census 2011
commuting volume data aggregated from Middle Super Output Area (MSOA) onto our
380 LAD-based spatial units. These data provide a matrix $G$ with
$g_{ij}\;i,j=1,\dots,380$ being the number of journeys made from “Residence”
$j$ to “Workplace” $i$. $G$ is non-symmetric, reflecting commuting behaviour
rather than the reciprocity of disease transmission. We therefore calculated a
symmetric matrix $C=G+G^{T}$ of the daily number of journeys between each LAD
assuming that commuters return to their Residence each day, and go from their
Residence to their Workplace and back at most once per day. We also set
$C_{ii}=0,\;\mbox{for all}\;i$ as within-LAD infection transmissibility is
delegated to another part of our model (Section 2).

Finally, we introduce LAD-level population estimates and LAD polygon area to
inform how transmission varies with population size. In all datasets, two
pairs of LADs (Cornwall and Scilly, and City of Westminster and City of
London) are merged to allow mapping of MSOAs onto LADs resulting in data
geolocated to 380 spatial units.

Figure 2: Cases of COVID-19 in the UK, as determined by daily numbers of
Pillar 2 positive tests. Left: COVID-19 daily case counts between 7th June and
29th August 2021 inclusive, showing marked weekly variation and long-term
drift; Right: Spatial distribution of COVID-19 cases as of 29th August 2021,
expressed as incidence per 100000 people per day. Concentrations of high
incidence are observed in Scotland, Northern Ireland, and Wales.

Data availability: All data used in this analysis were public at the time of
publication. COVID-19 case data is available from the UK Government
Coronavirus website (Public Health England, 2020, 2021). Snapshots of the
covariate data used for the analysis may be found in the accompanying software
archive at https://doi.org/10.5281/zenodo.7715657.

##  2 Model

In the following sections, we construct a model to represent the epidemic in
terms of the number of individuals testing positive on each day
$t=0,\dots,T-1$ within spatial unit $i=1,\dots,M$ defined by UK Local
Authority Districts. Since our motivation is to measure quantities such as the
degree of disease transmission between spatial units and reproduction number,
and predict the ongoing disease trajectory, we are interested in capturing
spatiotemporal dependence between disease prevalence and incidence.

Below, we first describe a linear state transition model that captures the
natural history of the COVID-19 disease progression per individual, accounting
for segregation of the population at risk into discrete spatial units. We then
describe how this model is set up as discrete time Markov process,
establishing a notation allowing us to describe our inference methodology in
the next section.

###  2.1 State transition model

Within each spatial unit, we model the COVID-19 epidemic by assuming that at
any time $t$ individuals exist in one of 4 mutually exclusive disease states:
susceptible, exposed (infected but not yet infectious), infectious, and
removed which we denote by $\mathrm{S}$, $\mathrm{E}$, $\mathrm{I}$, and
$\mathrm{R}$ respectively. We assume that individuals progress sequentially
between the disease states, such that the allowed transitions are
$[\mathrm{S}\mathrm{E}]$, $[\mathrm{E}\mathrm{I}]$, and
$[\mathrm{I}\mathrm{R}]$ as shown in Figure 1. Each transition is associated
with a hazard rate, $\lambda^{\scaleto{\mathrm{SE}}{0.75ex}}_{i}(t)$,
$\lambda^{\scaleto{\mathrm{EI}}{0.75ex}}_{i}(t)$, and
$\lambda^{\scaleto{\mathrm{IR}}{0.75ex}}_{i}(t)$ respectively giving the
transition rate between states for a single individual, and allowing us
flexibility to specify how they evolve over space and time.

####  2.1.1 Infection rate

The infection hazard rate $\lambda^{\scaleto{\mathrm{SE}}{0.75ex}}_{i}(t)$ is
space- and time-dependent, parameterised as a product of a log-linear term
describing the susceptibility of an individual in $i$ to infection, and a
function of infectious challenge assuming homogeneous mixing of individuals
within each spatial unit $i$ and mixing due to commuting flows between $i$ and
all other spatial units $j\neq i$. We scale the hazard rate by the local
population size $n_{i}$ following the usual assumption of frequency-dependent
disease transmission. Thus

$\lambda^{\scaleto{\mathrm{SE}}{0.75ex}}_{i}(t)=\frac{\exp(u_{t}+s_{i})}{n_{i}}\left[x^{\scaleto{\mathrm{I}}{0.75ex}}_{it}+\psi\sum_{j\neq
i}\frac{c_{ij}x^{\scaleto{\mathrm{I}}{0.75ex}}_{jt}}{n_{j}}\right]$ (1)

where $u_{t}$ is a temporally-correlated random effect, $s_{i}$ a spatially-
correlated random effect, $X^{\scaleto{\mathrm{I}}{0.75ex}}_{it}$ is the
number of infectious individuals in $i$ at time $t$ (see below), and $c_{ij}$
represents the commuting flow to spatial unit $i$ from spatial unit $j$. The
coefficient $\psi$ is assumed unknown and is to be estimated given the data.

Due to the considerable variability in the infection dynamics of the COVID-19
outbreak we assume a random walk for the temporally-correlated random effect
$\bm{u}$, such that

$u_{t}=\begin{cases}\alpha_{0}&\mbox{if}\;t=0\\\
u_{t-1}+\alpha_{t}&\mbox{if}\;t>0\end{cases}$ (2)

with

$\displaystyle\alpha_{0}$ $\displaystyle\sim$
$\displaystyle\mbox{Normal}(0,10)$ (3) $\displaystyle\alpha_{t}$
$\displaystyle\overset{\mbox{\tiny iid}}{\sim}$
$\displaystyle\mbox{Normal}(0,\sigma_{u}^{2}),\;t>0$ (4)

and $\sigma_{u}=0.005$ chosen heuristically.

Furthermore, we assume that the spatially-correlated random effect follows a
Conditional Autoregressive (CAR) process with joint distribution

$\bm{s}\sim\mbox{MVN}\left(\bm{0},\Sigma\right).$ (5)

The covariance matrix $\Sigma$ is specified in terms of an adjacency matrix
$W$ with elements $(W)_{ij}=1\iff i\sim j$ and 0 otherwise, and a diagonal
matrix $D_{w}$ with element $(D_{w})_{ii}=W\cdot\mathbf{1}$ (i.e. the number
of spatial units neighbouring $i$) such that

$\Sigma=\sigma_{s}^{2}(D_{w}-\rho W)^{-1}$ (6)

with $\sigma_{s}$ an unknown parameter and correlation parameter $\rho=0.25$
chosen heuristically to give a smooth posterior surface across LADs in an
effort to improve reliability of model fitting by keeping the number of free
parameters to a minimum. We note that in a less operational context, a Besag-
York-Mollie model might be preferred, though would require the development of
more complex and bespoke MCMC methodology to fit in the presence of latent
event times.

The following prior distributions are imposed on the remaining free parameters

$\displaystyle\psi$ $\displaystyle\sim$ $\displaystyle\mbox{Gamma}(0,100)$ (7)
$\displaystyle\sigma_{s}$ $\displaystyle\sim$
$\displaystyle\mbox{HalfNormal}(0.1)$ (8) $\displaystyle\gamma_{1}$
$\displaystyle\sim$ $\displaystyle\mbox{Normal}(0,100)$ (9)

####  2.1.2 Latent and infectious periods

The latent period (sojourn in E) and infectious period (sojourn in I) are
parameterised through the transition rates
$\lambda^{\scaleto{\mathrm{EI}}{0.75ex}}_{i}(t)$ and
$\lambda^{\scaleto{\mathrm{IR}}{0.75ex}}_{i}(t)$. Evidence from the literature
suggests a mean latent period of 5 - 7 days (Quesada et al., 2021), with a
mean infectious period of 5 days starting 2 days before the onset of symptoms
(Tian et al., 2020). In the context of our model, we therefore assume

$\lambda^{\scaleto{\mathrm{EI}}{0.75ex}}_{i}(t)=0.25$

giving an _effective_ 4 day mean latent period irrespective of time and space,
and assume

$\lambda^{\scaleto{\mathrm{IR}}{0.75ex}}_{i}(t)=e^{\gamma_{0}+\gamma_{1}g_{t}}$

irrespective of space, where $\gamma_{0}=\log 0.25$ is the (known) mean log
hazard rate for removal ($[\mathrm{I}\mathrm{R}]$) events, $\gamma_{1}$ the
(unknown) log hazard ratio for being removed on a weekday versus the weekend
and

$g_{t}=\begin{cases}-\frac{2}{7}&\mbox{if}\;t\in\\{\mbox{mon},\dots,\mbox{fri}\\}\\\
\frac{5}{7}&\mbox{otherwise}\end{cases}$

to stabilise the inference.

###  2.2 Discrete-time Markov chain implementation

In this section, we implement the SEIR state transition model described in
Section 2.1 as a discrete-time stochastic process, more specifically a
discrete-time finite state-space Markov chain using the so-called “chain
binomial” setup (e.g. Becker, 1981).

For timepoint $t$ and spatial unit $i$, let $x^{q}_{it}$ represent the number
of individuals resident in state
$q\in\\{\mathrm{S},\mathrm{E},\mathrm{I},\mathrm{R}\\}$, and $y^{qr}_{it}$ for
transition
$[qr]\in\\{[\mathrm{S}\mathrm{E}],[\mathrm{E}\mathrm{I}],[\mathrm{I}\mathrm{R}]\\}$
represent the number of transitions occurring between each state. Given an
initial state $\bm{x}_{0}$, we evolve the epidemic by iterating

$y^{qr}_{it}\sim\mbox{Binomial}\left(x^{q}_{it},p^{qr}_{i}(t)\right),\qquad[qr]\in\left\\{[\mathrm{S}\mathrm{E}],[\mathrm{E}\mathrm{I}],[\mathrm{I}\mathrm{R}]\right\\}$
(10)

and

$\displaystyle x^{\scaleto{\mathrm{S}}{0.75ex}}_{i,t+1}$
$\displaystyle=x^{\scaleto{\mathrm{S}}{0.75ex}}_{it}-y^{\scaleto{\mathrm{SE}}{0.75ex}}_{it}$
(11) $\displaystyle x^{\scaleto{\mathrm{E}}{0.75ex}}_{i,t+1}$
$\displaystyle=x^{\scaleto{\mathrm{E}}{0.75ex}}_{it}+y^{\scaleto{\mathrm{SE}}{0.75ex}}_{it}-y^{\scaleto{\mathrm{EI}}{0.75ex}}_{it}$
(12) $\displaystyle x^{\scaleto{\mathrm{I}}{0.75ex}}_{i,t+1}$
$\displaystyle=x^{\scaleto{\mathrm{I}}{0.75ex}}_{it}+y^{\scaleto{\mathrm{EI}}{0.75ex}}_{it}-y^{\scaleto{\mathrm{IR}}{0.75ex}}_{it}$
(13) $\displaystyle x^{\scaleto{\mathrm{R}}{0.75ex}}_{i,t+1}$
$\displaystyle=x^{\scaleto{\mathrm{R}}{0.75ex}}_{it}+y^{\scaleto{\mathrm{IR}}{0.75ex}}_{it}$
(14)

where the transition probabilities
$p^{qr}_{i}(t)=1-e^{-\lambda^{qr}_{i}(t)\delta t}$ are assumed constant
throughout each timestep of size $\delta t$.

Considering (10) as the data generating model, a realisation of the Markov
chain results in the $T\times M$ state matrices $\bm{x}^{q}$ and event
matrices $\bm{y}^{qr}$ giving the state of the Markov chain at each timepoint
and transition events responsible for evolving the state.

##  3 Inference

In order to fit our model to the spatial timeseries of positive COVID-19 tests
across the UK, we make the assumption that the occurrence of a positive test
is equivalent to observing a $[\mathrm{I}\mathrm{R}]$ event. Furthermore,
since the $[\mathrm{S}\mathrm{E}]$ or $[\mathrm{E}\mathrm{I}]$ events are
censored, we adopt the conventional statistical notation
$z^{\scaleto{\mathrm{SE}}{0.75ex}}_{it}=y^{\scaleto{\mathrm{SE}}{0.75ex}}_{it}$
and
$z^{\scaleto{\mathrm{EI}}{0.75ex}}_{it}=y^{\scaleto{\mathrm{EI}}{0.75ex}}_{it}$
to denote the presence of censored data $\bm{z}$.

Given the data generating model given in Equations (10) and (14), and
conditional on unknown parameter vector
$\bm{\theta}=\\{\psi,\sigma_{s},\gamma_{1},\bm{\alpha},\bm{s}\\}$ and initial
conditions $\bm{x}_{0}$, the log likelihood of observing transitions $\bm{y}$,
censored transitions $\bm{z}$, and states $\bm{x}_{1:T}$ given a set of
initial conditions $\bm{x}_{0}$ and $\bm{\theta}$ is

$\displaystyle\ell(\bm{y},\bm{z},\bm{x}|\bm{x}_{0},\bm{\theta})=\sum_{t=0}^{T-1}\sum_{i=1}^{M}[$
$\displaystyle z^{\scaleto{\mathrm{SE}}{0.75ex}}_{it}\log
p^{\scaleto{\mathrm{SE}}{0.75ex}}_{i}(t)+(x^{\scaleto{\mathrm{I}}{0.75ex}}_{it}-z^{\scaleto{\mathrm{SE}}{0.75ex}}_{it})\log(1-p^{\scaleto{\mathrm{SE}}{0.75ex}}_{i}(t))+$
(15) $\displaystyle z^{\scaleto{\mathrm{EI}}{0.75ex}}_{it}\log
p^{\scaleto{\mathrm{EI}}{0.75ex}}_{i}(t)+(x^{\scaleto{\mathrm{E}}{0.75ex}}_{it}-z^{\scaleto{\mathrm{EI}}{0.75ex}}_{it})\log(1-p^{\scaleto{\mathrm{EI}}{0.75ex}}_{i}(t))+$
(16) $\displaystyle y^{\scaleto{\mathrm{IR}}{0.75ex}}_{it}\log
p^{\scaleto{\mathrm{IR}}{0.75ex}}_{i}(t)+(x^{\scaleto{\mathrm{I}}{0.75ex}}_{it}-y^{\scaleto{\mathrm{IR}}{0.75ex}}_{it})\log(1-p^{\scaleto{\mathrm{IR}}{0.75ex}}_{i}(t))]$
(17)

We estimate the joint posterior of the parameters and missing data
$\pi(\theta,\bm{z}|\bm{y},\bm{x}_{0})$ by using a Metropolis-within-Gibbs MCMC
scheme in which we draw alternately from
$\pi(\bm{\theta}|\bm{z},\bm{y},\bm{x}_{0})$ by adaptive Hamiltonian Monte
Carlo, and $\pi(\bm{z}|\bm{\theta},\bm{y},\bm{x}_{0})$ using a discrete-space
Metropolis Hastings step. A high-level description of our approach is shown in
Algorithm 1.

The adaptive Hamiltonian Monte Carlo algorithm for drawing from
$\pi(\bm{\theta}|\bm{z},\bm{y},\bm{x}_{0})$ is well described elsewhere (see
for example Carpenter et al., 2017), and therefore in the following sections
we focus on describing the methodology for drawing from
$\pi(\bm{z}|\bm{\theta},\bm{y},\bm{x}_{0})$ and
$\pi(\bm{x}_{0},\bm{z}|\bm{\theta},\bm{y})$.

Initialise $\bm{\theta}^{(0)},\bm{z}^{(0)}$

for  _$k$ in $1,\dots,K$_ do

Update $\bm{\theta}$ by HMC

for  _$l$ in $1,\dots,L$_ do

Update partially-censored events in $\bm{z}^{\scaleto{\mathrm{SE}}{0.75ex}}$

Update partially-censored events in $\bm{z}^{\scaleto{\mathrm{EI}}{0.75ex}}$

Update occult events in $\bm{z}^{\scaleto{\mathrm{SE}}{0.75ex}}$

Update occult events in $\bm{z}^{\scaleto{\mathrm{EI}}{0.75ex}}$

Update initial conditions $\bm{x}_{0}$ for events
$\bm{z}^{\scaleto{\mathrm{SE}}{0.75ex}}$

Update initial conditions $\bm{x}_{0}$ for events
$\bm{z}^{\scaleto{\mathrm{EI}}{0.75ex}}$

Algorithm 1 Overall Metropolis-within-Gibbs MCMC algorithm for sampling from
the joint posterior.

###  3.1 Drawing samples from $\pi(\bm{z}|\bm{\theta},\bm{y},\bm{x}_{0})$

In this section we describe the discrete-space Metropolis-Hastings algorithms
used to draw samples from the conditional posterior of the censored events and
initial conditions given state transition parameters and data,
$\pi(\bm{z},\bm{x}_{0}|\bm{\theta},\bm{y})$. The approach follows the
principles for data-augmentation in individual level models (see for example
Jewell et al., 2009), which involves two separate Metropolis-Hastings kernels
to explore the space of partially-censored and fully-censored events
respectively. Unlike in the continuous-time setting, for our discrete model we
also introduce a third algorithm which explores the initial conditions space.

In the context of our model, the challenge to successful updating of censored
event times is to create a proposal mechanism that respects the constraint
that $x^{q}_{it}\geq 0$ for all
$q\in\\{\mathrm{S},\mathrm{E},\mathrm{I},\mathrm{R}\\}$ and $t=0,\dots,T-1$
and $i=1,\dots,M$.

####  3.1.1 Updating partially-censored events

The proposed Metropolis Hastings kernel operates on a $M\times T$ censored
event matrix
$\bm{z}^{qr},\;[qr]\in\\{[\mathrm{S}\mathrm{E}],[\mathrm{E}\mathrm{I}]\\}$.
The algorithm proceeds by proposing to move a number of events $w$ from
$z^{qr}_{it}$ to $z^{qr}_{it+d}$ where we draw

$\displaystyle i,t$
$\displaystyle\sim\mbox{Discrete}\left(i,t:z^{qr}_{it}>0\right)$ (19)
$\displaystyle d$ $\displaystyle\sim\mbox{Discrete}\left(0\vee
t-d_{\mathrm{max}},\dots,-1,1,\dots,T\wedge t+d_{\mathrm{max}}\right)$ (20)
$\displaystyle w$
$\displaystyle\sim\mbox{Discrete}\left(1,\dots,B_{1}(z^{qr}_{it})\right)$ (21)

where $d_{max}>0$ is a tuning constant, and with bounding function

$B_{1}(z^{qr}_{it})=z^{qr}_{it}\wedge
w_{\max}\wedge\begin{cases}\min(x^{r}_{it},\dots,x^{r}_{i,t+d-1})&\mbox{if}\;d>0\\\
\min(x^{q}_{i,t-d},\dots,x^{q}_{t-1})&\mbox{if}\;d<0\end{cases}$ (22)

where $w_{\max}>0$ is also a tuning constant.

We then propose

$\displaystyle z^{qr\star}_{it}$ $\displaystyle=z^{qr}_{it}-w$ (23)
$\displaystyle z^{qr\star}_{t+\delta,i}$ $\displaystyle=z^{qr}_{t+\delta}+w$
(24)

and accept the proposal with probability

$\alpha(\bm{z}^{qr},\bm{z}^{qr\star})=\frac{\pi(\bm{z}^{qr\star}|\bm{\theta},\bm{y},\bm{x}_{0})}{\pi(\bm{z}^{qr}|\bm{\theta},\bm{y},\bm{x}_{0})}\cdot\frac{\mbox{nnz}\left(\bm{z}^{qr\star}\right)B_{1}(z^{qr\star}_{it};\bm{x}^{\star})}{\mbox{nnz}\left(\bm{z}^{qr}\right)B_{1}(z^{qr}_{it};\bm{x})}\wedge
1$ (25)

where $\mbox{nnz}\left(\bm{z}^{qr}\right)$ denotes the number of non-zero
elements in $\bm{z}^{qr}$.

####  3.1.2 Updating occult events

To explore the space of occult events, we again operate on $M\times T$
censored event matrices
$\bm{z}^{qr},\;[\mathrm{q}\mathrm{r}]\in\\{[\mathrm{S}\mathrm{E}],[\mathrm{E}\mathrm{I}]\\}$,
proposing to add or delete a number of events to randomly chosen elements. In
our model, occult events are overwhelmingly likely to occur close to the end
of the analysis time-window, and so we limit our choice of elements to those
in the last 21 days of the timeseries.

To add events, we choose an element $z^{qr}_{it}$ to update via

$\displaystyle i$ $\displaystyle\sim\mbox{Discrete}(1,\dots,M)$ (26)
$\displaystyle t$ $\displaystyle\sim\mbox{Discrete}(T-21,\dots,T)$ (27)

and propose to add $v$ events by

$v\sim\mbox{Discrete}(1,\dots,B_{2}(\bm{z}^{qr}_{it})),$ (28)

where bounding function

$B_{2}(z^{qr}_{it})=v_{\max}\wedge\min(x^{q}_{it},\dots,x^{q}_{iT})$ (29)

with $v_{\max}$ a tuning constant.

We then update $z^{qr\star}_{it}=z^{qr}_{it}+v$, and accept the proposal with
probability

$\alpha(\bm{z}^{qr},\bm{z}^{qr\star})=\frac{\pi(\bm{z}^{qr\star}|\bm{\theta},\bm{y},\bm{x}_{0})}{\pi(\bm{z}^{qr}|\bm{\theta},\bm{y},\bm{x}_{0})}\cdot\frac{21MB_{2}(\bm{z}^{qr})}{\mbox{nnz}\left(\bm{z}^{qr\star}\right)B_{3}(\bm{z}^{qr\star})}$
(30)

with $B_{3}(\bm{z}^{qr\star})$ defined below.

To delete events, we restrict our choice of $i$ and $t$ to positive elements
of $\bm{z}^{qr}$ as in Equation (19), and propose

$\displaystyle i,t$ $\displaystyle\sim\mbox{Discrete}(i,t:z^{qr}_{it}>0,t\geq
T-21)$ (31) $\displaystyle v$
$\displaystyle\sim\mbox{Discrete}(1,\dots,B_{3}(z^{qr}_{it};\bm{x}))$ (32)

where bounding function

$B_{3}(z^{qr}_{it})=z^{qr}_{it}\wedge
v_{\max}\wedge\min(x^{r}_{it},\dots,x^{r}_{iT}).$ (33)

We then accept the proposal $z^{qr\star}_{it}=z^{qr}_{it}-v$ with probability

$\alpha(z^{qr},z^{qr\star})=\frac{\pi(\bm{z}^{qr\star}|\bm{\theta},\bm{y},\bm{x}_{0})}{\pi(\bm{z}^{qr}|\bm{\theta},\bm{y},\bm{x}_{0})}\cdot\frac{\mbox{nnz}\left(\bm{z}^{qr}\right)B_{3}(\bm{z}^{qr})}{21MB_{2}(\bm{z}^{qr\star})}.$
(34)

####  3.1.3 Updating initial conditions

In the third data-augmentation kernel, we update the initial conditions matrix
$\bm{x}_{0}$. Given that we assume a closed population, such that
$\bm{x}^{\scaleto{\mathrm{S}}{0.75ex}}+\bm{x}^{\scaleto{\mathrm{E}}{0.75ex}}+\bm{x}^{\scaleto{\mathrm{I}}{0.75ex}}+\bm{x}^{\scaleto{\mathrm{R}}{0.75ex}}=\bm{N}$
the size vector of the population, we have two options for updating
$\bm{x}_{0}$. The first would be to swap individuals between pairs of adjacent
epidemiological states. However, since $\bm{x}_{0}$ is _a posteriori_
conditional on the events $\bm{z}$, the small conditional variance
$Var(\bm{x}_{0}|\bm{z},\bm{y},\bm{\theta})$ results in slow exploration of the
censored data space. Instead, we employ a _joint_ update of $\bm{x}_{0}$ and
$\bm{z}$ in which we consider a set of left-censored events occurring prior to
timepoint $t=0$ that have given rise to $\bm{x}_{0}$.

For a given transition $[\mathrm{q}\mathrm{r}]$, we begin as before by
choosing an element $z^{qr}_{it}$ to update

$\displaystyle i$ $\displaystyle\sim\mbox{Discrete}(1,\dots,M)$ $\displaystyle
t$ $\displaystyle\sim\mbox{Discrete}(0,\dots,6).$

noting that we restrict the choice of $t$ to a window $[0,7)$ at the beginning
of the epidemic – this helps to reduce the possible dimensionality of the
discrete random walk, since it is unlikely that left-censored events could be
moved to later timepoints without large changes in the posterior leading to a
rejected move.

We now choose to either move events forwards ($d=1$) or backwards ($d=-1$) in
time with equal probability such that

$d\sim\mbox{Discrete}(-1,1)$

and a number of events to move such that

$h\sim\mbox{Discrete}(1,\dots,B_{4}(\bm{x}_{0}))$ (35)

with

$B_{4}(\bm{x}_{0})=h_{\max}\wedge\begin{cases}z^{qr}_{it}\wedge\min(x^{q}_{i1},\dots,x^{q}_{it})&\mbox{if}\;d=-1\\\
\min(x^{r}_{i0},\dots,x^{r}_{it-1})&\mbox{if}\;d=1)\end{cases}$

and where $h_{\max}>1$ is a tuning constant. We then let

$\displaystyle x^{q\star}_{i0}=x^{q}_{i0}+dh$ $\displaystyle
x^{r\star}_{i0}=x^{r}_{i0}-dh$ $\displaystyle z^{qr\star}_{it}=z^{qr}_{it}+dh$

and accept the proposal with

$\alpha(\bm{x}_{0},\bm{z},\bm{x}_{0}^{\star},\bm{z}^{\star})=\frac{\pi(\bm{x}_{0}^{\star},\bm{z}^{\star}|\bm{y},\bm{\theta})}{\pi(\bm{x}_{0},\bm{z}|\bm{y},\bm{\theta})}\cdot\frac{B_{4}(\bm{x}_{0}^{\star})}{B_{4}(\bm{x}_{0})}.$

We note that if $B_{4}=0$, then $f_{k}(h)=0\;\forall\;h$ in Equation 35 and
the move will be rejected.

###  3.2 Definition of $R_{t}$

We define an approximate stratified temporal reproduction number $R_{jt}$ as
the expected number of further individuals that an individual in stratum $j$
will go on to infect given the state at time $t$. We define this by first
considering the pairwise force of infection defined by $\bm{K}$ exerted by an
individual in $j$ on a susceptible individual in $i$ such that

$\bm{R}_{t}\approx\frac{1-\exp\left(-\left(\bm{X}^{\scaleto{\mathrm{S}}{0.75ex}}_{\cdot
t}\right)^{T}\cdot\bm{K}\right)}{1-\exp(-\gamma_{0})}$ (36)

We remark that this is an approximation since both
$\bm{X}^{\scaleto{\mathrm{S}}{0.75ex}}_{\cdot t}$ (and therefore
$\lambda^{\scaleto{\mathrm{SE}}{0.75ex}}(t)$) is assumed constant over the
course of an individual’s infectious period.

###  3.3 Software implementation

The model and sampler code were implemented in Python 3.8 using TensorFlow and
TensorFlow Probability computational and probabilistic programming libraries
for GPU acceleration (Abadi et al., 2015; Dillon et al., 2017). Version 1.0 of
the Python code implementing this analysis is freely available under the MIT
license at https://gitlab.com/chicas-covid19/covid19uk, with the version
snapshot available at https://doi.org/10.5281/zenodo.7715657.

##  4 Sampler optimisation

Before addressing the real-world application, we demonstrate that the
discrete-space samplers described above are optimised at the conventional
accept/reject ratio of 0.23. Here we provide results for tuning $m$ the number
of metapopulations and $w$ the number of events for the partially-observed
event time moved described in Section 3.1.1, as applied to the
$[\mathrm{S}\mathrm{E}]$ and $[\mathrm{E}\mathrm{I}]$ transitions
respectively.

We first define the root mean squared jumping distance to be

$MSJD=\frac{1}{K}\sum_{k=1}^{K}\sum_{m,t,s}\left(x^{(k)}_{mtl}-x^{(k-1)}_{mts}\right)^{2}$

where $\bm{x}^{(k)}$ denotes the $M\times T\times L$ state tensor at iteration
$k=1,\dots,K$, for $M$ metapopulations, $T$ timepoints, and $L$
epidemiological compartments.

We then run the partially-observed event time sampler separately for the
$[\mathrm{S}\mathrm{E}]$ and $[\mathrm{E}\mathrm{I}]$ transitions
respectively, considering all other censored data and parameters fixed at
values taken from a randomly-chosen iteration of the converged MCMC chain
described in Section 5. For each transition, we run the sampler for
$(m,x_{max})\in\\{m:1,\dots,10\\}\times\\{x_{max}:1,\dots,100\\}$ with $m$ and
$x_{max}$ as defined in Section 3.

For each transition and $(m$, $x_{max})$ tuple, we plot the MSJD against the
acceptance ratio in Figure 3. This confirms that the MSJD is maximised at an
acceptance ratio of approximately 0.234 irrespective of the transition.
However, the way in which $(m$, $x_{max})$ affects the MSJD differs depending
on the transition selected. Moreover, we observe that the the total number of
events moved across all selected metapopulations is conserved, such that as
long as $mx_{max}\approx 30$ we achieve the maximum MSJD.

Figure 3: Mean squared jumping distance versus acceptance probability for the
$[\mathrm{S}\mathrm{E}]$ (left) and $[\mathrm{E}\mathrm{I}]$ (right)
optimisation studies respectively (Section 4). Blue dots represent each
combination of $m$ and $x_{max}$, and red vertical line indicates
theoretically optimal acceptance rate. Figure 4: Mean squared jumping distance
(left) and Metropolis Hastings acceptance probability (right) for the _move_
algorithm for both the $[\mathrm{S}\mathrm{E}]$ (top) and
$[\mathrm{E}\mathrm{I}]$ (bottom) transitions.

##  5 Application to the UK COVID-19 epidemic

In this section we apply the model described in Section 2 to the UK COVID-19
spatial case timeseries shown in Section 1.1. The MCMC algorithm were run for
$k=40000$ iterations, with $l=380$ censored event updates per iteration
(Algorithm 1) with tuning constants as in Table 1. Three independent MCMC
algorithms were run, initialising each Markov chain using a random draw from
the prior distribution, and discarding the first 10000 iterations for any
calculated quantities ($R_{t}$, riskmaps, and predictive distributions).

Table 1: Tuning constants used for MCMC algorithm to fit the spatial stochastic meta-population model to COVID-19 case timeseries data. Transition | Update | Tuning constant | Value
---|---|---|---
$[\mathrm{S}\mathrm{E}]$ | Partially-censored | $d_{max}$ | 84
$w_{max}$ | 22
Fully-censored | $v_{max}$ | 20
Initial conditions | $h_{max}$ | 19
$[\mathrm{E}\mathrm{I}]$ | Partially-censored | $d_{max}$ | 84
$w_{max}$ | 18
Fully-censored | $v_{max}$ | 20
Initial conditions | $h_{max}$ | 17

Traceplots of the scalar quantities $\exp(\alpha_{0})$, $\exp(\gamma_{1})$,
$\psi$ are shown in Figure 5, superimposing the three independent chains and
calculating the Brookes-Gelman-Rubin statistics ($\hat{R}$, not to be confused
with the time-varying reproduction number $R_{t}$) for each parameter (Brooks
and Gelman, 1998). The algorithm exhibits satisfactory convergence for all
three chains, after a conservative 10000 iteration burn-in which is removed to
compute the following now-casting and predictive results.

Figure 5: Traceplots of uni-dimensional parameters $\exp(\alpha_{0})$ (i.e.
$\alpha_{0}$ expressed as absolute infection risk per individual per day),
$\exp(\gamma_{1})$ (i.e. $\gamma_{1}$ expressed as relative risk), and $\psi$
(infection rate per commute visit). Figure 6: Posterior distribution of
$R_{t}$ over the analysis period.

The temporal trend in the global time-varying reproduction number is shown in
Figure 6, demonstrating marked variation through time either side of $R_{t}=1$
as expected given the case timeseries trajectory in Figure 2. The credible
intervals around $R_{t}$ are narrow, reflecting the choice of case observation
model and the large number of daily cases occurring in the UK as a whole, with
variation due to spatial connectivity and risk accounted for by the model.

Figure 7: Posterior mean spatial reproduction number $R_{it}$ (left) and
Population Attributable Fraction of the infection hazard due to between-LAD
mobility (right) as of 31st August 2021.

At the LAD-level, Figure 7 (left) shows the mean posterior local reproduction
number $R_{it}$ as of the 31st August 2021. Considerable spatial heterogeneity
is seen, driven only in part by the spatial distribution of cases on the same
date (Figure 2). This reflects the dynamics of the spatial epidemic over the
whole time-window, with areas such as East Anglia exhibiting higher $R_{it}$
values than might be concluded simply by looking at the most recent case data.

A convenient way to study the importance of spatial transmission is to plot
the attributable fraction of the total infection risk on an individual in LAD
$i$ which it due to between-LAD transmission, defined as

$AF_{it}=\frac{\phi\sum_{j\neq
i}\frac{c_{i}jx^{\scaleto{\mathrm{I}}{0.75ex}}_{jt}}{n_{j}}}{x^{\scaleto{\mathrm{I}}{0.75ex}}_{it}+\phi\sum_{j\neq
i}\frac{c_{i}jx^{\scaleto{\mathrm{I}}{0.75ex}}_{jt}}{n_{j}}}.$ (37)

For a particular LAD of interest, the importance of spatial transmission
therefore varies with both connectivity and within-LAD disease prevalence. For
the most recent timepoint in our case timeseries, we plot the posterior mean
AF per LAD. Highly connected LADs, such as those comprising or close to the
major cities, typically have a higher AF than more rural LADs, driven by
overall population migration.

Whilst our model captures the national-level baseline transmission rate and
inter-LAD connectivity, it also allows local variation in baseline
transmission rate through our spatially-correlated term $\bm{s}$. The
posterior mean of $\bm{s}$ and exceedance probability $Pr(s>0)$ is plotted
geographically in Figure 8. In general, these results indicate a trend towards
higher baseline transmission in less densely populated, rural LADs compared to
more populated, urban LADs.

Figure 8: Posterior mean of $\bm{s}$ (left), and exceedance probability
$q=Pr(s>0)$ (right) for each Local Authority District.

Posterior predictive checking of our model was performed by comparing the 90%
credible interval of the in-sample smoothing distribution of
$[\mathrm{I}\mathrm{R}]$ transitions against observed case numbers for the
last two weeks of our analysis window. Out-of-sample checking was performed by
comparing the 90% prediction interval of the 2-week-ahead predictive
distribution against the subsequent 2 weeks’ worth of case data. These
comparisons are plotted for the Lancaster LAD in Figure 9. We see that our
model captures both the modest increase and weekly periodicity of the case
numbers well, accommodating the apparent “catch-up” in cases seen on a Monday
after the propensity to test-and-report less at weekends. In practice, these
plots are useful to detecting LADs departing from the expected epidemic
trajectory early, so as to alert public health authorities to either case
surges or unexpectedly low incidence in individual LADs.

Figure 9: Smoothing distribution (left of vertical line) and predictive
distribution (right of vertical line) of number of $[\mathrm{I}\mathrm{R}]$
events. The smoothing distribution is compared against the observed case data
for the Lancaster LAD for the last 2 weeks of the analysis period, 18th June
2021 – 31st August 2021, whilst the predictive distribution is shown for a 2
week prediction horizon with subsequently observed cases superimposed.

##  6 Discussion

In this paper we have developed an MCMC-based method for fitting a discrete-
time spatial stochastic SEIR model to geolocated case timeseries data for the
COVID-19 pandemic in the UK. We address particularly the challenge of
providing local now-casts and short-term forecasts of epidemic spread, based
on knowledge of population structure and mobility. Our approach differs to
“classic” epidemic models (e.g. Birrell et al., 2021), in that we make no
attempt to fit our model to the entire epidemic time-series (in the case of
COVID-19 from the earliest cases in February 2020), but choose to analyse only
the most recent 12 weeks’ worth of data. This refects that the epidemic in the
UK was highly non-stationary: analysing a longer time-series yields little
further posterior information and considerably increases the computational
cost.

Our model represents a trade-off between complexity and the ability to fit it
to our 380-dimensional timeseries. It was implemented rapidly in response to
the pandemic, with the innovation being the constrained-space samplers
operating on the censored event data as in Section 3.1. To our knowledge, this
is the first time that such an MCMC scheme has been attempted for a discrete-
time model such as this. Nevertheless, it is clear that further research is
badly needed to address three key limitations of our approach as follows.

With individual-level epidemic models incorporating a high degree of
population heterogeneity, it is common to attempt to estimate the baseline
$[\mathrm{I}\mathrm{R}]$ transition rate $\gamma_{0}$ (e.g. Jewell et al.,
2009; McKinley et al., 2009; Chis Ster and Ferguson, 2007). However, it is
known that estimating $\gamma_{0}$ in the presence of censored event times is
problematic, requiring an exceptionally efficient sampler and re-
parameterisation of the model to orthogonalise $\gamma_{0}$ with respect to
the censored data (Neal and Roberts, 2005). With interpretability of the
parameters driving our model construction for this application, however, we
chose to fix $\gamma_{0}=0.25$ based on clinically-derived data. We conjecture
this is a weak assumption compared to other modelling approaches (e.g. Scott
et al., 2020), though is unable to detect if and when this transition rate
changes as a result of improved case detection or changes in the host-pathogen
interaction.

Our temporal and spatial random effects, $\bm{\alpha}$ and $\bm{s}$
respectively, were introduced following the observation that the disease
transmission rate appeared to fluctuate across time and space more than could
be explained by human mobility and the size of the population at risk. In
principle, a spatiotemporal random effect could be used, similar to empirical
methods commonly used for disease mapping. However, such methods increase the
dimensionality of the latent surface from $T+M$ to $TM$ and typically require
fast approximate methods to compute the posterior (Schrödle and Held, 2011;
Zammit-Mangion et al., 2012). Even if such a method were to be applied to this
model, the discrete nature of the censored event space would still require a
Metropolis-Hastings approach such as ours, though and extremely efficient
proposal would be required to ensure an adequate effective sample size.

A further pragmatic assumption adopted by our method is that the case
timeseries is a perfect observation of the number of $[\mathrm{I}\mathrm{R}]$
events per day. Though this assumption might hold incontravertably for highly
pathogenic disease in small populations (such as foot-and-mouth disease in
cattle), for large human populations in which social and demographic factors
affect the propensity of an individual to take a test, a stochastic
observation process would be preferable. Empirical testing showed that whilst
our method could be implemented if the $[\mathrm{I}\mathrm{R}]$ events were
treated as latent – but informed by for example a Binomial model – the
resulting algorithm was extremely slow to converge due to the extra censored
data and the interactions within the epidemic event space this engendered.

The corollary to these limitations is that although our approach certainly has
utility in informing disease control policy during an outbreak, the field of
epidemic modelling is in urgent need of improved methods for inference in the
presence of high-dimensional, correlated, and discrete censored data. Whilst
recent developments in particle filtering methods have offered promise for
improving inference in individual-level models (e.g. (Ju et al., 2021; Rimella
et al., 2023)), dependent stochastic metapopulation models such as ours still
present a challenge for high-dimensional importance sampling. In MCMC
methodology, recent advances in non-centering to orthogonalise state
transition events with respect to the model parameters, these are based on
individual-level continuous time models which fail to scale to national-level
human population sizes (Pooley et al., 2015, 2019).

Finally, all modelling approaches are contingent on reliable and timely data
to maximise their utility in prediction and policy advice. In this study, we
based our analysis on publicly available case incidence data which as argued
above is fraught with uncertainty due to changes in individuals’ preferences
towards testing uptake. However, unbiased prevalence sampling data was
collected in the UK during the pandemic, though was unavilable at the LAD-
level due to concerns over data privacy (Eales et al., 2022; House et al.,
2022). Nevertheless, the ability to use prevalence _as well as_ incidence data
for inference conditional on an epidemic model offers the possibility of
improving parameter identifiability through observations not only of the 1st-
order process (transitions), but also of the 0th-order process (epidemic
compartments). Additionally, our approach would have been greatly improved
with the addition of real-time human mobility data obtained indirectly through
methods such geolocation via cellular telephony. Such telephony data has been
shown to be highly effective as a source of covariate data in epidemic models
(Grantz et al., 2020). As it stands, our model is capable of explaining
spatial patterns of disease spread via established commuting routes, but
cannot separate spatially-homogeneous varition in human mobility due to social
determinants (holiday periods, media announcements, etc) from the underlying
transmissibility of the virus – the inclusion of accurate time-varying
mobility data in place of our fixed commuting data would offer the opportunity
to surmount this limitation. We therefore conclude by recommending that due
consideration be made to appropriate sharing of these data during a future
outbreak emergency, such that the full potential of spatial epidemic models be
realised as a resource for policy information and evidence.

## Competing interests

The authors declare no competing interests.

## Acknowledgments

ACH and CPJ were supported through the Wellcome Trust ‘GEM: translational
software for outbreak analysis’ UNS73114. ACH, CPJ, JR were supported by MRC
through the JUNIPER modelling consortium (grant number MR/V038613/1).

We are indebted to the Google Research team, who helped us to implement our
model, and fixed bugs in the underlying software libraries swiftly and
effectively.

The authors would like to thank The High End Computing facility at Lancaster
University for providing the facilities required for fitting the models in
this paper.

The views expressed in this paper are those of the authors and not necessarily
those of their respective funders or institutions.

## Data availability

The COVID-19 data were obtained from the UK Health Security Agency. These data
contain confidential information, with public data deposition non-permissible
for socioeconomic reasons. Requests for this data should be made to the UK
Health Security Agency (www.gov.uk/government/organisations/uk-health-
security-agency).

## References

  * Abadi et al. (2015) Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G. S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y. and Zheng, X. (2015) TensorFlow: Large-scale machine learning on heterogeneous systems.  URL: https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Bartlett (1964) Bartlett, M. S. (1964) The relevance of stochastic models for large-scale epidemiological phenomena.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 13, 2–8.  URL: https://rss.onlinelibrary.wiley.com/doi/abs/10.2307/2985217.
  * Becker (1981) Becker, N. (1981) A general chain binomial model for infectious diseases.  Biometrics, 251–258.
  * Birrell et al. (2021) Birrell, P., Blake, J., van Leeuwen, E., Gent, N. and De Angelis, D. (2021) Real-time nowcasting and forecasting of COVID-19 dynamics in England: the first wave.  Phil. Trans. R. Soc. B, 376, 20200279.
  * Brand et al. (2015) Brand, S., Tildesley, M. J. and Keeling, M. J. (2015) Rapid simulation of spatial epidemics : a spectral method.  Journal of Theoretical Biology, Volume 370, 121–134.  URL: https://wrap.warwick.ac.uk/67531/.
  * Brooks and Gelman (1998) Brooks, S. and Gelman, A. (1998) General methods for monitoring convergence of iterative simulations.  Journal of Computational and Graphical Statistics, 7, 434–455.
  * Carpenter et al. (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P. and Riddell, A. (2017) Stan: A probabilistic programming language.  Journal of Statistical Software, Articles, 76, 1–32.  URL: https://www.jstatsoft.org/v076/i01.
  * Chapman et al. (2018) Chapman, L., Jewell, C., Spencer, S., Pellis, L., Datta, S., Chowdhury, R., Bern, C., Medley, G. and Hollingsworth, T. (2018) The role of case proximity in transmission of visceral leishmaniasis in a highly endemic village in bangladesh.  PLOS Neglected Tropical Diseases, 12, 1–29.  URL: https://doi.org/10.1371/journal.pntd.0006453.
  * Chis Ster and Ferguson (2007) Chis Ster, I. and Ferguson, N. M. (2007) Transmission parameters of the 2001 foot and mouth epidemic in great britain.  PLOS ONE, 2, 1–12.  URL: https://doi.org/10.1371/journal.pone.0000502.
  * Deardon et al. (2010) Deardon, R., Brooks, S., Grenfell, B., Keeling, M., Tildesley, M., Savill, N., Shaw, D. and Woolhouse, M. (2010) Inference for individual-level models of infectious diseases in large populations.  Statistica Sinica, 20, 239–261.
  * (11) Diekmann, O., Othmer, H., Planqué, R., and Bootsma, M. () The discrete-time Kermack–McKendrick model: A versatile and computationally attractive framework for modeling epidemics.  Proc. Nat. Acad. Sci., 18, e2106332118.
  * Dillon et al. (2017) Dillon, J. V., Langmore, I., Tran, D., Brevdo, E., Vasudevan, S., Moore, D., Patton, B., Alemi, A., Hoffman, M. D. and Saurous, R. A. (2017) Tensorflow distributions.  CoRR, abs/1711.10604.  URL: http://arxiv.org/abs/1711.10604.
  * Eales et al. (2022) Eales, O., Wang, H., Haw, D., Ainslie, K. E. C., Walters, C. E., Atchison, C., Cooke, G., Barclay, W., Ward, H., Darzi, A., Ashby, D., Donnelly, C. A., Elliott, P. and Riley, S. (2022) Trends in sars-cov-2 infection prevalence during england’s roadmap out of lockdown, january to july 2021.  PLOS Computational Biology, 18, 1–16.  URL: https://doi.org/10.1371/journal.pcbi.1010724.
  * Grantz et al. (2020) Grantz, K., Meredith, H., Cummings, D., Metcalf, C., Grenfell, B., Giles, J., Mehta, S., Solomon, S., Labrique, A., Kishore, N., Buckee, C. and Wesolowski, A. (2020) The use of mobile phone data to inform analysis of covid-19 pandemic epidemiology.  Nature Communications, 11, 4961.
  * Grenfell et al. (2002) Grenfell, B. T., Bjørnstad, O. N. and Finkenstädt, B. F. (2002) Dynamics of measles epidemics: Scaling noise, determinism, and predictability with the tsir model.  Ecological Monographs, 72, 185–202.  URL: http://www.jstor.org/stable/3100024.
  * House et al. (2022) House, T., Riley, H., Pellis, L., Pouwels, K. B., Bacon, S., Eidukas, A., Jahanshahi, K., Eggo, R. M. and Walker, A. S. (2022) Inferring risks of coronavirus transmission from community household data.  Statistical Methods in Medical Research, 31, 1738–1756.  URL: https://doi.org/10.1177/09622802211055853.  PMID: 36112916.
  * Jewell et al. (2009) Jewell, C., Kypraios, T., Neal, P. and Roberts, G. (2009) Bayesian analysis for emerging infectious diseases.  Bayes. Anal., 4, 465–496.
  * Ju et al. (2021) Ju, N., Heng, J. and Jacob, P. E. (2021) Sequential Monte Carlo algorithms for agent-based models of disease transmission.
  * Keeling and Eames (2005) Keeling, M. J. and Eames, K. T. (2005) Networks and epidemic models.  J. R. Soc. Interface, 2, 295–307.
  * Keeling et al. (2001) Keeling, M. J., Woolhouse, M. E. J., Shaw, D. J., Matthews, L., Chase-Topping, M., Haydon, D. T., Cornell, S. J., Kappey, J., Wilesmith, J. and Grenfell, B. T. (2001) Dynamics of the 2001 uk foot and mouth epidemic: Stochastic dispersal in a heterogeneous landscape.  Science, 294, 813–817.  URL: https://www.science.org/doi/abs/10.1126/science.1065973.
  * Levins (1969) Levins, R. (1969) Some demographic and genetic consequences of environmental heterogeneity for biological control.  Bulletin of the Entomological Society of America, 15, 237–240.
  * Lipshtat et al. (2021) Lipshtat, A., Alimi, R. and Ben-Horin, Y. (2021) Commuting in metapopulation epidemic modeling.  Scientific Reports, 11, 15198.  URL: https://doi.org/10.1038/s41598-021-94672-w.
  * McKinley et al. (2009) McKinley, T., Cook, A. R. and Deardon, R. (2009) Inference in epidemic models without likelihoods.  The International Journal of Biostatistics, 5.  URL: https://doi.org/10.2202/1557-4679.1171.
  * Minh et al. (2011) Minh, P., Stevenson, M., Jewell, C., French, N. and Schauer, B. (2011) Spatio–temporal analyses of highly pathogenic avian influenza h5n1 outbreaks in the mekong river delta, vietnam, 2009.  Spatial and Spatio-temporal Epidemiology, 2, 49–57.  URL: https://www.sciencedirect.com/science/article/pii/S1877584510000572.
  * Neal and Roberts (2005) Neal, P. and Roberts, G. (2005) A case study in non-centering for data augmentation: Stochastic epidemics.  Statistics and Computing, 15, 315–327.
  * Neal and Roberts (2004) Neal, P. J. and Roberts, G. O. (2004) Statistical inference and model selection for the 1861 hagelloch measles epidemic.  Biostatistics, 5, 249–261.
  * Palomo-Briones et al. (2022) Palomo-Briones, G., Siller, M. and Grignard, A. (2022) An agent-based model of the dual causality between individual and collective behaviors in an epidemic.  Comput Biol Med, 141, 104995.
  * Pooley et al. (2019) Pooley, C., Bishop, S., Doeschl-Wilson, A. and Marion, G. (2019) Posterior-based proposals for speeding up markov chain monte carlo.  Royal Society Open Science, 6, 190619.
  * Pooley et al. (2015) Pooley, C., Bishop, S. and Marion, G. (2015) Using model-based proposals for fast parameter inference on discrete state space, continuous-time Markov processes.  J. R. Soc. Interface, 12, 20150225.
  * Probert et al. (2018) Probert, W. J. M., Jewell, C. P., Werkman, M., Fonnesbeck, C. J., Goto, Y., Runge, M. C., Sekiguchi, S., Shea, K., Keeling, M. J., Ferrari, M. J. and Tildesley, M. J. (2018) Real-time decision-making during emergency disease outbreaks.  PLOS Computational Biology, 14, 1–18.  URL: https://doi.org/10.1371/journal.pcbi.1006202.
  * Public Health England (2020) Public Health England (2020) GOV.UK Coronavirus (COVID-19) in the UK.  https://coronavirus.data.gov.uk.  Accessed: 2021-06-15.
  * Public Health England (2021) — (2021) UK LAD-level new cases by specimen date.  https://api.coronavirus.data.gov.uk/v2/data?areaType=ltla\&metric=newCasesBySpecimenDate\&format=csv\&release=2021-06-15.  Accessed: 2021-06-15.
  * Quesada et al. (2021) Quesada, J., López-Pineda, A., Gil-Guillén, V., Arriero-Marín, J., Gutiérrez, F. and Carratala-Munuera, C. (2021) Incubation period of COVID-19: A systematic review and meta-analysis.  Revista Clínica Española (English Edition), 221, 109–117.  URL: https://www.sciencedirect.com/science/article/pii/S2254887420301466.
  * Rimella et al. (2023) Rimella, L., Jewell, C. and P, F. (2023) Approximating Optimal SMC Proposal Distributions in Individual-Based Epidemic Models.  Statistica Sinica.
  * Schrödle and Held (2011) Schrödle, B. and Held, L. (2011) Spatio-temporal disease mapping using inla.  Environmetrics, 22, 725–734.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/env.1065.
  * Scott et al. (2020) Scott, J. A., Gandy, A., Mishra, S., Unwin, J., Flaxman, S. and Bhatt, S. (2020) epidemia: Modeling of epidemics using hierarchical bayesian models.  URL: https://imperialcollegelondon.github.io/epidemia/.  R package version 1.0.0.
  * Sellman et al. (2018) Sellman, S., Tsao, K., Tildesley, M. J., Brommesson, P., Webb, C. T., Wennergren, U., Keeling, M. J. and Lindström, T. (2018) Need for speed : an optimized gridding approach for spatially explicit disease simulations.  PLOS Computational Biology, 14.  URL: https://wrap.warwick.ac.uk/100972/.  ** From Crossref via Jisc Publications Router. ** Licence for AM version of this article starting on 06-04-2018: http://creativecommons.org/licenses/by/4.0/.
  * Smieszek et al. (2011) Smieszek, T., Balmer, M., Hattendorf, J., Axhausen, K., Zinsstag, J. and Scholz, R. (2011) Reconstructing the 2003/2004 H3N2 influenza epidemic in Switzerland with a spatially explicit, individual-based model.  BMC infectious diseases, 11, 115.
  * Swallow et al. (2022) Swallow, B., Birrell, P., Blake, J., Burgman, M., Challenor, P., Coffeng, L. E., Dawid, P., De Angelis, D., Goldstein, M., Hemming, V., Marion, G., McKinley, T. J., Overton, C. E., Panovska-Griffiths, J., Pellis, L., Probert, W., Shea, K., Villela, D. and Vernon, I. (2022) Challenges in estimation, uncertainty quantification and elicitation for pandemic modelling.  Epidemics, 38, 100547.  URL: https://www.sciencedirect.com/science/article/pii/S1755436522000093.
  * Tian et al. (2020) Tian, S., Hu, N., Lou, J., Chen, K., Kang, X., Z, X., H, C., Wang, D., Liu, N., Liu, D., Chen, G., Zhang, Y., Li, D., Li, J., Lian, H., Niu, S., Zhang, L. and Zhang, J. (2020) Characteristics of covid-19 infection in beijing.  Journal of Infection, 80, 401–406.  URL: https://www.sciencedirect.com/science/article/pii/S0163445320301018.
  * Tildesley et al. (2006) Tildesley, M., Savill, N., Shaw, D., Deardon, R., Brooks, S., Woolhouse, M., Grenfell, B. and Keeling, M. (2006) Optimal reactive vaccination strategies for a foot-and-mouth outbreak in the UK.  Nature, 440, 83–86.
  * Zammit-Mangion et al. (2012) Zammit-Mangion, A., Sanguinetti, G. and Kadirkamanathan, V. (2012) Variational estimation in spatiotemporal systems from continuous and point-process observations.  IEEE Transactions on Signal Processing, 60, 3449–3459.