# Inferring epidemics from multiple dependent data via pseudo-marginal methods

Alice Corbellalabel=e1 [<EMAIL_ADDRESS>Anne M
Presanislabel=e2 [<EMAIL_ADDRESS>Paul J
Birrelllabel=e3 [<EMAIL_ADDRESS>Daniela De Angelislabel=e4 [
<EMAIL_ADDRESS>Department of Statistics, University
of Warwick,  MRC Biostatistics Unit, University of Cambridge,  UK Health
Security Agency

###### Abstract

Health-policy planning requires evidence on the burden that epidemics place on
healthcare systems. Multiple, often dependent, datasets provide a noisy and
fragmented signal from the unobserved epidemic process including transmission
and severity dynamics.

This paper explores important challenges to the use of state-space models for
epidemic inference when multiple dependent datasets are analysed. We propose a
new semi-stochastic model that exploits deterministic approximations for
large-scale transmission dynamics while retaining stochasticity in the
occurrence and reporting of relatively rare severe events. This model is
suitable for many real-time situations including large seasonal epidemics and
pandemics. Within this context, we develop algorithms to provide exact
parameter inference and test them via simulation.

Finally, we apply our joint model and the proposed algorithm to several
surveillance data on the 2017-18 influenza epidemic in England to reconstruct
transmission dynamics and estimate the daily new influenza infections as well
as severity indicators as the case-hospitalisation risk and the hospital-
intensive care risk.

Epidemics,

State-Space Models,

Dependence,

Pseudo-Marginal Methods,

Evidence Synthesis,

Data Fusion,

###### keywords:

, and

##  1 Introduction

SSM have often been used to model the spread of a virus in a population, since
they provide a practical and simple framework to describe the generation of
epidemic time series (Dureau, Ballesteros and Bogich, 2013). state-space
models are a latent variable representation of a dynamical systems evolving
over specific domains (e.g. time, space, generations) (Schön et al., 2018).
The definition of an SSM involves the specification of a system of stochastic
or deterministic equations that define the behaviour of the latent process,
$X_{t}$, where $t$ indexes the domain (e.g. time), as well as the emission of
observable quantities throughout the domain, denoted by $Y_{t}$ (Birrell, De
Angelis and Presanis, 2018). Thanks to their flexibility, SSM representations
of observable phenomena have been widely used in many fields, from indoor
positioning problems in engineering (Solin et al., 2018), to environmental
studies (Anderson, 1996), to epidemic models (Magpantay et al., 2016).

Modelling epidemics as SSMs, the latent process generally encapsulates: (i) a
transmission process, describing the spread of infection through a population,
(ii) a severity process, modelling the likelihood of infected individuals to
become severe cases (e.g. symptomatic, in need of hospitalization, etc.), and
(iii) a detection process that links processes (i) and (ii) to available
observations. Moreover, (iv) other sources of randomness might affect the
generation of epidemic data, including seasonal fluctuations in parameters and
misclassification of cases unrelated to the pathogen of interest.

While formulating a SSM to describe epidemic data, there are several model
choices to make; a key choice is the level of stochasticity of the overall
system. As highlighted in the four points above, there are many possible
sources of randomness, each of them could be accounted for in several
different ways or could be overlooked through the use of deterministic
systems, a reasonable approximation when numbers are sufficiently large. A
model with a higher level of stochasticity attempts to make only few
deterministic approximations and, the more the noise is ignored, the lower the
level of stochasticity of the model.

Some models include complex dynamical components dynamic systems: Birrell et
al. (2011), for example, assumes that an age-structured deterministic
discrete-time system could approximate the spread of the 2009 pandemic
influenza; it accounts for delays between infection and detection of
symptomatic cases and places most of the randomness in the observation model.
A similar approach, with even more complex layers, is considered in other
works including Weidemann et al. (2014), analysing rotavirus data, and in
recent works analysing COVID-19 data, e.g. Keeling et al. (2020) and Birrell
et al. (2021). Conversely, other models adopt higher levels of stochasticity
and rely on simplifying assumptions on the dynamics. An example is the work of
Lekone and Finkenstädt (2006): here a discrete-time stochastic model is fitted
to Ebola data and substantial, sometimes unrealistic, assumptions are made,
including exponentially distributed time to events and full ascertainment of
the cases. These kinds of models have been particularly utilised for small-
size epidemics where the large-numbers deterministic approximation does not
hold (e.g. Nishiura (2011), Funk et al. (2018)).

Hybrid models that exploit the strength from both approaches could potentially
prove useful in many contexts. The transmission dynamics of large seasonal
epidemics and of pandemics can be reliably approximated by deterministic
systems, which, compared to their correspondent stochastic systems, allow
greater model complexity (mixing heterogeneity, many compartments, etc.). In
many contexts, however, the number of cases that become severe and detected
might be smaller and subject to significant stochastic variation. Stochastic
systems would be required to describe the severity process, marking a
difference from many models listed above, where most of the randomness lies in
the observational process.

While a model of this type (with all stochastic relationships except for the
deterministic transmission dynamics) has greater realism, it might be more
challenging under an inferential perspective, especially when multiple data
sources are considered. If individuals are detected at different levels of
severity in multiple datasets, an evidence synthesis framework (Presanis et
al., 2013; Presanis et al., 2014) should be employed, where correlations and
delays among these layers of severity (and their related data sources) are
accounted for.

This paper addresses the problem of inference of epidemics from multiple
dependent data by presenting a general framework for the inferential problem;
it identifies and tackles the challenges that arise in this context,
considering the specific case of a deterministic-transmission model with
multiple layers of stochasticity in the severity and detection stage. Section
2 provides a background to the methods used for modelling and inference,
specifically SSMs and pseudo-marginal methods. Section 3 considers their
application to the setting of epidemics, specifying a model that mirrors
realistic situations; presents its problems and proposes an inferential
routine. The aim of Section 4 is to compare the model presented in Section 3
to a similar model that, through a simplifying assumption does not account for
the intrinsic dependence between datasets, inducing under-estimation of
parameter uncertainty. Lastly Section 5 presents the analysis of real data on
cases of seasonal influenza at different levels of severity, appropriately
attributing stochasticity to each process.

##  2 State-space models and their inference

SSMs provide a flexible and general modelling framework for, among other
problems, the inference of data arising from an epidemic. They are briefly
outlined below, together with some useful tools to perform inference.

###  2.1 Notation

In this paper, Greek letters identify parameters (e.g. $\theta,\lambda$) and
their transformations ($\xi=f(\theta)$). Parameters are random variables, here
distinguished from the system’s r.v.s, usually dependent on the parameters,
denoted by uppercase Latin letters (e.g. $X,Y$) and their instances which are
lowercase ($x,y$). Lower indexes indicate the domain (e.g. $X_{t}$ is the r.v.
$X$ at time $t$) and upper indices denote the states of the system
($X^{\textsc{k}}$, ${}^{\textsc{j}}X^{\textsc{k}}$ ). Lastly, the notation
$1:r$ is used to define a running index (e.g. $y_{1:t}$ is the set of
observation $y_{1},y_{2},\dots,y_{t}$).

###  2.2 State-space models

A SSM is a stochastic process that makes use of a latent variable
representation to describe dynamical phenomena (Schön et al., 2018). It has
two components: a latent process, $\left\\{{X_{t}}\right\\}_{t\geq 0}$
representing the underlying states and their dynamics; and an observed process
$\left\\{{Y_{t}}\right\\}_{t\geq 1}$. Here, the domain of the process is
assumed to be discrete time, with intervals indexed by $t$, although the
methods reported are easily extended to other domains. A parameter-driven SSM
can be defined through the state and the observation equations, as in
Equations 1, 2 and 3, respectively

$\displaystyle X_{0}|\boldsymbol{\theta}$ $\displaystyle\sim$ $\displaystyle
p(x_{0}|\boldsymbol{\theta})$ (1)
$\displaystyle{X_{t}}|({X_{t-1}},\boldsymbol{\theta})$ $\displaystyle\sim$
$\displaystyle p({x_{t}}|{x_{t-1}},\boldsymbol{\theta})$ (2) $\displaystyle
Y_{t}|({X_{t}},\boldsymbol{\theta})$ $\displaystyle\sim$ $\displaystyle
p(y_{t}|{x_{t}},\boldsymbol{\theta})$ (3)

for $t=1,2,\dots,T$. These express the joint distribution of the initial state
of the system, $X_{0}$, of the latent states, $X_{t}$s, and of the
observations, $Y_{t}$s, parametrised by a vector $\boldsymbol{\theta}$
(Brockwell et al., 2016). Equations 1 to 3 define a Markovian model.
Markovianity is a common assumption for epidemic models and many other SSMs
that fall in the category of partially observed Markov processes (King, Nguyen
and Ionides, 2015) or hidden Markov models (Churchill, 2005).

A Markovian SSM defines the following full probability model:

$p({x_{0:T}},y_{1:T}|{\boldsymbol{\theta}})=\prod_{t=1}^{T}p(y_{t}|{x_{t}},{\boldsymbol{\theta}})\prod_{t=1}^{T}p({x_{t}}|{x_{t-1}},{\boldsymbol{\theta}})p({x_{0}}|{\boldsymbol{\theta}})$
(4)

where, thanks to Markovianity and conditional independence, the joint density
can be factorised into the state and observation densities.

A SSM can also be represented through a graphical model where a graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ represents the conditional
independence structure (edges $\mathcal{E}$) between r.v.s (nodes
$\mathcal{V}$). Some graphical models representing the state and observational
process of an SSM are illustrated in Figure 1; the dependence on the parameter
$\boldsymbol{\theta}$ is omitted for simplicity.

This framework could account for many different modelling assumptions,
including multidimensional state processes or deterministic relations, as
exemplified in Figure 1.

$\xi_{0}$$\xi_{1}$$\xi_{2}$…$\xi_{T}$$y_{1}$$y_{2}$…$y_{T}$ (a) SSM assuming
deterministic underlying relation: the states $\xi_{1:T}$ are evolving without
stochasticity.

$x^{\textsc{a}}_{0}$$x^{\textsc{a}}_{1}$$x^{\textsc{a}}_{2}$…$x^{\textsc{a}}_{T}$$x^{\textsc{b}}_{1}$$x^{\textsc{b}}_{2}$…$x^{\textsc{b}}_{T}$$y^{\textsc{b}}_{1}$$y^{\textsc{b}}_{2}$…$y^{\textsc{b}}_{T}$$x^{\textsc{c}}_{1}$$x^{\textsc{c}}_{2}$…$x^{\textsc{c}}_{T}$$y^{\textsc{c}}_{1}$$y^{\textsc{c}}_{2}$…$y^{\textsc{c}}_{T}$
(b) SSM for multiple data streams: data $y^{\textsc{b}}_{1:T}$ and
$y^{\textsc{c}}_{1:T}$ are related to processes $x^{\textsc{b}}_{1:T}$ and
$x^{\textsc{c}}_{1:T}$, which themselves are influenced by the Markov process
$x^{\textsc{a}}_{0:T}$.

Figure 1: Some examples of graphical models for SSMs. Grey nodes correspond to
observed variables and white nodes are latent variables. Solid arrows express
stochastic dependence among r.v.s, while dashed arrows express deterministic
links.

###  2.3 Inference in state-space models

Observations $\left\\{{Y_{t}}\right\\}_{t\geq 1}$ can be used to infer the
state process $\left\\{{X_{t}}\right\\}_{t\geq 0}$, usually conditional on
specific values of the static parameter $\boldsymbol{\theta}$, or to estimate
the static parameter $\boldsymbol{\theta}$, usually marginally on the
distribution of the state process.

####  2.3.1 State inference

State inference could have multiple perspectives, for a full review of these
see Lindsten and Schön (2013); in this paper the inference relies on the
estimation of the filtering distribution,
$p(x_{t}|y_{1:t},\boldsymbol{\theta})$, for $t=1,\dots,T$, which is the state
distribution at time $t$, conditional on the data and parameters.

The filtering distribution, except for few cases (e.g. linear Gaussian SSMs
(Kalman, 1960)), can be approximated recursively by a two-step procedure that
alternates approximation of the state distribution given the observations
(measurement update) and given the previous states (prediction update).

The sequential simulation of samples from the distributions
$X_{t}|y_{1:t},\boldsymbol{\theta}$ and $X_{t}|y_{1:t-1},\boldsymbol{\theta}$
allows the approximation of the likelihood of the data $y_{1:T}$, given a
parameter value $\boldsymbol{\theta}$ and marginally w.r.t. the state
distribution $X_{1:T}$ in Equation 5:

$p(y_{1:T}|\boldsymbol{\theta})=\prod_{t=1}^{T}p(y_{t}|y_{1:t-1},\boldsymbol{\theta})=\prod_{t=1}^{T}\int_{X_{0:t}}p(y_{t},x_{0:t}|y_{1:t-1},\boldsymbol{\theta})\text{
d}x_{0:t}$ (5)

The most common among these simulation-based procedures is the bootstrap
particle filter (BPF) (Arulampalam et al., 2002), which uses the state and
observation equations to iterative propose and weight samples to approximate
the filtering distribution. The BPF provides an unbiased estimator
$\widehat{p}(y_{1:T}|\boldsymbol{\theta})$ of the likelihood of the data given
a parameter value $\boldsymbol{\theta}$ which can be used to perform
inference.

A full derivation of Equation 5 and of the BPF is not necessary for the
understanding of this paper, we leave curious readers to consult Chapter 4 of
Corbella (2019) and Lindsten and Schön (2017).

####  2.3.2 Parameter inference

There is not a unique way in which this approximated likelihood can be used to
drive inference on the parameter $\boldsymbol{\theta}$. Within the Bayesian
framework, assumed here to enable a complex synthesis of available evidence,
many iterative algorithms have been developed to sample from the posterior
distribution of interest, $\boldsymbol{\Theta}|y_{1:t}$, when only an
approximation of the likelihood is available. Andrieu, Doucet and Holenstein
(2010) review and summarise the algorithms used more frequently for parameter
inference in SSMs. Among the algorithms listed, pseudo-marginal approaches,
previously introduced in Andrieu and Roberts (2009), provide a simple way to
integrate simulation-based approximation of the likelihood (such as the BPF)
into Monte Carlo Markov chain (MCMC) algorithms for Bayesian inference.

Pseudo-marginal algorithms are aimed at exploring only the posterior
distribution of the parameter, marginally from the distribution of the states,
and they are based on the classical Metropolis Hastings (MH) algorithm
(Metropolis et al., 1953; Hastings, 1970). Differently from the original MH
algorithm, here the unnormalised posterior distribution is approximated by the
product of the prior and a simulation-based approximation of the likelihood in
the acceptance ratio (e.g. the BPF). Two pseudo-marginal algorithms are
employed throughout this paper: grouped independence Metropolis Hastings
(GIMH) (Beaumont, 2003) and Monte Carlo within Metropolis (MCWM) (Andrieu and
Roberts, 2009).

###  2.4 State-space models for epidemics

The SSM methodology marries well with epidemic models: available data are only
a partially-observed signal of latent variables which encapsulate all the
processes involved in the data generation (e.g. transmission, severity,
background noise). More specifically, transmission models usually consists of
stochastic or deterministic systems of equations that describe the flow of
individuals between disjoint compartments according to their disease status:
susceptible, infectious, recovered and immune, etc.

Even the first epidemic models formulated (Kermack and McKendrick, 1927) could
be seen from a SSM perspective: the state system describes the deterministic
transmission dynamics via differential equations and the observational process
consists of the detection likelihood. This approach has endured over time: the
deterministic transmission dynamics enable the modelling of several complex
aspects of an epidemic (e.g. heterogeneous populations, non-exponential
transition time, see Keeling and Rohani (2011)) and, when large populations
are considered, they approximate well the correspondent stochastic system
(Diekmann, Heesterbeek and Britton, 2012). However this kind of model relies
on the unrealistic assumption that the only source of noise observed in the
data is caused by the case ascertainment/detection randomness.

Another stream of literature models explicitly the transitions among disease
status as random variables, assuming a stochastic state process. Before the
advent of simulation-based inference for SSMs, the inference for these models
heavily relied on simplifying model assumptions to derive closed forms of the
likelihood (Britton, 2010; Andersson and Britton, 2012).

Finally in recent decades, SSM methodology has been used for epidemic dynamics
more explicitly. Some notable works include Bretó et al. (2009); Dukic, Lopes
and Polson (2012); Dureau, Ballesteros and Bogich (2013); McKinley et al.
(2014); Shubin et al. (2016). Even if, as illustrated in Section 2.3, SSM
inference is possible (both with pseudo-marginal methods and other
algorithms), many simplifying assumptions are needed to allow computation
under time constraints. These simplifications often include assuming
independence between data stream, dropping transmission compartments, assuming
full ascertainment of the cases, and discretizing time in large intervals.
Moreover, some but not all parts of the state process of these models involve
large numbers, and could be be reliably approximated by its deterministic
counterpart.

Within the application of SSMs to epidemics, there seem to be a lack of hybrid
models that exploits the deterministic approximations of transmission dynamics
that characterise large epidemics and pandemics while proper accounting for
stochasticity in the latent states that cannot be approximated, e.g. the
states that involve severe cases. In the remainder of the paper we will
propose and discuss such a semi-stochastic model, with a particular focus on
its use when multiple dependent data are available.

##  3 Multiple data on epidemics

When a virus spreads in a population, multiple signals of its presence could
be available in the form of time-series counts. Each of these counts is likely
to be related to a specific level of severity (e.g. mild symptomatic cases,
patients that require hospitalization) and could be affected by specific
sources of noise and, possibly bias. Nevertheless, all data are linked to the
underlying system that describes the transmission process in the population.
Therefore, the joint analysis of multiple data allows better knowledge than
separate analyses of single datasets, but introduces several challenges.

###  3.1 Dependency between data streams at time $t$

Let $X^{\textsc{0}}_{t}$ be the process describing the number of new
infections at time $t$ of a specific epidemic in a given population. Infected
individuals could experience a series of increasingly severe events in a
pyramid perspective (Presanis et al., 2009), e.g. they might experience mild
symptoms, be hospitalised, require intensive care, die; with cases in each
severity layer being a subset of the cases in the previous layer as
represented in Figure 2(a). The process describing the number of individuals
experiencing an event j among these severity layers is denoted by
$X^{\textsc{j}}_{t}$, experiencing event k by $X^{\textsc{k}}_{t}$, etc. for
$t=1,\dots,T$, with $\textsc{j}=\textsc{0}$ representing infection.

$0$ j k l (a) Pyramid representation of the severity states.

${}^{0}x_{t}$$x^{\textsc{j}}_{t}$$x^{\textsc{k}}_{t}$$x^{\textsc{l}}_{t}$$y^{\textsc{k}}_{t}$$y^{\textsc{j}}_{t}$$y^{\textsc{l}}_{t}$${}^{\textsc{j}}\theta^{\textsc{k}}$${}^{0}\theta^{\textsc{j}}$${}^{\textsc{k}}\theta^{\textsc{l}}$$\zeta^{\textsc{j}}$$\zeta^{\textsc{k}}$$\zeta^{\textsc{l}}$$t=2,\dots,T$
(b) DAG representation of a stochastic-severity model.

Figure 2: Pyramid representation and DAG of a severity process.

A simple model would assume that the number of people in a severity category
would be a subsample (e.g. a Binomial sample) of the number of people in the
less-severe category with severity parameters
${}^{\textsc{j}}\theta^{\textsc{k}}$ encapsulating the risk of moving from
severity state j to k. The observational processes would describe the number
of detected cases in each specific severity layer $Y^{\textsc{j}}_{1:t}$, as,
again a subsample of the cases in a specific layer j with a layer-specific
detection probability $\zeta^{\textsc{j}}$ . Figure 2(b) illustrates a model
of this type comprising 3 events of increasing severity and the connection
among multiple data can be easily spotted: i.e. not only they are linked
through the underlying infection process, but also via the progression through
severity states. In fact, while distributional properties could be used to
write the likelihood of each of multiple datasets given the infection process
$X^{\textsc{0}}_{1:t}$ and some parameters, these data are not independent,
conditionally on the parameters: they share common severity processes
$X^{\textsc{j}},X^{\textsc{k}},X^{\textsc{l}}$, and therefore cannot be simply
multiplied to obtain a joint likelihood of the data.

####  3.1.1 Dependency between data streams across time

Another challenge is given by delays between events: severe events usually do
not happen simultaneously. More likely, there will be some time elapsing
between events of increasing severity; this could lead to, for example, cases
detected in severity layer j at time $t$, be present in data on layer k at
some time $t+s,s\geq 0$ with some chance. This introduces dependence not only
between observations at different levels of severity at the same time, but
also across time.

To include delays, a new notation is introduced: denote by
${}^{\textsc{j}}_{t}X^{\textsc{k}}$ the number of people that will eventually
experience severe event k, having already experienced event j at $t$; and by
${}^{\textsc{j}}_{t}X^{\textsc{k}}_{s}$ the number of people that experience
event k at time $s$, having already experienced event j at $t$. Thus, the
subscript and superscript on the right side denote the time and type of final
events and the ones on the left side denote the time and type of a previous
event. When only the right superscript/subscript is reported, e.g.
$X^{\textsc{k}}_{t}$, as in the previous section, this denotes the number of
people experiencing event k at $t$ irrespectively from the time of previous
event, which can be obtained summing over the number of people experiencing
event $j$ at the previous event times: $X^{\textsc{k}}_{s}=\sum_{t\leq
s}{}^{\textsc{j}}_{t}X^{\textsc{k}}_{s}$.

A stochastic model for the time series ${}^{\textsc{j}}_{1:T}X^{\textsc{k}}$
of the number of people that have had event j at time $t$ ($t=1,\dots,T$) and
will experience event k can be specified by, for example, a Binomial sample:

$\left({}^{\textsc{j}}_{t}X^{\textsc{k}}|x^{\textsc{j}}_{t}\right)\sim\text{Bin}(x^{\textsc{j}}_{t},{}^{\textsc{j}}\theta^{\textsc{k}})$
(6)

for $t=1,\dots,T$. Denote by
${}^{\textsc{j}}f^{\textsc{k}}_{d}({}^{\textsc{j}}\vartheta^{\textsc{k}})$ the
probability that the delay experienced between event j and k is in the $d$th
interval of length $\delta$, $[\delta d;\delta d+\delta)$ for $d=0,1,\dots,D$,
with $D$ being the largest interval index for which the delay is relevant
(i.e.
${}^{\textsc{j}}f^{\textsc{k}}_{d}({}^{\textsc{j}}\vartheta^{\textsc{k}})\approx
0$ for $d>D$).
${}^{\textsc{j}}f^{\textsc{k}}_{d}({}^{\textsc{j}}\vartheta^{\textsc{k}})$ is
often derived from the discretization of a parametric distribution with
appropriate parameter vector
${}^{\textsc{j}}\boldsymbol{\vartheta}^{\textsc{k}}$. To make the notation
lighter, ${}^{\textsc{j}}\vartheta^{\textsc{k}}$ is dropped, with
${}^{\textsc{j}}f^{\textsc{k}}_{d}$ representing both the function and the
parameters used to describe the delay from j to k.

Under this discrete definition of the distribution of the time to event, and
conditionally on ${}^{\textsc{j}}_{t}X^{\textsc{k}}$, the introduction of
stochastic delays could be allowed by defining
${}^{\textsc{j}}_{t}X^{\textsc{k}}_{s}$, as a component of the Multinomial
r.v.:

${}^{\textsc{j}}_{t}X^{\textsc{k}}_{t:t+D}|{}^{\textsc{j}}_{t}x^{\textsc{k}}\sim\text{Multi}\left({}^{\textsc{j}}_{t}x^{\textsc{k}},{}^{\textsc{j}}f^{\textsc{k}}_{0:D}\right).$
(7)

with
${}^{\textsc{j}}f^{\textsc{k}}_{0:D}=({}^{\textsc{j}}f^{\textsc{k}}_{0},{}^{\textsc{j}}f^{\textsc{k}}_{1},\dots,{}^{\textsc{j}}f^{\textsc{k}}_{D})$
the vector containing the probabilities for the waiting time between 0 and j
as described above. The number of people that experience event k at each time
$t=1,\dots T$ can then be obtained by summing these stochastic terms, i.e.:

$X^{\textsc{k}}_{t}=\sum_{d=0}^{D}{}_{t-d}{{}^{\textsc{j}}X^{\textsc{k}}_{t}},$
(8)

which can be recognised as a typical stochastic convolution to describe delays
in epidemic models (Brookmeyer et al., 1994). The counts of events at a higher
level of severity can be modelled likewise.

This form of delay structure, as many other stochastic delay formulations,
introduces a dependence over time. This is evident in Equation 7, where the
number of people experiencing k at time $t$ depends on r.v.s defined on the
previous $D$ intervals.

###  3.2 A model for multiple data on severe influenza

Influenza is monitored by many surveillance schemes in the UK, one of these is
the Severe Acute Respiratory Infection (SARI)-Watch scheme (UK Health Security
Agency (UKHSA), 2021), which has evolved from the UK Severe Influenza
Surveillance System (USISS), the main source of data on severe influenza cases
in the UK prior to the COVID-19 pandemic. According to the USISS protocol
(Health Protection Agency, 2011a), all National Health Service (NHS) trusts in
England report, among other things, the weekly number of laboratory-confirmed
influenza cases admitted to Intensive Care (IC) units. In addition to this
mandatory scheme, a sentinel subgroup of NHS trusts in England is recruited
every year to participate in the USISS sentinel scheme (Health Protection
Agency, 2011b; Boddington, Verlander and Pebody, 2017), which reports weekly
numbers of laboratory-confirmed influenza cases hospitalised at all levels of
care. Some individuals might be detected in both datasets, leading to a
dependence.

####  3.2.1 Parametrization and data generation

Denote by $\boldsymbol{\theta}$ the set of parameters, composed of
$\boldsymbol{\theta}^{T}$, the parameters of a $SEIR$ transmission model,
tracking the number of susceptible ($S$), exposed ($E$), infected ($I$) and
removed ($R$), individuals, and $\boldsymbol{\theta}^{S}$, the parameters of
the severity and detection model.

$\boldsymbol{\theta}^{T}=\left\\{\pi,\iota,\sigma,\gamma,\beta\right\\}$
consists of the transmission rate $\beta$; the exit rates from compartments
$E$ and $I$, $\sigma$ and $\gamma$ respectively; and the initial proportions
immune, $\pi$, and of infected/infectious, $\iota$. The parameters $\pi$ and
$\iota$, together with $\sigma$, $\gamma$ and the known constant $N$, the
total size of the population, contribute to the formulation of the initial
state of the epidemic. Seasonal influenza is a large epidemic that takes place
every winter with high likelihood, hence a deterministic transmission model is
assumed: the information contained in $\boldsymbol{\theta}^{T}$, together with
the known constants, provides the full time series of the number of new
infections. Let the time be discretised in intervals of length $\delta$ so
that the $t$-th interval covers the time $[t\delta,t\delta+\delta)$ and the
intervals are indexed by $t=0,1,2,\dots,T$, where $t=0$ coincides with the
beginning of the data collection period and $t=T$ is the end of the data
collection period. Denote the number of susceptible individuals at the
beginning of interval $t$ by $S_{t}$ and likewise for the other compartments
$E,I,R$. Denote by $\xi^{0}_{1:T}$ the vector of the number of new infections
in interval $t=0,1,2,\dots T$.

$\boldsymbol{\theta}^{S}=\left\\{{}^{0}\theta^{\textsc{h}},{}^{\textsc{h}}\theta^{\textsc{ic}},{}^{0}f^{\textsc{h}},{}^{\textsc{h}}f^{\textsc{ic}},\zeta_{t}^{\textsc{h}},\zeta_{t}^{\textsc{ic}}\right\\}$
includes severity and detection parameters: ${}^{0}\theta^{\textsc{h}}$ is the
probability of hospitalization given infection;
${}^{\textsc{h}}\theta^{\textsc{ic}}$ is the probability of IC admission given
hospitalization; ${}^{0}f^{\textsc{h}}$ and ${}^{\textsc{h}}f^{\textsc{ic}}$
denote probabilities for the times from infection to hospitalisation and from
hospital admission to IC admission, respectively; $\zeta_{t}^{\textsc{h}}$ and
$\zeta_{t}^{\textsc{ic}}$ are the probability of detecting an hospitalised and
IC case, respectively.

The transmission and severity compartments are reported in Figure 3.

$S$$E$$I$$R$$\xi^{0}$$\sigma$$\gamma$HospitalH sentinelH obsH non sentinelIC
sentinelIC non sentinelIC totalIC obs
${}^{0}\theta^{\textsc{h}}$${}^{\textsc{h}}\theta^{\textsc{ic}}$${}^{\textsc{h}}\theta^{\textsc{ic}}$$\zeta_{t}^{\textsc{h}}$$\zeta_{t}^{\textsc{ic}}$
Figure 3: Flowchart of the possible severe events recorded in USISS.

####  3.2.2 SSM formulation

Describing the model using SSM notation, the state process is composed of the
distributions of: ${}_{t}^{0}X^{\textsc{h}}$ , the number of hospitalizations
that were infected at each interval $t$; $X^{\textsc{h}}_{t}$, the number of
hospitalizations at $t$, obtained via a convolution of
${}_{t}^{0}X^{\textsc{h}}$; the number of hospitalizations that eventually
will be admitted to IC and have been hospitalised at $t$,
${}_{t}^{\textsc{h}}X^{\textsc{ic}}$; the number of IC admissions at $t$,
$X^{\textsc{ic}}_{t}$ obtained by the convolution of
${}_{t}^{\textsc{h}}X^{\textsc{ic}}$, for $t=0,1,2,\dots,T$. The state process
is described as:

$\begin{split}\left({}_{t}^{0}X^{\textsc{h}}\right)&\sim\text{Pois}({}^{0}\xi\cdot{}^{0}\theta^{\textsc{h}})\\\
\left({}_{t}^{0}X^{\textsc{h}}_{t:t+D}|_{t}^{0}X^{\textsc{h}}={{}_{t}^{0}x^{\textsc{h}}}\right)&\sim\text{Multi}(_{t}^{0}x^{\textsc{h}},{}^{0}f^{\textsc{h}}_{0:D}),\qquad
X^{\textsc{h}}_{t}=\sum_{s=0}^{S}{{}_{t-s}{{}^{0}X}^{\textsc{h}}_{t}}\\\
\left({}_{t}^{\textsc{h}}X^{\textsc{ic}}|X^{\textsc{h}}_{t}={x^{\textsc{h}}_{t}}\right)&\sim\text{Bin}(x^{\textsc{h}}_{t},{}^{\textsc{h}}\theta^{\textsc{ic}})\\\
\left({}_{t}^{\textsc{h}}X^{\textsc{ic}}_{t:t+D}|_{t}^{\textsc{h}}X^{\textsc{ic}}={{}_{t}^{\textsc{h}}x^{\textsc{ic}}}\right)&\sim\text{Multi}(_{t}^{\textsc{h}}x^{\textsc{ic}},{}^{\textsc{h}}f_{0:D}^{\textsc{ic}}),\qquad
X^{\textsc{ic}}_{t}=\sum_{s=0}^{S}{{}_{t-s}{{}^{\textsc{h}}X}^{\textsc{ic}}_{t}}\\\
\end{split}$ (9)

for $t=0,1,\dots,T$. Here ${}^{0}f^{\textsc{h}}_{0:D}$ and
${}^{\textsc{h}}f^{\textsc{ic}}_{0:D}$ are vectors containing elements
${}^{0}f^{\textsc{h}}_{d}$ and ${}^{\textsc{h}}f^{\textsc{ic}}_{d}$ denoting
the probability of experiencing a delay of $d$ weeks between infection and
hospitalization and hospitalization and IC admission, respectively, for
$d=0,\dots,D$. These are considered known and fixed.

The observational process, consists of the distributions of two datasets,
$y_{1:T}^{\textsc{h}}$, the count of hospitalizations, and
$y_{1:T}^{\textsc{ic}}$, the count of IC admissions, conditional on the hidden
states $X_{1:T}^{\textsc{h}}$ and $X_{1:T}^{\textsc{ic}}$. These are assumed
to be Binomial with detection probabilities $\zeta^{\textsc{h}}_{t}$ and
$\zeta^{\textsc{ic}}_{t}$ respectively:

$\begin{split}\left(Y_{t}^{\textsc{h}}|X^{\textsc{h}}_{t}=x^{\textsc{h}}_{t}\right)&\sim\text{Bin}(x^{\textsc{h}}_{t},\zeta^{\textsc{h}}_{t})\\\
\left(Y_{t}^{\textsc{ic}}|X^{\textsc{ic}}_{t}=x^{\textsc{ic}}_{t}\right)&\sim\text{Bin}(x^{\textsc{ic}}_{t},\zeta^{\textsc{ic}}_{t})\end{split}$
(10)

for $t=0,1,2,\dots,T$.

###  3.3 Inference

Two inferential methods are proposed here: the first does not account for the
dependence among data while the second properly encapsulates the layers of
stochasticity present in the model.

####  3.3.1 The approximation under an independence assumption

Thanks to the Poisson properties (Kingman, 1992), several hidden states and
the data distribute marginally according to a Poisson distribution:

$\begin{split}X^{\textsc{h}}_{t}&\sim\text{Pois}\left({}^{0}\theta^{\textsc{h}}\cdot\sum_{d=0}^{D}\xi^{0}_{t-d}\cdot{}^{0}f^{\textsc{h}}_{d}\right)\\\
X^{\textsc{ic}}_{t}&\sim\text{Pois}\left({}^{\textsc{h}}\theta^{\textsc{ic}}\cdot{}^{0}\theta^{\textsc{h}}\cdot\sum_{d=0}^{D}\sum_{g=0}^{d}\xi^{0}_{t-d}\cdot{}^{0}f^{\textsc{h}}_{d}\cdot{}^{\textsc{h}}f^{\textsc{ic}}_{g}\right)\\\
Y_{t}^{\textsc{h}}&\sim\text{Pois}\left(\zeta^{\textsc{h}}_{t}\cdot{}^{0}\theta^{\textsc{h}}\cdot\sum_{d=0}^{D}\xi^{0}_{t-d}\cdot{}^{0}f^{\textsc{h}}_{d}\right)\\\
Y_{t}^{\textsc{ic}}&\sim\text{Pois}\left(\zeta^{\textsc{ic}}_{t}\cdot{}^{\textsc{h}}\theta^{\textsc{ic}}\cdot{}^{0}\theta^{\textsc{h}}\cdot\sum_{d=0}^{D}\sum_{g=0}^{d}\xi^{0}_{t-d-g}\cdot{}^{0}f^{\textsc{h}}_{d}\cdot{}^{\textsc{h}}f^{\textsc{ic}}_{g}\right)\end{split}$
(11)

for $t=0,1,2,\dots T$.

Ignoring the dependence between the two data streams, the joint likelihood is

$p(y_{1:T}^{\textsc{h}}|\xi^{0}_{1:T},{}^{0}\theta^{\textsc{h}},\zeta^{\textsc{h}}_{1:T})\times
p(y_{1:T}^{\textsc{ic}}|\xi^{0}_{1:T},{}^{\textsc{h}}\theta^{\textsc{ic}},{}^{0}\theta^{\textsc{h}},\zeta^{\textsc{ic}}_{1:T}),$
(12)

where each factor is the Poisson density of Equation 11.

The independence assumption would substantially simplify inference since, with
a likelihood available in closed form, there would be no need of a simulation-
based estimation of the likelihood and, indeed, of pseudo-marginal methods.
The aim of Section 4 is to assess if this approximation induces any error in
terms of bias or uncertainty quantification.

####  3.3.2 Exact inference via pseudo-marginal methods

To properly account for dependence and stochasticity, a simulation algorithm
is proposed to approximate the joint likelihood of the hospitalization and IC
data. The joint probability distribution can be decomposed in two ways:

$\begin{split}p(y_{1:T}^{\textsc{h}},y_{1:T}^{\textsc{ic}}|\boldsymbol{\theta})&=p(y_{1:T}^{\textsc{h}}|y_{1:T}^{\textsc{ic}},\boldsymbol{\theta})p(y_{1:T}^{\textsc{ic}}|\boldsymbol{\theta})\\\
&=p(y_{1:T}^{\textsc{ic}}|y_{1:T}^{\textsc{h}},\boldsymbol{\theta})p(y_{1:T}^{\textsc{h}}|\boldsymbol{\theta})\end{split}$

where, in both cases, the second of the two factors is available in closed
form (Equation 11).

To approach the estimation of the other factor, state-inference methods for
SSMs can be used. The methods for inference described in Section 2 address the
Markovian dependence across time. Here however, the time-dependence of the
transmission process disappears thanks to the deterministic approximation.
However, the dependence over the severity domain remains: the distributional
assumptions of Equation 9 can be used to construct a simulation-based
estimator of the likelihood that sequentially approximates severity states.

Algorithm 1 exploits the first decomposition, where
$p(y_{1:T}^{\textsc{ic}}|\boldsymbol{\theta})$ is available in closed form and
a solution is needed for
$p(y_{1:T}^{\textsc{h}}|y_{1:T}^{\textsc{ic}},\boldsymbol{\theta})$, which is
obtained by approximating the $T$-dimensional integral:

$\begin{split}p(y_{1:T}^{\textsc{h}}|y_{1:T}^{\textsc{ic}},\boldsymbol{\theta})=&\int_{X^{\textsc{h}}_{1}}\dots\int_{X^{\textsc{h}}_{T}}p(y_{1:T}^{\textsc{h}},X^{\textsc{h}}_{1:T}|y_{1:T}^{\textsc{ic}},\boldsymbol{\theta})\text{d}X^{\textsc{h}}_{1}\dots\text{d}X^{\textsc{h}}_{T}\\\
=&\int_{X^{\textsc{h}}_{1}}\dots\int_{X^{\textsc{h}}_{T}}p(y_{1:T}^{\textsc{h}}|X^{\textsc{h}}_{1:T},y_{1:T}^{\textsc{ic}},\boldsymbol{\theta})p(X^{\textsc{h}}_{1:T}|y_{1:T}^{\textsc{ic}},\boldsymbol{\theta})\text{d}X^{\textsc{h}}_{1}\dots\text{d}X^{\textsc{h}}_{T}\\\
=&\int_{X^{\textsc{h}}_{1}}\dots\int_{X^{\textsc{h}}_{T}}p(y_{1:T}^{\textsc{h}}|X^{\textsc{h}}_{1:T},\boldsymbol{\theta})p(X^{\textsc{h}}_{1:T}|y_{1:T}^{\textsc{ic}},\boldsymbol{\theta})\text{d}X^{\textsc{h}}_{1}\dots\text{d}X^{\textsc{h}}_{T}\end{split}$

The simulation of the hidden states is made simple by the distribution chosen
for the severity and detection process.

Result:
$\widehat{p}(y_{1:T}^{\textsc{h}},y_{1:T}^{\textsc{ic}}|\boldsymbol{\theta})$

Input: fixed parameter $\boldsymbol{\theta}$, number of particels ${N}$, data
$y_{1:T}^{\textsc{h}}$, $y_{1:T}^{\textsc{ic}}$

compute

$p(y_{1:T}^{\textsc{ic}}|\boldsymbol{\theta})$ =
$f(y_{1:T}^{\textsc{ic}}|\zeta^{\textsc{ic}}_{t}\cdot{}^{\textsc{h}}\theta^{\textsc{ic}}\cdot{}^{0}\theta^{\textsc{h}}\cdot\sum_{d\text{=}0}^{D}\sum_{g\text{=}0}^{d}\xi^{0}_{t\text{-}d\text{-}g}\cdot{}^{0}f^{\textsc{h}}_{d}\cdot{}^{\textsc{h}}f^{\textsc{ic}}_{g})$
with $f(\cdot)$ being a Poisson density

for  _$n=1,\dots,{N}$_ do

for  _$t=0,1,\dots T$_ do

sample : ${x^{\textsc{ic}}_{t}}^{(n)}\sim$ Pois
$\left(\left(1-\zeta^{\textsc{ic}}_{t}\right)[{}^{\textsc{h}}\theta^{\textsc{ic}}\cdot{}^{0}\theta^{\textsc{h}}\sum_{d=0}^{D}\sum_{g=0}^{d}\xi^{0}_{t-d-g}\cdot{}^{0}f^{\textsc{h}}_{d}\cdot{}^{\textsc{h}}f^{\textsc{ic}}_{g}]\right)+y_{T}^{\textsc{ic}}$

sample : ${{}_{t-1}^{\textsc{h}}x^{\textsc{ic}}_{1}}^{(n)}\text{,
}\dots\text{,
}{{}_{t-S}^{\textsc{h}}x^{\textsc{ic}}_{S}}^{(n)}\sim\text{Multi}\left({x^{\textsc{ic}}_{t}}^{(n)},{}^{\textsc{h}}f^{\textsc{ic}}_{1:S}\right)$

compute :
${{}_{t}^{\textsc{h}}x^{\textsc{ic}}}^{(n)}=\sum_{s=1}^{S}{{}_{t}^{\textsc{h}}x^{\textsc{ic}}_{s}}^{(n)}$

sample :
${x^{\textsc{h}}_{t}}^{(n)}|{{}_{t}^{\textsc{h}}x^{\textsc{ic}}}^{(n)},\boldsymbol{\theta}\sim\text{Pois}\left(\left(1-{}^{\textsc{h}}\theta^{\textsc{ic}}\right)\left[{}^{0}\theta^{\textsc{h}}\sum_{d=0}^{D}\xi^{0}_{t-d}\cdot{}^{0}f^{\textsc{h}}_{d}\right]\right)+{{}_{t}^{\textsc{h}}x^{\textsc{ic}}}^{(n)}$

end for

compute :
$p\left(y_{1:T}^{\textsc{h}}|{x^{\textsc{h}}_{1:T}}^{(n)},\boldsymbol{\theta}\right)$
=
$g\left(y_{1:T}^{\textsc{h}}|{x^{\textsc{h}}_{1:T}}^{(n)},\zeta^{\textsc{h}}_{t}\right)$
with $g(\cdot)$ being a Binomial density

end for

$\widehat{p}\left(y_{1:T}^{\textsc{h}},y_{1:T}^{\textsc{ic}}|\boldsymbol{\theta}\right)=p\left(y_{1:T}^{\textsc{ic}}|\boldsymbol{\theta}\right)\cdot\frac{1}{N}\sum_{n=1}^{{N}}p\left(y_{1:T}^{\textsc{h}}|{x^{\textsc{h}}_{1:T}}^{(n)},\boldsymbol{\theta}\right)$

Algorithm 1 Approximation of the likelihood
${p}(y_{1:T}^{\textsc{h}},y_{1:T}^{\textsc{ic}}|\boldsymbol{\theta})$

A second algorithm uses the alternative factorization of the joint likelihood.
While still feasible, this algorithm performs more poorly than Algorithm 1,
mainly because the IC admissions are more-completely observed and smaller in
numbers, hence particle degradation is more likely to take place (Brooks et
al., 2011). See the Supplementary information for the derivation of the two
algorithms.

##  4 Relevance of the dependence

A full simulation study is set up to assess whether (and in which situations)
accounting for the dependence makes any difference compared to assuming the
two datasets to be independent. Intuitively, a miss-specified model that does
not account for dependencies would assume more independent information than is
truly present in the dependent data, hence would achieve overly-confident
results. Conversely, truly accounting for dependence should reflect more
properly the unknowns and noise of the model and also help the estimation of
those parameters that effectively link the dependent data streams.

###  4.1 Simulation study set-up

The simulated data are formulated to reflect a situation similar to the
motivating USISS data on a smaller population, chosen to reduce the
computation time. The datasets are generated with some common parameters and
some scenario-specific parameters.

The common parameters are:

$\left\\{N=10000,\beta=0.63,\pi=0.3,\iota=0.0001,\sigma=\frac{1}{4},\gamma=\frac{1}{3.5},{}^{0}f^{\textsc{h}}\sim\text{Exp}(0.3),{}^{\textsc{h}}f^{\textsc{ic}}\sim\text{Exp}(0.4)\right\\}$

while the scenario specific parameters are reported in Table 1. Each of the
severity and detection parameters can take either a small or a large value.
The smaller leads to situations where the probability of being observed in
both datasets is low, therefore data are less dependent; while when parameters
take larger values, there is more overlap between datasets and therefore more
dependence.

Table 1: Parameters used to generate the datasets. | small dependence | large dependence
---|---|---
${}^{0}\theta^{\textsc{h}}$ | 0.1 | 0.5
${}^{\textsc{h}}\theta^{\textsc{ic}}$ | 0.1 | 0.9
$\zeta_{t}^{\textsc{h}}=\zeta^{\textsc{h}}$ | 0.1 | 0.3
$\zeta_{t}^{\textsc{ic}}=\zeta^{\textsc{ic}}$ | 0.1 | 0.9

The values have been chosen to be at the extremes of a realistic range: e.g.,
the probability of hospitalization and the detection of hospitalization in the
large dependence case are smaller than the respective quantities for IC unit
admission because usually more severe cases are better monitored, and hence
have a higher detection.

The aim of the comparison is to assess whether a misspecified independent
likelihood (Equation 10) would affect the inference of the parameters, leading
to different posterior distributions from the ones obtained with the Monte
Carlo (MC) approximation of the joint likelihood assuming dependent data
(Algorithm 1). For this reason the results presented here are to be compared
within each scenario: between the ones obtained with the independent
misspecified model (abbreviated with miss ind) and the ones obtained using the
dependent joint model (abbreviated with joint dep).

###  4.2 Results of the simulation study

As outlined in detail below, the misspecified model leads to overly precise
results in the estimation of the transmission parameters when the dependence
is large. Furthermore, the parameter connecting the two severity states to
which the data refer is estimated with less precision when the misspecified
model is assumed.

####  4.2.1 Comparison for transmission parameters

For the parameter inference with the approximated joint likelihood, a MCWM
algorithm with likelihood approximation via Algorithm 1 with $N=2000$ was
chosen. Resulting estimates were compared with those from the misspecified
independent Poisson likelihood over 1000 synthetic datasets. 500 datasets were
simulated using the smaller values of the parameters (left column of Table 1)
and 500 datasets using the larger values (right column of Table 1). The only
parameters inferred are the transmission parameters $\beta,\iota$ and $\pi$,
with the severity parameters being fixed at their true, scenario-specific,
value.

In the case of small dependence, the results show that the posterior
distributions obtained with the misspecified independent likelihood are very
similar to the ones obtained with the approximated joint likelihood. To
measure any discrepancy between the two estimation methods, the pairwise
difference (PWD) in variance Var$(\hat{\alpha}^{m}\mid y^{d})$ for
$m\in\\{\textsc{joint dep},\textsc{miss ind}\\}$ and $\alpha$ being one of the
parameters, and in the width R${}_{95}(\hat{\alpha}^{m}\mid y^{d})$ of the 95%
Credible Intervals (CrIs) was computed as:

$\begin{split}\text{PWD}(\text{Var}(\alpha))^{d}&=\text{Var}(\widehat{\alpha}^{\textsc{joint
dep}}|\boldsymbol{y}^{d})-\text{Var}(\widehat{\alpha}^{\textsc{miss
ind}}|\boldsymbol{y}^{d})\qquad\alpha=\beta,\pi,\dots;d=1,2,\dots,500\\\
\text{PWD}(\text{R}_{95}(\alpha))^{d}&=\text{R}_{95}(\widehat{\alpha}^{\textsc{joint
dep}}|\boldsymbol{y}^{d})-\text{R}_{95}(\widehat{\alpha}^{\textsc{miss
ind}}|\boldsymbol{y}^{d})\qquad\alpha=\beta,\pi,\dots;d=1,2,\dots,500.\end{split}$

These quantities are reported in Figure 4 (top) which shows an imperceptible
difference in the posterior distributions and their precision-summaries
between the two models.

The same analysis is run on the 500 datasets with a large dependence, with
results reported in Figure 4 (bottom).

Figure 4: Histogram of the pairwise differences in variance and in 95% CrI
length of the posterior distribution of the transmission parameters $\beta$
for the small dependence scenario (left) and big dependence scenario (right).
The bottom plots report the posterior distribution estimated with the joint
dep (solid) and miss ind (dashed) model for 5 randomly-selected simulated data
(colors).

Here there is a notable difference between the results from the two models:
the posterior distributions from the misspecified model that assumes
independent data are less variable than the ones derived using the MC
approximation of the joint dependent likelihood.

This result was expected, since the misspecified model, by assuming
independent data, accounts for more information than is contained in the data.
This leads to an overconfidence that can be detected in the underestimation of
the posterior variance. Results are confirmed by the proportion of datasets
for which the pairwise differences are less than or equal to 0 for each of the
parameters (Table 2). When this quantity is close to 0.5, the variances of the
estimates obtained with the two methods are similar within datasets; when this
quantity is close to 1 it suggests that the variance of the estimates obtained
using the misspecified independent likelihood is systematically larger than
the variance of the estimates obtained with the joint likelihood; and when
this quantity is close to 0 it highlights that the former variance is
systematically smaller then the latter.

Table 2: Proportion of datasets in which the pairwise difference of variance is smaller or equal to 0 for the three transmission parameters. Parameter | Proportion of ${\textsc{pwd}(\text{Var})}\leq 0$
---|---
| Small dependence | Large dependence
${\beta}$ | 0.392 | 0
${\pi}$ | 0.390 | 0
${\iota}$ | 0.378 | 0

The results strongly suggest a systematic difference in variability between
the two methods in the large dependence scenario. In all the simulated
datasets the variances of the posterior distributions of the parameters are
smaller in the analysis using the misspecified independent model than in the
approximate joint model.

####  4.2.2 Comparison for transmission and severity parameters

The same comparison is carried out in a context where inference is drawn both
for the transmission and the severity parameters. Here, since more quantities
are estimated and due to the high correlation of the parameters of epidemic
models, a difference between the results from the two models may be more
difficult to spot. Moreover, in this multi-parameter context, convergence is
sometimes compromised, particularly in the large-dependence scenario. In the
distribution of the pairwise differences, neither for the transmission
parameters nor for the newly estimated severity parameters, can a large
difference be seen (figures reported in the Supplementary Information).

For the large dependence scenario, the only notable difference can be seen in
the distribution of ${}^{\textsc{h}}\theta^{\textsc{ic}}$: the parameter that
links the two datasets, since it defines the probability of IC admission
conditional on hospitalization. When the two datasets are jointly analysed,
they both contribute to the estimation of
${}^{\textsc{h}}\theta^{\textsc{ic}}$, with hospital data informing the
Binomial size in Equation 10 and IC data informing the proportion of people in
the more-severe state. When the two datasets are considered independently, the
hospital data do not play any role in the inference of
${}^{\textsc{h}}\theta^{\textsc{ic}}$.The proportions of pairwise differences
less than or equal to 0 confirm this: the variance of the posterior sample of
the parameter ${}^{\textsc{h}}\theta^{\textsc{ic}}$ is always lower when
inference is drawn with the approximation to the joint dependent likelihood
compared to when the misspecified independent model is adopted (Table 3).

Table 3: Proportion of datasets in which the pairwise difference of variance is smaller or equal to 0 for the transmission and severity parameters.  Parameter | Proportion of ${\textsc{pwd}(\text{Var})}\leq 0$
---|---
| Small dependence | Large dependence
${\beta}$ | 0.498 | 0.554
${\pi}$ | 0.464 | 0.546
${\iota}$ | 0.348 | 0.202
${}^{0}\theta^{\textsc{h}}$ | 0.466 | 0.566
${}^{\textsc{h}}\theta^{\textsc{ic}}$ | 0.778 | 1

####  4.2.3 Influential parameters

As a final comparison, a further investigation into the main cause of the
difference is undertaken. Starting from the small-dependence scenario, one at
a time, each parameter of Table 1 is allowed to take the larger value in
simulating the 500 datasets.

Estimates of the five parameters are then obtained according to the
misspecified independent and the joint dependent model. The posterior
distributions and the plots of the precision statistics are reported in the
Supplementary Information. While a detectable difference in the results is
observed when all the parameters affecting the level of dependence vary, the
same cannot be said when each parameter increases alone. Differences are less
evident, with the probability of detection in IC being the most influential
parameter, as shown in Table 4, where each column corresponds to a scenario
where all the parameters but the header of the column are assumed small.

Table 4: Proportion of datasets in which the pairwise difference of variance is smaller or equal to 0 for the transmission and severity parameters in the scenario with small dependence except for the respective column-name parameter. Increased Parameter |  ${}^{0}\theta^{\textsc{h}}$ |  ${}^{\textsc{h}}\theta^{\textsc{ic}}$ |  $\zeta^{\textsc{h}}$ |  $\zeta^{\textsc{ic}}$
---|---|---|---|---
Parameter | Proportion of ${\textsc{pwd}(\text{Var})}\leq 0$
$\beta$ | 0.468 | 0.454 | 0.296 | 0.490
$\pi$ | 0.450 | 0.454 | 0.214 | 0.496
$\iota$ | 0.342 | 0.082 | 0.052 | 0.032
${{}^{0}\theta^{\textsc{h}}}$ | 0.476 | 0.458 | 0.290 | 0.464
${{{}^{\textsc{h}}\theta^{\textsc{ic}}}}$ | 0.682 | 0.940 | 0.072 | 0.994

##  5 Case-study: influenza during the 2017/18 season

The UK Health Security Agency (UKHSA), formerly Public Health England (PHE),
routinely collects several sources of data to monitor influenza cases. Many of
these datasets have been used separately to provide information on the
transmission (Birrell et al. (2011), Baguelin et al. (2013), Corbella et al.
(2018)) and severity (Presanis et al., 2014) of influenza; nevertheless, a
joint analysis of all the sources has never been performed. Here a joint model
is fitted to multiple data from the 2017/18 influenza season, with the aim of
retrospectively characterising both severity and transmission, while
appropriately attributing stochasticity to the various processes. Among other
datasets, the USISS is analysed and the methodology proposed in Section 3 is
used to jointly exploit hospitalisation and IC unit admissions data.

###  5.1 Data

The following datasets are included in the analysis:

  * •

Confirmed influenza cases collected from week 40 of 2017 to week 20 of 2018
via the USISS comprising the weekly count of IC admissions from, in principle,
all trusts in England (Health Protection Agency, 2011a) and the weekly count
of hospitalizations at all levels of care in a stratified sentinel sample of
the trusts (Health Protection Agency, 2011b);

  * •

Daily counts of General Practitioners consultations for influenza-like illness
(ILI) from a sample of GP monitored by EMIS (Harcourt et al., 2012) or The
Phoenix Partnership (The Phoenix Partnership, 2013);

  * •

Observed proportion of respiratory swabs taken on a selected subset of GP
patients consulting for ILI that test positive for influenza (virological
positivity, Health Protection Agency (2014)).

  * •

Results from cross-sectional serological surveys that inform on the presence
of antibodies in the population (Charlett, 2018).

###  5.2 Model Specification

Figure 5 provides an illustration of the data-generating processes which
includes spread of the virus (transmission), probability of mild symptoms and
severe outcomes (severity), background ILI cases and detection.

The model is specified in discrete time; intervals of length 1 day are denoted
by the index $u=1,2,\dots$, while 1-week intervals are indexed by
$t=1,2,\dots$ (e.g. daily GP consultations are denoted by
$y^{\textsc{g}}_{u}$, while weekly hospitalizations are denoted by
$y^{\textsc{h}}_{t}$). The remaining notation is similar to that of Section 4,
with, e.g., ${}^{\textsc{j}}_{r}X_{s}^{\textsc{k}}$ being the number of people
that experienced event j at time ${r}$ and subsequently event k at time $s$.

Figure 5:  The transmission model classifies the population into susceptible
($S$), exposed ($E_{1},E_{2}$), infectious ($I_{1},I_{2}$) and removed ($R$).
The severity model defines the occurrence of: mild flu cases
${}^{0}X^{\textsc{f}}$ with probability ${}^{0}\theta^{\textsc{f}}$, hospital
cases ${}^{0}X^{\textsc{h}}$ with probability ${}^{0}\theta^{\textsc{h}}$, and
of IC cases conditional on hospitalization ${}^{\textsc{h}}X^{\textsc{ic}}$
with probability ${}^{\textsc{h}}\theta^{\textsc{ic}}$. The detection process
links cases to data defining the probability of reporting GP consultations
($\zeta_{u}^{\textsc{g}}$), hospitalizations ($\zeta_{u}^{\textsc{h}}$) and IC
admissions ($\zeta_{u}^{\textsc{ic}}$), respectively. The background process
models the non-influenza ILI cases, $X^{\textsc{b}}_{t}$.

Some essential elements of the model are outlined here; the complete
description of the model and the derivations of the data-distributions are
provided in the Supplementary Information.

####  5.2.1 Transmission and first severity layer

Denote by $\xi_{u}^{0}$ the number of new infections generated during day $u$.
A deterministic $SEIR$ transmission model, is assumed so that $\xi_{u}^{0}$ is
a function of the parameters $\pi,\iota,\beta,\sigma,\gamma,\kappa$ ,
representing the proportion of individuals initially immune; the proportion of
initially infected/infectious individuals; the transmission rate; the rate of
becoming infectious; the recovery rate; and the school-closure effect,
respectively.

The infection processes of individuals who will experience hospital
admissions, ${}^{0}_{u}X^{\textsc{h}}$, and influenza-related GP
consultations, ${}^{0}_{u}X^{\textsc{f}}$, are assumed to follow a time non-
homogeneous Poisson process, i.e.:

$\begin{split}\left({}^{0}_{u}X^{\textsc{h}}\bigg{|}\xi^{0}_{u},{}^{0}\theta^{\textsc{h}}\right)&\sim\text{Pois}\left({}^{0}\theta^{\textsc{h}}\cdot\xi^{0}_{u}\right)\qquad\text{for
}u=0,1,\dots,U\\\
\left({}^{0}_{u}X^{\textsc{f}}\bigg{|}\xi^{0}_{u},{}^{0}\theta^{\textsc{f}}\right)&\sim\text{Pois}\left({}^{0}\theta^{\textsc{f}}\cdot\xi^{0}_{u}\right)\qquad\text{for
}u=0,1,\dots,U\end{split}$ (13)

with ${}^{0}\theta^{\textsc{f}}$ and ${}^{0}\theta^{\textsc{h}}$ denoting the
probability of being visited by a GP and being admitted to hospital,
respectively.

####  5.2.2 GP-consultations

The inhomogeneous Poisson process of the people who will become mildly
symptomatic and consult a GP, ${}^{0}_{u}X^{\textsc{f}}$ in Equation 13, is
the main component of the likelihood of the GP consultations data. This
process too is affected by delays that can be modelled by a discrete r.v.
modelling the days elapsing between infection and consulting a practice.

In addition to these individuals, background, non-influenza cases appear in GP
consultation data; these endemic cases of other respiratory viruses and
bacterial infections often follow a yearly seasonality, peaking around the
same time as the seasonal influenza epidemic (Paul, Held and Toschke, 2008).
The background seasonality pattern is modelled by a weekly-varying sine-cosine
oscillation, similar to Held, Höhle and Hofmann (2005).

The general process describing the number of people consulting a GP with ILI
symptoms is then obtained by adding the endemic and epidemic processes.
Virological data are used to disentangle the proportion of the former process
out of the total.

Lastly, a Binomial emission is again chosen to model the observational
process. However, the probability of attending a GP practice is subject to
weekly fluctuations, caused by the weekend closure of GP practices, hence a
day-of the week distortion effect is included in the probability of detecting
a GP consultation.

####  5.2.3 Hospitalization and IC admissions

A model for dependent data on hospitalizations and IC-admissions data is
illustrated in Section 3. Here too the joint likelihood is factorised in the
marginal distribution of the IC admissions $Y_{1:T}^{\textsc{ic}}$, from the
Poisson distribution of Equation 11, and the distribution of the
hospitalizations $Y_{1:T}^{\textsc{h}}$ conditionally on IC data approximated
via MC integration as proposed in Algorithm 1.

###  5.3 Results

The joint distribution of the unknown parameter vector (including
transmission, severity, background and detection parameters) is derived via
GIMH. Weakly-informative priors are assumed for most of the parameters of
interest, all the prior distributions used and the estimated posteriors are
reported in the Supplementary Information.

The model presents a fair fit to the data (Figure 6) with all the observations
being included in the 95% CrIs of the posterior predictive distributions of
the GP consultation, virology and hospitalization data. The model, however,
struggles to fit IC data which might be in conflict with other sources of
information in the model. These data, unlike GP and hospital data, don’t
suggest a second peak in infection in the latest weeks of winter.

Figure 6: Median and 95% CrIs (green) for the posterior predicted distribution
of: GP data (a); virological data (b); hospitalizations (c) and IC admissions
(d). Red points are the observed data.

The model provides a comprehensive picture of influenza season 2017/18.
Transmission, in terms of daily number of new infections is inferred and
background ILI cases are also described (Figure 7). Moreover, all the key
severity parameters of interest are well identified and informed by the joint
use of the data, including the case-hospitalization risk (Median 0.0032, 95%
CrI 0.0022-0.0049) and the hospital-IC admission risk (Median 0.0667, 95% CrI
0.0574-0.078).

(a) Median (red) and 95% CrIs (green) of the daily number of new infections.
The grey areas corresponding to periods of lower-transmission (school
holidays). 20 randomly selected trajectories are also plotted as thin green
lines.

(b) Median (solid line) and 95% CrIs (shaded area) of the prior (red) and
posterior (green) for $b_{u}$, the mean of the rate of the daily number of
non-influenza ILI GP consultations.

Figure 7: Posterior results from the joint analysis of 2017/18 data:
transmission and background.

The model in all presents many levels of sophistication. Firstly, ILI cases
attributable to influenza are disentangled from the background ILI, which,
while not being directly related to influenza, contributes to pressure on GP
practices. Secondly, the model accounts for peculiarities of the GP data, such
as the day-of-the-week variation, uncovering the true underlying transmission
and severity process. Lastly and more importantly, exact inference from
dependent data via MC methods enables the simultaneous use of information from
hospital and IC surveillance.

##  6 Discussion

This work contributes to the state-of-the-art literature on the analysis of
epidemic data, from a methodological/computational perspective, providing new
tools, applicable to a wide class of problems that includes dependent data;
and from an application perspective, providing innovative results on the
analysis of influenza from multiple sources.

###  6.1 Methodological and computational perspective

Section 2 contains a concise review of SSMs and their analysis, particularly
focussed on the analysis of multiple epidemic data. There are other reviews of
these methods in the literature (Kantas et al. (2015); Schön et al. (2018)),
some of which even target epidemic applications (McKinley et al. (2014);
Dureau, Ballesteros and Bogich (2013); Ionides, Bretó and King (2006)).
However, these reviews mainly address temporal dependences and, to the
knowledge of the authors, little is known about the dependence across other
domains. In the context of epidemic analysis, where temporal dependence can be
easily overcome through realistic deterministic approximations of transmission
dynamics, it is key to consider the sources of stochasticity and dependence
linked to other processes (in this case the severity process) and its
implications.

Sections 3 and 4 embrace this challenge and highlight the possible
complexities of such an analysis. Starting from a specific example, where
there is dependence between data sources, an approach to the estimation of the
static parameters of the system is proposed. Key aspects such as overlaps
between datasets and delays between consecutive events are then revealed to
affect results. It is shown via simulation that simplifying assumptions, which
assume independence among data and avoid expensive integration, lead to
overconfident results and misleadingly narrow interval estimates.

Our general and comprehensive way to set up epidemic models are paired with a
specific estimation routine for their static parameters: pseudo-marginal
methods. These methods are robust and simple to set up since they rely on the
standard MH algorithms, however they are expensive in two respects: firstly, a
high number of simulations of the hidden states is required at each iteration
of the static parameter to approximate its likelihood; and, secondly, the
parameter space might be explored inefficiently due to their random walk
behaviour. The first of these aspects might be addressed, for example, by
converting to a data-augmented (DA) MCMC perspective (O’Neill and Roberts,
1999), which only needs a single sample from the hidden states per MCMC
iteration. DA approaches however, require tailored implementations for each
different system considered, and therefore do not provide a general recipe
applicable to a wide class of models. In terms of exploration of the parameter
space, the unavailability of closed-form likelihood prevents the use of
popular gradient based methods (e.g. Hamiltonian Monte Carlo (HMC) (Neal et
al., 2011) or the more recent Piecewise Deterministic Markov Process (PDMP)
sampler (Bierkens, Fearnhead and Roberts, 2019)). Alternatively, a natural way
to improve the estimation routine would be to consider sequential Monte Carlo
(SMC) methods on the parameter space, which in this context can be thought of
as SMC2 (Chopin, Jacob and Papaspiliopoulos, 2013).

###  6.2 Application perspective

The ability of estimating key transmission and severity parameters from
epidemic data has never been more crucial. With the COVID-19 epidemic
affecting almost every aspect of people’s life, evidence-based decision-making
has become an imperative aspect of public health policy. Uncertainty
quantification around central estimates is key to understand both the extreme
possible scenarios that could present in the near future and what further
information is needed to increase knowledge and inform decisions.

Section 5 of the paper proposes a complex joint analysis of multiple data on
influenza where both mild and severe cases inform transmission and severity
parameters. The model presents multiple levels of sophistication including:
the use of dependent data; the use of data from both non-confirmed and
confirmed cases; the disentangling of background ILI cases; the presence of
delays between events. The key innovation of this model is that it matches a
deterministic transmission dynamic with a stochastic severity and reporting
process. In this sense it differentiates from other works that make joint use
of multiple data but assume deterministic severity dynamics (e.g. Keeling et
al. (2020)).

Shubin et al. (2016) proposes a similar analysis of the 2009 pandemic in
Finland. Even though his model presents many interesting aspects (e.g. the
inclusion of both environmental and demographic stochasticity), unrealistic
simplifications are made to allow for estimation in reasonable time, including
the exclusion of an exposed compartment, the discretization of time into
extremely-large intervals (one week), and the assumption of no time elapsing
between events of increasing severity. We consider our approach as a valid
alternative to Shubin et al. (2016), especially on a large population such as
that of England, where the deterministic assumption is even more likely to
hold.

Going forward, the model could be extended in, at least, two interesting ways:
firstly, population heterogeneity could be considered, through the formulation
of an age-specific model that makes use of contact patterns between age
groups; and secondly, more overdispersion could be allowed in the severity
process by, for example, replacing Poisson r.v.s with Negative Binomial r.v.s.
These changes would almost certainly improve the fitting of the model to data.

Lastly, other interesting analyses could be performed starting from the
proposed model. The analysis presented here is retrospective and aimed at a
posterior estimation of key epidemiological parameters; it would be
interesting to test the predictive ability of such a model so that it can be
used in real time for decision-making on resource allocation for different
healthcare facilities. Moreover, more insights could be gained from the
quantification of the contribution to inference provided by each separate
dataset from a value of information perspective Jackson et al. (2019).

SSM

    state-space model
r.v.

    random variable
POMP

    partially observed Markov process
HMM

    hidden Markov model
BPF

    bootstrap particle filter
MC

    Monte Carlo
HMC

    Hamiltonian Monte Carlo
SMC

    sequential Monte Carlo
MH

    Metropolis Hastings
MCMC

    Monte Carlo Markov chain
GIMH

    grouped independence Metropolis Hastings
MCWM

    Monte Carlo within Metropolis
DAG

    directed acyclic graph
USISS

    UK Severe Influenza Surveillance System
NHS

    National Health Service
IC

    Intensive Care
CrI

    Credible Interval
PHE

    Public Health England
ILI

    influenza-like illness
GP

    General Practitioner
RCGP

    Royal College of General Practitioners
ONS

    Office of National Statistics
UKHSA

    UK Health Security Agency
DA

    data-augmented
PDMP

    Piecewise Deterministic Markov Process

## Acknowledgements

The author would like to thank colleagues at Public Health England,
particularly Dr Richard Pebody, Dr Andre Charlett and Dr Nikolaos
Panagiotopoulos, for providing the data and prior information for the analysis
in Section 6. We are grateful to Dr Michail Shubin for insightful discussions
on the earlier work on dependent data streams. Lastly, Dr Trevelyan J McKinley
and Dr Paul Kirk provided useful feedback on an earlier version of this work
for which we are thankful.

AC is supported by Bayes4Health EPSRC Grant EP/R018561/1. All four authors
were supported by the MRC, programme grant MC_UU_00002/11.

Supplementary Information This report contains further work on the study of
dependence between data and supplementary information to the study of Section
5. The former includes further comments on SSM inference, the derivations of
Algorithm 1 and the results form the simulation study of Section 4. The latter
includes extensive explanation of the model, inferential methods and results
of the analysis. Supplementary code and data Code for the simulation study in
Section 4 and the analysis of Section 5 are available at
https://github.com/alicecorbella/EpiDependentData. The data analysed in
Section 5 are based on routine health-care data, which cannot be made
available to others by the study authors. Requests to access these non-
publicly available data are handled by the Public Health England Office for
Data Release and its successor at the UK Health Security Agency.

## References

  * Health Protection Agency (2011a) [author] Health Protection Agency (2011a). Sources of UK flu data—Influenza Surveillance in the United Kingdom. Available at the following link.
  * Health Protection Agency (2011b) [author] Health Protection Agency (2011b). UK Severe Influenza Surveillance System (USISS) - Protocol for sentinel Acute NHS Trusts 2011-12. Available at the following link.
  * Health Protection Agency (2014) [author] Health Protection Agency (2014). Sources of UK Flu Data: Influenza Surveillance in the United Kingdom. Available at the following link.
  * Anderson (1996) [author] Anderson, Jeffrey LJ. L. (1996). A method for producing and evaluating probabilistic forecasts from ensemble model integrations. Journal of climate 9 1518–1530. doi: 10.1175/1520-0442(1996)009. 10.1175/1520-0442(1996)009¡1518:AMFPAE¿2.0.CO;2
  * Andersson and Britton (2012) [author] Andersson, HakanH. and Britton, TomT. (2012). Stochastic epidemic models and their statistical analysis 151\. Springer Science & Business Media.
  * Andrieu, Doucet and Holenstein (2010) [author] Andrieu, ChristopheC., Doucet, ArnaudA. and Holenstein, RomanR. (2010). Particle markov chain monte carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72 269–342.
  * Andrieu and Roberts (2009) [author] Andrieu, ChristopheC. and Roberts, Gareth OG. O. (2009). The pseudo-marginal approach for efficient Monte Carlo computations. The Annals of Statistics 37 697–725. doi: 10.1214/07-AOS574.
  * Arulampalam et al. (2002) [author] Arulampalam, M SanjeevM. S., Maskell, SimonS., Gordon, NeilN. and Clapp, TimT. (2002). A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on signal processing 50 174–188. doi: 10.1109/78.978374.
  * Baguelin et al. (2013) [author] Baguelin, MarcM., Flasche, StefanS., Camacho, AntonA., Demiris, NikolaosN., Miller, ElizabethE. and Edmunds, W JohnW. J. (2013). Assessing optimal target populations for influenza vaccination programmes: an evidence synthesis and modelling study. PLoS medicine 10 e1001527.
  * Beaumont (2003) [author] Beaumont, Mark AM. A. (2003). Estimation of population growth or decline in genetically monitored populations. Genetics 164 1139–1160.
  * Bierkens, Fearnhead and Roberts (2019) [author] Bierkens, JorisJ., Fearnhead, PaulP. and Roberts, GarethG. (2019). The zig-zag process and super-efficient sampling for Bayesian analysis of big data. The Annals of Statistics 47 1288–1320.
  * Birrell, De Angelis and Presanis (2018) [author] Birrell, Paul JP. J., De Angelis, DanielaD. and Presanis, Anne MA. M. (2018). Evidence synthesis for stochastic epidemic models. Statistical science: a review journal of the Institute of Mathematical Statistics 33 34\. doi: 10.1214%2F17-STS631.
  * Birrell et al. (2011) [author] Birrell, Paul JP. J., Ketsetzis, GeorgiosG., Gay, Nigel JN. J., Cooper, Ben SB. S., Presanis, Anne MA. M., Harris, Ross JR. J., Charlett, AndréA., Zhang, Xu-ShengX.-S., White, Peter JP. J., Pebody, Richard GR. G. et al. (2011). Bayesian modeling to unmask and predict influenza A/H1N1pdm dynamics in London. Proceedings of the National Academy of Sciences 108 18238–18243. doi: 10.1073/pnas.1103002108.
  * Birrell et al. (2021) [author] Birrell, PaulP., Blake, JoshuaJ., Van Leeuwen, EdwinE., Gent, NickN. and De Angelis, DanielaD. (2021). Real-time nowcasting and forecasting of COVID-19 dynamics in England: the first wave. Philosophical Transactions of the Royal Society B 376 20200279\. doi: 10.1098/rstb.2020.0279.
  * Boddington, Verlander and Pebody (2017) [author] Boddington, NLN., Verlander, NQN. and Pebody, RGR. (2017). Developing a system to estimate the severity of influenza infection in England: findings from a hospital-based surveillance system between 2010/2011 and 2014/2015. Epidemiology & Infection 145 1461–1470.
  * Bretó et al. (2009) [author] Bretó, CarlesC., He, DaihaiD., Ionides, Edward LE. L. and King, Aaron AA. A. (2009). Time series analysis via mechanistic models. The Annals of Applied Statistics 319–348.
  * Britton (2010) [author] Britton, TomT. (2010). Stochastic epidemic models: a survey. Mathematical biosciences 225 24–35.
  * Brockwell et al. (2016) [author] Brockwell, Peter JP. J., Brockwell, Peter JP. J., Davis, Richard AR. A. and Davis, Richard AR. A. (2016). Introduction to time series and forecasting. Springer.
  * Brookmeyer et al. (1994) [author] Brookmeyer, RonR., Gail, Mitchell HM. H., Gail, Mitchell HM. H. et al. (1994). AIDS epidemiology: a quantitative approach. Oxford University Press on Demand.
  * Brooks et al. (2011) [author] Brooks, SteveS., Gelman, AndrewA., Jones, GalinG. and Meng, Xiao-LiX.-L. (2011). Handbook of markov chain monte carlo. CRC press.
  * Charlett (2018) [author] Charlett, AndreA. (2018). Personal communication on the level of immunity priori 2017/18 Influenza season.
  * Chopin, Jacob and Papaspiliopoulos (2013) [author] Chopin, NicolasN., Jacob, Pierre EP. E. and Papaspiliopoulos, OmirosO. (2013). SMC2: an efficient algorithm for sequential analysis of state space models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75 397–426.
  * Churchill (2005) [author] Churchill, Gary AG. A. (2005). Hidden M arkov Models. Encyclopedia of Biostatistics 4\.
  * Corbella (2019) [author] Corbella, AliceA. (2019). Statistical inference in stochastic/deterministic epidemic models to jointly estimate transmission and severity, PhD thesis, University of Cambridge doi: 10.17863/CAM.41539.
  * Corbella et al. (2018) [author] Corbella, AliceA., Zhang, Xu-ShengX.-S., Birrell, Paul JP. J., Boddington, NickiN., Pebody, Richard GR. G., Presanis, Anne MA. M. and De Angelis, DanielaD. (2018). Exploiting routinely collected severe case data to monitor and predict influenza outbreaks. BMC public health 18 1–11.
  * Diekmann, Heesterbeek and Britton (2012) [author] Diekmann, OdoO., Heesterbeek, HansH. and Britton, TomT. (2012). Mathematical tools for understanding infectious disease dynamics. Princeton University Press.
  * Dukic, Lopes and Polson (2012) [author] Dukic, VanjaV., Lopes, Hedibert FH. F. and Polson, Nicholas GN. G. (2012). Tracking epidemics with Google flu trends data and a state-space SEIR model. Journal of the American Statistical Association 107 1410–1426.
  * Dureau, Ballesteros and Bogich (2013) [author] Dureau, JosephJ., Ballesteros, SébastienS. and Bogich, TiffanyT. (2013). SSM: Inference for time series analysis with State Space Models. arXiv preprint arXiv:1307.5626.
  * Funk et al. (2018) [author] Funk, SebastianS., Camacho, AntonA., Kucharski, Adam JA. J., Eggo, Rosalind MR. M. and Edmunds, W JohnW. J. (2018). Real-time forecasting of infectious disease dynamics with a stochastic semi-mechanistic model. Epidemics 22 56–61. doi: 10.1016/j.epidem.2016.11.003.
  * Harcourt et al. (2012) [author] Harcourt, SES., Smith, GEG., Elliot, AJA., Pebody, RichardR., Charlett, AndreA., Ibbotson, SS., Regan, MM. and Hippisley-Cox, JJ. (2012). Use of a large general practice syndromic surveillance system to monitor the progress of the influenza A (H1N1) pandemic 2009 in the UK. Epidemiology & Infection 140 100–105.
  * Hastings (1970) [author] Hastings, W KeithW. K. (1970). Monte Carlo sampling methods using Markov chains and their applications.
  * Held, Höhle and Hofmann (2005) [author] Held, LeonhardL., Höhle, MichaelM. and Hofmann, MathiasM. (2005). A statistical framework for the analysis of multivariate infectious disease surveillance counts. Statistical modelling 5 187–199.
  * Ionides, Bretó and King (2006) [author] Ionides, Edward LE. L., Bretó, CarlesC. and King, Aaron AA. A. (2006). Inference for nonlinear dynamical systems. Proceedings of the National Academy of Sciences 103 18438–18443.
  * Jackson et al. (2019) [author] Jackson, ChristopherC., Presanis, AnneA., Conti, StefanoS. and De Angelis, DanielaD. (2019). Value of information: Sensitivity analysis and research design in Bayesian evidence synthesis. Journal of the American Statistical Association. doi: 10.1080/01621459.2018.1562932.
  * Kalman (1960) [author] Kalman, Rudolph EmilR. E. (1960). A new approach to linear filtering and prediction problems.
  * Kantas et al. (2015) [author] Kantas, NikolasN., Doucet, ArnaudA., Singh, Sumeetpal SS. S., Maciejowski, JanJ. and Chopin, NicolasN. (2015). On particle methods for parameter estimation in state-space models. Statistical science 30 328–351.
  * Keeling and Rohani (2011) [author] Keeling, Matt JM. J. and Rohani, PejmanP. (2011). Modeling infectious diseases in humans and animals. Princeton university press.
  * Keeling et al. (2020) [author] Keeling, Matt JM. J., Dyson, LouiseL., Guyver-Fletcher, GlenG., Holmes, AlexA., Semple, Malcolm GM. G., Investigators, ISARIC4CI., Tildesley, Michael JM. J. and Hill, Edward ME. M. (2020). Fitting to the UK COVID-19 outbreak, short-term forecasts and estimating the reproductive number. Statistical Methods in Medical Research. doi: 10.1101/2020.08.04.20163782.
  * Kermack and McKendrick (1927) [author] Kermack, William OgilvyW. O. and McKendrick, Anderson GA. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character 115 700–721.
  * King, Nguyen and Ionides (2015) [author] King, Aaron AA. A., Nguyen, DaoD. and Ionides, Edward LE. L. (2015). Statistical inference for partially observed Markov processes via the R package pomp. arXiv preprint arXiv:1509.00503.
  * Kingman (1992) [author] Kingman, John Frank CharlesJ. F. C. (1992). Poisson processes 3\. Clarendon Press.
  * Lekone and Finkenstädt (2006) [author] Lekone, Phenyo EP. E. and Finkenstädt, Bärbel FB. F. (2006). Statistical inference in a stochastic epidemic SEIR model with control intervention: Ebola as a case study. Biometrics 62 1170–1177. doi: 10.1111/j.1541-0420.2006.00609.x.
  * Lindsten and Schön (2013) [author] Lindsten, FredrikF. and Schön, Thomas BT. B. (2013). Backward simulation methods for Monte Carlo statistical inference. Foundations and Trends® in Machine Learning 6 1–143.
  * Lindsten and Schön (2017) [author] Lindsten, FredrikF. and Schön, Thomas BT. B. (2017). Learning of dynamical systems. Particle filters and Markov chain methods.
  * Magpantay et al. (2016) [author] Magpantay, FMGF., De Celles, M DomenechM. D., Rohani, PP. and King, AAA. (2016). Pertussis immunity and epidemiology: mode and duration of vaccine-induced immunity. Parasitology 143 835–849.
  * McKinley et al. (2014) [author] McKinley, Trevelyan JT. J., Ross, Joshua VJ. V., Deardon, RobR. and Cook, Alex RA. R. (2014). Simulation-based Bayesian inference for epidemic models. Computational Statistics & Data Analysis 71 434–447.
  * Metropolis et al. (1953) [author] Metropolis, NicholasN., Rosenbluth, Arianna WA. W., Rosenbluth, Marshall NM. N., Teller, Augusta HA. H. and Teller, EdwardE. (1953). Equation of state calculations by fast computing machines. The journal of chemical physics 21 1087–1092.
  * Neal et al. (2011) [author] Neal, Radford MR. M. et al. (2011). MCMC using Hamiltonian dynamics. Handbook of markov chain monte carlo 2 2\.
  * Nishiura (2011) [author] Nishiura, HiroshiH. (2011). Real-time forecasting of an epidemic using a discrete time stochastic model: a case study of pandemic influenza (H1N1-2009). Biomedical engineering online 10 1–16. doi: 0.1186/1475-925X-10-15.
  * O’Neill and Roberts (1999) [author] O’Neill, Philip DP. D. and Roberts, Gareth OG. O. (1999). Bayesian inference for partially observed stochastic epidemics. Journal of the Royal Statistical Society: Series A (Statistics in Society) 162 121–129.
  * The Phoenix Partnership (2013) [author] The Phoenix Partnership (2013). Real-time Syndromic Surveillance.
  * Paul, Held and Toschke (2008) [author] Paul, MichaelaM., Held, LeonhardL. and Toschke, André MA. M. (2008). Multivariate modelling of infectious disease surveillance data. Statistics in medicine 27 6250–6267.
  * Presanis et al. (2009) [author] Presanis, Anne MA. M., De Angelis, DanielaD., New York City Swine Flu Investigation Team, Hagy, AngelaA., Reed, CarrieC., Riley, StevenS., Cooper, Ben SB. S., Finelli, LynL., Biedrzycki, PaulP. and Lipsitch, MarcM. (2009). The severity of pandemic H1N1 influenza in the United States, from April to July 2009: a Bayesian analysis. PLoS medicine 6 e1000207.
  * Presanis et al. (2013) [author] Presanis, Anne MA. M., Ohlssen, DavidD., Spiegelhalter, David JD. J. and De Angelis, DanielaD. (2013). Conflict diagnostics in directed acyclic graphs, with applications in Bayesian evidence synthesis. Statistical Science 28 376–397. doi: 10.1214/13-STS426.
  * Presanis et al. (2014) [author] Presanis, Anne MA. M., Pebody, Richard GR. G., Birrell, Paul JP. J., Tom, Brian DMB. D., Green, Helen KH. K., Durnall, HayleyH., Fleming, DouglasD. and De Angelis, DanielaD. (2014). Synthesising evidence to estimate pandemic (2009) A/H1N1 influenza severity in 2009–2011. The Annals of Applied Statistics 2378–2403. doi: 10.1214/14-AOAS775.
  * Schön et al. (2018) [author] Schön, Thomas BT. B., Svensson, AndreasA., Murray, LawrenceL. and Lindsten, FredrikF. (2018). Probabilistic learning of nonlinear dynamical systems using sequential Monte Carlo. Mechanical systems and signal processing 104 866–883. doi: 10.1016/j.ymssp.2017.10.033.
  * Shubin et al. (2016) [author] Shubin, MikhailM., Lebedev, ArtemA., Lyytikäinen, OutiO. and Auranen, KariK. (2016). Revealing the true incidence of pandemic A (H1N1) pdm09 influenza in Finland during the first two seasons—An analysis based on a dynamic transmission model. PLoS computational biology 12 e1004803.
  * Solin et al. (2018) [author] Solin, ArnoA., Kok, ManonM., Wahlström, NiklasN., Schön, Thomas BT. B. and Särkkä, SimoS. (2018). Modeling and interpolation of the ambient magnetic field by Gaussian processes. IEEE Transactions on robotics 34 1112–1127. doi: 10.1109/TRO.2018.2830326.
  * UK Health Security Agency (UKHSA) (2021) [author] UK Health Security Agency (UKHSA) (2021). Sources of UK flu data—Influenza Surveillance in the United Kingdom. Available at the following link.
  * Weidemann et al. (2014) [author] Weidemann, FelixF., Dehnert, ManuelM., Koch, JudithJ., Wichmann, OleO. and Höhle, MichaelM. (2014). Bayesian parameter inference for dynamic infectious disease modelling: rotavirus in Germany. Statistics in Medicine 33 1580–1599. doi: 10.1002/sim.6041.

  *[SSM]: state-space model
  *[SSMs]: state-space model
  *[r.v.s]: random variable
  *[r.v.]: random variable
  *[BPF]: bootstrap particle filter
  *[MCMC]: Monte Carlo Markov chain
  *[MH]: Metropolis Hastings
  *[GIMH]: grouped independence Metropolis Hastings
  *[MCWM]: Monte Carlo within Metropolis
  *[DAG]: directed acyclic graph
  *[USISS]: UK Severe Influenza Surveillance System
  *[NHS]: National Health Service
  *[IC]: Intensive Care
  *[MC]: Monte Carlo
  *[CrI]: Credible Interval
  *[UKHSA]: UK Health Security Agency
  *[PHE]: Public Health England
  *[ILI]: influenza-like illness
  *[GP]: General Practitioner
  *[CrIs]: Credible Interval
  *[DA]: data-augmented
  *[HMC]: Hamiltonian Monte Carlo
  *[PDMP]: Piecewise Deterministic Markov Process
  *[SMC]: sequential Monte Carlo


same map $f_\infty = f$ that realizes the lamination $L$. Since
$\mathcal{MPS}(S, K)$ is compact, this means that $f_i \to f$ in the
topology on $\mathcal{MPS}(S, K)$. This means that $f_i \to f$ in the
compact-open topology, and also that the hyperbolic metrics $d_i$
induced by $f_i$ converge in $\teich(S)$ to the hyperbolic metric $d$
induced by $f$.

Let $N \cong F \times \RR$ be a cusped quasi-Fuchsian $3$–manifold,
as in the statement of the theorem. Let $L$ be the pleating lamination
on the lower core boundary $\bdy_- \core(N)$. The bending measure
$\mu$ on $\bdy_- \core(N)$ is a transverse measure of full support, so
$L$ is a measurable lamination. By Lemma <ref>,
there is a sequence of ideal triangulations $\tau_i \to L'$ in the
Chabauty topology, where $L \subset L'$. Now, by Lemma
<ref>, the pleating maps $f_i \from F \to N$ that
pleat along $\tau_i$ induce hyperbolic metrics on $F$ that converge in
$\teich(F)$ to the induced metric on $\bdy_- \core(N)$. By the definition of
$\mathcal{MPS}(S, N)$, these pleating maps also converge in the
compact–open topology to a pleating map for $\bdy_- \core(N)$.

§ A LEMMA IN POINT-SET TOPOLOGY

Let $X$ be a topological space, and let $\mathcal{S}$ be a cover of
$X$ by closed sets. Define a discrete walk of length $k$
through sets of $\mathcal{S}$ to be a sequence of points $x_0, x_1,
\ldots, x_k$, such that $x_0, x_1 \in S_1$, $x_1, x_2 \in S_2$, and so
on, for sets $S_1, \ldots, S_k$ that belong to $\mathcal{S}$. We say
this sequence is a walk from $x_0$ to $x_k$.

The following observation is needed in the proof of
Lemma <ref>.

Let $X$ be a connected topological space, and let $S_1, \ldots, S_n$
be closed sets whose union is $X$. Then, for any pair of points $x,y
\in X$, there is a discrete walk from $x$ to $y$ through sets in the
collection $\{ S_1, \ldots, S_n \}$.

Define an equivalence relation $\equiv$ on $X$, where

\[
x \equiv y \quad \Leftrightarrow \quad \mbox{there exists a discrete
walk from $x$ to $y$ through $\{ S_1, \ldots, S_n \}$.}
\]

Reversing and concatenating walks proves this is an equivalence
relation.  Furthermore, each closed set $S_i$ must be entirely
contained in an equivalence class, because its points are connected by
a walk of length $1$. Thus each equivalence class is closed. Since the
connected space $X$ cannot be expressed as a disjoint union of
finitely many closed sets, all of $X$ must be in the same equivalence