Bayesian state space modelling for COVID-19: with Tennessee and New York case studies

12/30/2020
by   Adam Spannaus, et al.
0

We develop a Bayesian inferential framework for the spread of COVID-19 using mechanistic epidemiological models, such as SIR or SEIR, and allow the effective contact rate to vary in time. A novel aspect of our approach is the incorporation of a time-varying reporting rate accounting for the initial phase of the pandemic before testing was widely available. By varying both the reporting rate and the effective contact rate in time, our models can capture changes in the data induced by external influences, such as public health intervention measures, for example. We view COVID-19 incidence data as the observed measurements of a hidden Markov model, with latent space represented by the underlying epidemiological model, and employ a particle Markov chain Monte Carlo (PMCMC) sampling scheme for Bayesian inference. Parameter inference is performed via PMCMC on incidence data collated by the New York Times from the states of New York and Tennessee from March 1, 2020 through August 30, 2020. Lastly, we perform Bayesian model selection on the different formulations of the epidemiological models, make predictions from our fitted models, and validate our predictions against the true incidence data for the week between August 31, 2020 and September 7, 2020.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 11

page 13

06/10/2020

Semiparametric Bayesian Inference for the Transmission Dynamics of COVID-19 with a State-Space Model

The outbreak of Coronavirus Disease 2019 (COVID-19) is an ongoing pandem...
06/18/2020

Sparse HP Filter: Finding Kinks in the COVID-19 Contact Rate

In this paper, we estimate the time-varying COVID-19 contact rate of a S...
03/10/2021

Bayesian sequential data assimilation for COVID-19 forecasting

We introduce a Bayesian sequential data assimilation method for COVID-19...
07/21/2021

Tracking the Transmission Dynamics of COVID-19 with a Time-Varying Coefficient State-Space Model

The spread of COVID-19 has been greatly impacted by regulatory policies ...
09/12/2020

Tracking disease outbreaks from sparse data with Bayesian inference

The COVID-19 pandemic provides new motivation for a classic problem in e...
06/04/2018

A Bayesian Penalized Hidden Markov Model for Ant Interactions

Interactions between social animals provide insights into the exchange a...
05/25/2021

Efficient Bayesian model selection for coupled hidden Markov models with application to infectious diseases

Performing model selection for coupled hidden Markov models (CHMMs) is h...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

During a disease outbreak, epidemiological forecasts help inform public health decisions by making predictions about how widely the disease will spread (Astolfi et al., 2012). This information is critical for public health officials trying to understand the dynamics through which a disease is transmitted and to slow its propagation (Anderson et al., 1992). For example, specific understanding of the transmission dynamics helps to quantify a quarantine period for those exposed to the disease. Moreover, in the case of an emerging health crisis, there are uncertainties in the process of disease transmission such as the rate of infection or incubation period of a novel disease. These considerations are especially relevant during the current SARS-CoV-2 (COVID-19) pandemic as the uncertainty surrounding the disease’s transmission and incubation rates have made coordinating the global response to this disease a substantial challenge to public health officials.

Real-time infectious disease forecasts are important tools utilized by public health officials to analyze, characterize, and predict disease propagation (Myers et al., 2000). One of the most widely employed mathematical models for general disease transmission is the Susceptible-Infected-Removed (SIR) model of Kermack and McKendrick (1927). The SIR model describes the interactions between three population groups: those who are are susceptible to an infection, those who are infected, and those who have been removed from the population through recovery or death. Over time, the SIR model has been augmented to include other compartments such as an exposed state between the susceptible and infected compartments, a compartment for asymptomatic carriers, or a compartment for individuals immune to the infection (Blackwood and Childs, 2018). In this vein and specifically related to the COVID-19 pandemic, it is thought that there is a latency period in the course of the disease (Rothe et al., 2020) during which an individual carrying the disease does not present symptoms, suggesting the possibility that a model of COVID-19 should include an ‘exposed’ state and be described using a Susceptible-Exposed-Infected-Resistant (SEIR) model.

Deterministic models cannot reflect the real-time implications of preventative measures, such as safer-at-home orders, social distancing, or mask mandates that directly impact the transmission rate (Andersson and Britton, 2012)

. To capture the effects of preventative measures when modelling the COVID-19 outbreak, we choose stochastic state space representations of the SIR and SEIR models. To this end, we adopt a Bayesian approach to modelling the progression of disease in a population through the different compartments in each epidemiological model, and allow the rate of transmission to vary as a function of time, thus capturing the spread of a disease throughout a community. The proposed stochastic models view the effective transmission rate as a time-dependent random variable and account for under-reporting of cases, heterogeneity of disease outcomes, i.e., asymptomatic cases and transmission, and mitigation strategies used to slow the spread of the disease, such as shelter-in-place orders or mask mandates. Furthermore, we consider two different scenarios in the reporting of active cases; such a parameter allows the model to account for differences in availability of testing and advances in testing accuracy.

To demonstrate our Bayesian approach, we study two regional outbreaks with different mitigation strategies using the COVID-19 case counts as collated by the New York Times (Times, 2020). Considering data from Tennessee and New York showcases the flexibility of our modelling strategy as the spread of COVID-19 in these states presents different modelling challenges. Indeed, each state handled the beginning of the pandemic differently, employed different testing strategies, and continue to have different reopening policies (Raifman et al., 2020)

. Notably, these two states experienced the initial six months of the pandemic in widely differing fashion. Such widely varying approaches are demonstrated by the COVID-19 data sets associated Tennessee and New York which exhibit drastically different dynamics of case counts. Our Bayesian state space model coupled with particle Markov chain Monte Carlo (PMCMC) sampling can readily quantify uncertainty in parameter estimates and predict short term the extent of disease spread when fitted to each of the two data sets. Our experiments highlight the capacity of this type of modelling approach to fit data representing different dynamics of disease spread, to provide a mechanistic framework to analyze the propagation of COVID-19, and to provide a principled methodology for post-hoc model analysis and criticism.

The outline of the paper is as follows. In Section 2 we detail the epidemiological and Bayesian models that we employ. The sequential Monte Carlo (SMC) and PMCMC methods are described in Section 3 with the associated pseudocode. Numerical results of fitting the epidemological models to the data, MCMC diagnostics, and short-term predictions are in Section 4. Lastly, we conclude and give discussion in Section 5.

2 Modelling

2.1 Epidemiological models

For the problem of modelling the COVID-19 pandemic, we will compare the SIR and SEIR epidemiological models. The salient difference between these is the inclusion of the ‘exposed’ state in the latter; individuals move into this the compartment once they have been exposed to an infected individual, but are not yet actively infected themselves. Moreover, such a compartment is pertinent when quantifying the breadth of the pandemic due to the challenges in correctly accounting for such individuals. The ‘exposed’ compartment not only acts as a delay between the ‘susceptible’ and ‘infected’ compartments, but individuals in the ‘exposed’ bin may transmit the disease before becoming actively infected themselves. We augment these epidemiological models with an additional parameter representing the rate at which active cases are reported. In the pandemic’s early phases, wide-spread access to testing was unavailable, and the disease prevalence within a community was widely under-reported (Fisher and Wilder-Smith, 2020; Lau et al., 2020). Thus we consider two different scenarios with the reporting rate: (i) a constant rate spanning the entire time duration from March 1, 2020 through August 30, 2020, or (ii) a reporting rate that varied between the first month of data that we consider, March 1 through March 28, and all subsequent reported incidence data, March 29 through Sept 7. The time period from March 1 to March 28 encompasses the initial phase of the pandemic, when testing was less reliable and less available (Arevalo-Rodriguez et al., 2020; Lau et al., 2020).

Defining the characteristics of an on-going pandemic presents a challenging problem for a deterministic epidemiological model. We take the approach similar to those of Del Moral and Murray (2015); Dureau et al. (2013); Funk et al. (2018) and allow the transmission rate to vary with time. Thus we may capture intervention measures enacted by public health officials, such as mask mandates or shelter-in-place orders, or ‘super-spreader’ that have a direct impact on disease prevalence. By accommodating the transmission via a time-varying transmission rate, we may better quantify the spread of the disease, which in turn, yields better information to public health officials concerned with slowing community spread.

The SIR epidemiological model is represented by the system of differential equations

(2.1)

where represents COVID-19 susceptible individuals, COVID-19 infected individuals, and removed individuals who are no longer able to become infected. Susceptible individuals move into the infected compartment at rate . Infected individuals become removed at rate . Lastly, we introduce a stochastic process via a stochastic differential equation (SDE) that defines the effective contract rate, thus controlling the extent to which a disease spreads throughout a population. Accordingly, for the SDE we define drift and diffusion terms, both parameterized by , and a function .

An alternative to the SIR model is the SEIR model which includes an additional ‘exposed’ compartment and corresponding parameter . This parameter controls the movement of those individuals who are in the disease’s latent period into an active infection, and thus capable of infecting other susceptible individuals. The SEIR model is written

(2.2)

In this study, we take and , thus assuming that the transmission rate on any day is likely to be the same as the previous day, and set . By choosing and , the resulting path , defined by the SDE , is Brownian motion. Allowing to vary in time defines the effective contact rate , which controls the extent to which a disease spreads between individuals within a population.

2.2 Hidden Markov model

In the present section we describe the hidden Markov model (HMM) methodology used in our Bayesian epidemiological modelling. Briefly, an HMM is a bivariate discrete time process such that the latent states represent a Markov chain. The observations , are independent random variables conditioned on , such that the conditional distribution of only depends on . An HMM may be written as a stochastic dynamical system

for real functions and sequences of i.i.d. random variables and that are also independent of . In the present context of an epidemiological model we discretize the system of ODEs in time and write or depending on the model under consideration. Making some regularity assumptions, described by Zhigljavsky and Zilinskas (2007) on the evolution dynamics of the system, given by eqs. 2.12.2, we may associate a state evolution density, parameterized by , with the evolution of the states in an epidemiological model. This density, coupled with an initial state, , fully characterizes the evolution of states in a HMM. Accounting for measurement errors, partial and modified observations, we introduce an observation density, parameterized by .

Denoting any sequence for as , we view the stochastic epidemiological model as an HMM, written

(2.3)

for an observation density . We define , or , depending on the model under investigation, and . We require the diffusion process to have a unique solution, i.e., an invariant distribution that may be viewed as a prior over . Conditions for existence of such a distribution are described in Oksendal (2013).

2.3 Bayesian formulation

To place the epidemiological models in eqs. 2.12.2 in the context of a discrete time HMM, we discretize the path of a Brownian motion, so at each time step we sample and set , where is the diffusion term of the SDE. Note that for a fixed , the state evolution density is a deterministic function of the model parameters and initial conditions , which may be computed from an approximation to the intractable integral solution of the ODE system. This point is important in the context of the sampling schemes that we present in the following section. For if this was not the case, we would be unable to estimate the marginal likelihood, a necessary step in the PMCMC method used to sample from our posterior.

For the observations in our HMM, we assume that not all instances of the disease are reported, due to lack of available testing, asymptomatic cases, or changes in the reporting of active cases. Consequently we consider two different modelling scenarios: (i) the rate at which active cases are reported remains constant for the duration of the collected data, or (ii) active cases were reported at differing rates during time intervals and

. Mathematically, the probability of successfully reporting an active case is

for case (i) and

for case (ii). For simplicity of notation, we will write as the reporting rate unless explicitly describing a particular modelling scenario. We define the observed weekly incidence rate as with for the SIR model, or with for the SEIR model. Precisely, denotes the number of new cases transitioning into the infected compartment in week .

As the reporting of each new case may be seen as an independent Bernoulli random variable, i.e., each case is reported with probability

, and the waiting time until the first successful report is geometrically distributed with the same parameter

. Thus we are interested in the probability of , conditional on the number of successful reported cases, . Since the latter is the sum of i.i.d. geometrically distributed random variables with parameter , it then follows that

. Invoking the Central Limit Theorem, these observations are approximately normal with mean

and variance

. The variance contains an additional parameter, , which describes the overdispersion within a population. It is a common occurrence in count data (Bolker, 2008; Breslow, 1984) and signals heterogeneity within a population. Furthermore, it accounts for large variances between individual outcomes, and can show the presence of ‘super-spreading’ events (Lloyd-Smith, 2007).

To fully specify our Bayesian model, we define and denote the prior of the parameters by . Then the posterior distribution of the transmission rate and model parameters , given the observed incidence data is

(2.4)

Following the standard decomposition we may write our posterior distribution, eq. 2.4, as

(2.5)

which suggests a sampling strategy of alternating between sampling from densities and via the PMCMC algorithm of Andrieu et al. (2010). Such an approach alleviates issues of convergence and insufficient exploration of the sample space that can arise due to correlations and dependencies between variables.

Sampling from the posterior defined in eq. 2.4, or factorization thereof eq. 2.5, allows us to infer the time-varying transmission rate and summary statistics about the reporting rate, and to make predictions about the future course of the pandemic. Moreover, by posing the epidemiological model as a Bayesian inference problem, we may readily quantify the parameter uncertainty and predictive uncertainty of the model.

2.4 Bayesian model selection

To compare the four different epidemiological models considered herein in a principled fashion, we employ Bayesian model selection. For each state we have the following models: (i) SIR with a constant reporting rate, (ii) SIR with a varying reporting rate, (iii) SEIR with a constant reporting rate, and (iv) SEIR with a varying reporting rate; we denote these by and

, respectively. Accordingly, we investigate the conditional posterior probability

of model given data .

Kass and Raftery (1995)

define the Bayes factor

associated with models and as

(2.6)

Computing the Bayes factors which are subsequently used to find the posterior probability of model conditioned on the observations , given by

(2.7)

where

is the ratio of prior probabilities of

against . In the present setting we set and therefore , due to not having prior knowledge about which model may be more preferable.

The question of model choice is then answered by choosing the model with the maximum posterior probability as computed by eq. 2.7. Indeed, each model’s posterior probability is a weighted average over all Bayes factors, and maximizing such a probability yields the most likely model to explain the data.

3 Markov chain Monte Carlo

3.1 Particle Markov chain Monte Carlo

To sample from the posterior density of eq. 2.4, we employ PMCMC sampling (Andrieu et al., 2010). We describe the algorithmic procedure and give the details necessary for the implementation of our model; for an in-depth discussion and theoretical results, see Andrieu et al. (2010); Del Moral et al. (2006). First, the sequential Monte Carlo (SMC) procedure is described, followed by PMCMC.

SMC algorithms (Del Moral et al., 2006) provide a methodology for sampling from distributions defined by state-space models. According to the decomposition of our posterior distribution in eq. 2.5, samples are first drawn from the conditional . The SMC methodology, outlined in algorithm 1, yields a sequence of densities that approximates and marginal distributions for a given and . The procedure first yields an approximation of and , by drawing samples from an importance density (Andrieu et al., 2010), for particles . The algorithm then approximates and , for all subsequent iterations by sampling from importance densities . By requiring these to be of the form

, one readily computes an unbiased estimate of the marginal likelihood

, which is necessary for the Metropolis-Hastings acceptance ratio in the PMCMC sampler (Del Moral et al., 2006). The SMC algorithm yields this approximation by noting that

which follows from the conditional independence between observations and we have defined

which is an MCMC estimate at time of

1:Initialize , , and sample particles
2:for  do
3:     for  do
4:         Sample
5:         Solve for from
6:         Set
7:         Compute the normalized weights:
8:     end for
9:     Sample
10:end for
11:Compute marginal likelihood
Algorithm 1 Sequential Mote Carlo (SMC)

Having obtained samples from via SMC, it remains to sample from . Recalling the posterior factorization in eq. 2.5, we look for a proposal distribution of the form . At each iteration of the PMCMC algorithm, a value of the parameter is proposed, followed by a sample generated by the SMC procedure. Thus the problem of sampling from is reduced to sampling from , as samples from are obtained via the SMC algorithm.

1:Initialize
2:Sample using algorithm 1 (SMC)
3:Compute marginal likelihood
4:for  do
5:     Sample from proposal density
6:     Sample
7:     Compute using algorithm 1 (SMC)
8:     Compute the acceptance probability
9:     Sample uniformly at random from the unit interval
10:     if  then
11:         Set
12:         Set
13:     else
14:         Set
15:         Set
16:     end if
17:end for
Algorithm 2 Particle Marginal Metropolis–Hastings (PMCMC)

3.2 MCMC diagnostics

We use two MCMC diagnostics to monitor the quality of convergence and of mixing, namely the potential scale reduction factor (PSRF) and the effective sample size (ESS). Multivariate versions of PSRF and ESS are employed, introduced by Brooks and Gelman (1998) and Vats et al. (2019), respectively. The multivariate PSRF and ESS require knowledge of the parameter posterior covariance matrix, which is estimated by the multivariate initial sequence estimator (MINSE) estimator of Dai and Jones (2017).

Multivariate PSRF and ESS estimators provide two advantages over their univariate counterparts. Firstly, they take into account the covariance structure of the parameter space instead of accounting for the parameter variance only. Secondly, they provide single number summaries across all parameters instead of yielding univariate estimates for each parameter. Typical thresholds for the univariate versions are and (Vehtari et al., 2020). Since there are no widely established thresholds for multivariate PSRF and ESS, we employ the empirical thresholds known from the univariate scenario.

3.3 Uncertainty quantification

For each model, we quantify the parameter and predictive uncertainty. To quantify the former, we utilize estimates of the parameter posterior covariance matrix, Monte Carlo standard error (MCSE), posterior correlation matrix and credible intervals. To quantify predictive uncertainty we make use of predictive credible intervals.

The posterior covariance matrix of the parameters of each model is estimated using the MINSE of Dai and Jones (2017). The MCSE of the parameters equals the square root of the diagonal of the MINSE estimator. The posterior correlation matrix of the parameters is estimated by standardizing the relevant entries of the MINSE estimator with the associated MCSEs.

To compute credible intervals for the parameter posteriors and predictive posteriors, we first estimate each posterior density at each time step via a kernel density estimator (KDE)

(Silverman, 1986)

. Having computed an estimate of the density, we then construct an empirical cumulative distribution function at each time step. The desired quantiles for the credible intervals are then readily calculated from the empirical distribution functions.

4 Experiments

4.1 Overview of data and of experimental setup

We have presented four different models for each state: (i) SIR with a constant reporting rate, (ii) SIR with a varying reporting rate, (iii) SEIR with a constant reporting rate, and (iv) SEIR with a varying reporting rate. Here we apply the PMCMC methodology to incidence data of the COVID-19 pandemic for the states of New York and Tennessee employing each of these four models. The data used in our experiments were obtained from the New York Times COVID data repository (Times, 2020)

, and fit our models on data collected through August 30. The data is the cumulative weekly case count for each state. We chose to model new cases on a weekly, as opposed to daily, basis for computational considerations. Indeed, each particle in the ensemble requires the numerical approximation of a system of non-linear ordinary differential equations with

time steps, at each MCMC iteration. This computational burden cannot be ignored, as a sufficiently large ensemble of particles and MCMC iterations are required to ensure that the diagnostic criteria we employ in Section 4.2 are satisfied. We then make predictions on the continued spread of the disease for the following week and compare our predictions with the reported case incidence obtained from the same database.

4.1.1 Data description

The differences between the initial phases of the pandemic in New York and Tennessee present challenges in modelling such diverse data sets. Indeed, the reported case data from these states present the spread of COVID-19 from two different perspectives. The pandemic went through a period of sustained exponential growth in New York, most notably in the New York City metropolitain region, then abated to a near constant level for subsequent months. The data for Tennessee mirrors that of most of the nation. We see a slow initial increase, followed by a first wave in April and a much larger increase in July. Furthermore, during the initial wave in New York, testing was not widely available and their 7-day rolling average of positive test case peaked at nearly 50% in early April (of Health, a). In fig. 1 we plot the incidence data from Tennessee and New York highlighting the variations in the two data sets.

Figure 1: Weekly new case summaries of Tennessee (blue) and New York (orange) over the first 6 months of the pandemic.

4.1.2 pMCMC setup

The model parameters , and or

, depending on the model, are given wide uninformative priors, due to the uncertainty about the ongoing pandemic and disparities in reporting data. We model the infection period as a Gaussian distribution with mean of 5.058 days and standard deviation of 1.51 that is truncated to have upper and lower bounds of 11.8 and 2.228 days respectively

(Lauer et al., 2020). The prior for the latent period is obtained from the study of Moghadas et al. (2020)

, and is modeled as a gamma distribution with shape and scale parameters, 1.058 and 2.174 respectively

(He et al., 2020). For the initial proportions of the population in states we chose a Dirichlet distribution, while constraining the mean of to be , and let the means of the other compartments be equal. By this choice, we ensure that the condition or is satisfied, depending on the model under consideration. Thus, the sum over all compartments in the epidemiological model at each time step is the same as the total population . Lastly, we ran PMCMC sampling with particles and obtain samples from the posterior after a burn-in period of iterations.

4.1.3 Software

For our modelling application, we use the Bayesian modelling software libbi, (Murray, 2013) and R packages rbi and rbi.helpers (Funk, 2016; Jacob et al., 2020). For MCMC diagnostics, we use the kanga Python package: https://github.com/papamarkou/kanga. Our models, data, and code for running the inference may be found at https://github.com/aspannaus/Covid-model.

4.2 MCMC diagnostics

The multivariate PSRF and ESS laid out in section 3.2 are used to monitor convergence and sampling effectiveness of the PMCMC simulations. These MCMC diagnostic results from fitting each of the four models to both data sets are presented in table 1. Four chains are run for each of the eight settings, and the reported PSRF for each setting is computed across the four chains. The reported ESS for each setting is the mean of the ESS, averaged across all four chains of the setting.

The multivariate PSRF is less than for all eight eight settings. Increasing the number of particles in PMCMC reduces the PSRF but increases the computational cost of PMCMC. ESS is above the threshold of in all eight settings with varying sampling effectiveness.

Model Diagnostic ESS PSRF TN SIR 3879 1.0004 2614 1.0102 SEIR 3113 1.0109 2526 1.0008 NY SIR 935 1.0005 1226 1.0275 SEIR 10026 1.0046 8060 1.0100

Table 1: MCMC diagnostics: ESS and PSRF.

4.3 Bayesian model selection

Bayesian model selection is performed based on the Bayes factors eq. 2.6 and associated model posterior probabilities eq. 2.7 presented in section 2.4. The results of model selection are presented in table 2.

SIR SIR SEIR SEIR New York 0.259 0.730 8.851-5 1.798-2 Tennessee 0.199 0.776 7.237-2 1.775-2

Table 2: Bayesian factors for the four candidate models.

Two main conclusions are drawn from the results of this anaylsis. Firstly, the model posterior probabilities indicate a strong preference for SIR over SEIR models for both Tennessee and New York data sets considered here. This is anticipated, since the used data sets contain only reported active cases without any information about the number of exposed individuals. Thus, the number of exposed individuals appearing in SEIR modelling is a latent variable to be learnt from limited amount of data rather than an observed variable from which information could had been elicited.

Secondly, we see that it is preferable to use a variable reporting rate both with the New York data and with the Tennessee data. The preference of models with variable reporting rate can be explained by two facts, namely the change in the size of the pandemic and the change of testing infrastructure. The early stage in New York was characterized by exponential growth, and relatively low transmission after the peak in April. Tennessee on the other hand, did not experience such a rapid increase in cases during the first months. During the first months of the pandemic, testing infrastructure was not well developed, and it was difficult for even actively sick individuals to get tested. Indeed, in Tennessee, only 42,000 tests were administered in the month of March, as compared with nearly 894,000 in July per the of Health (b).

4.4 Uncertainty quantification

4.4.1 Parameter uncertainty

(a) New York
(b) Tennessee
Figure 4: Estimates of the contact rate for (a) New York and (b) Tennessee from our model. The mean transmission rate is given by the black line, the 75%, and 95% credible intervals are denoted by the dark and light grey shaded areas respectively.

To quantify parameter uncertainty, the MCSE of the parameter posterior mean estimates is displayed in table 3 (see section 3.3 for its derivation). For each parameter in our posterior, the reported MCSE is the mean of MCSEs across all four chains.

Tennessee New York SIR SEIR SIR SEIR 34.88 35.16 19.27 27.81 11.76 14.07 12.12 13.83 21.53 17.14 6.42 9.52 0.15 0.21 0.73 0.67 0.27 0.25 0.48 0.71 15.47 16.34 9.60 16.72 9.16 9.63 15.71 16.85 16.54 17.37 7.78 11.51 1.55 1.64 1.54 1.44 2.31 2.54 0.58 1.08 4.09 4.61 5.09 4.35 5.62 5.61 1.67 2.55 1-9 1.93 0.65 1.57 0.67 3.39 1.01 0.94 0.58 0.75 0.94 1.54 0.60 1-9 0.49 0.63 0.70 0.65 0.42 0.48 0.11 0.17 3.28 3.58 3.31 3.07 3.27 3.58 0.54 0.77 4.16 4.56 3.90 4.15 6.26 6.71 1.26 2.21 5.82 6.42 5.13 5.57 10.50 10.51 2.38 4.02 8.22 8.98 10.50 10.86 13.90 15.41 5.49 16.50 11.33 12.17 13.70 14.36 18.63 20.22 26.35 94.68 8.99 9.71 7.97 8.78 20.48 19.89 420.69 410.67 7.45 8.26 7.19 7.73 25.02 23.50 652.01 642.43 11.13 12.34 20.64 19.32 32.12 28.61 543.20 839.92 12.37 13.59 13.95 14.53 31.27 27.55 427.42 487.86 10.04 11.06 9.96 11.05 28.52 25.03 496.65 455.42 10.28 11.19 11.08 12.16 27.12 23.79 312.24 604.22 8.95 9.74 8.19 8.91 25.77 23.47 69.57 116.02 9.13 10.02 10.11 10.89 27.37 24.39 75.91 52.96 9.89 10.79 9.87 10.81 26.37 23.57 13.92 25.26 8.62 9.45 8.09 9.01 22.87 20.71 8.25 11.39 7.38 8.02 6.57 7.27 19.74 18.10 6.12 8.15 7.66 8.51 7.40 8.07 18.74 17.64 6.92 9.41 8.55 9.28 8.25 8.95 17.75 16.66 5.15 7.50 8.91 9.63 8.92 9.45 17.21 16.62 6.27 8.32 10.09 11.01 12.09 12.34 19.57 18.18 9.28 11.88 10.97 12.02 13.15 13.67 19.36 18.12 6.94 9.41 11.82 12.73 15.21 15.67 18.31 17.23 6.15 8.29 12.19 12.97 13.99 14.63 18.94 17.99 8.25 10.68 11.20 12.02 12.29 12.84 18.64 17.85 6.09 8.74 10.84 11.58 15.70 15.46 17.56 16.44 16.53 20.57

Table 3: Monte Carlo standard errors (MCSEs) for the parameter posterior mean estimates across different models and data sets.

The MCSE of parameter , which controls the rate at which individuals move from the exposed to infected compartments, is high across both data sets and across SEIR models with constant and varying reporting rate. SIR modelling does not include and therefore does not incur an associated parameter uncertainty.

The MSCE of constant reporting rate is higher than the the MSCEs of the two reporting rates across both SIR and SEIR models and across both data sets. So, a constant reporting rate comes with higher Monte Carlo variability, and from this point of view induces higher uncertainty, than a time-varying reporting rate.

Figure 9 visualizes heatmaps of the posterior correlation matrices for the SIR models with varying reporting rate and for the SEIR models with constant reporting rate. Overall, common posterior correlation structure emerges across the two different data sets.

It appears that there is strong correlation among for various in all models and for both data sets. The SIR model with time-varying reporting rate, which is the most prevalent model arising from Bayesian model selection, appears to have nearly identical and very strong correlation structure in each of the two data sets (fig. (a)a and fig. (b)b). Strong correlation among for various is also present in the SEIR model with constant reporting rate, yet this correlation structure is less preserved across the two data sets (fig. (c)c and fig. (d)d).

(a) TN, SIR, varying
(b) NY, SIR, varying
(c) TN, SEIR, constant
(d) NY, SEIR, constant
Figure 9: Heat map of mean posterior correlation matrix among the parameters of each model and for each dataset.

4.4.2 Predictive uncertainty

The results of fitting SIR with varying reporting rate, the most probable model from our model selection process, are in figs. (b)b(a)a. These plots show the observations in blue, taken every 7 days, with 75% and 95% credible intervals. Note the decreases in cases in the model around day 30 or so, which is approximately one week after Tennessee’s ‘safer-at-home’ order was issued on March 31st. Correspondingly, the transmission rate also decreases at this time, but begins to increase around day 50. These plots corroborate the mobility data from Aliyu, Barkin, Buntin, Clayton, Edwards, Griffin, Grijalva, and McPheeters (Aliyu et al.); Unacast (2020), which show a lag between decreases/increases in mobility and a corresponding change in the incidence and transmission rate.

We ran our model with data through August 30, and made predictions about the number of cases one week into the future and plotted the mean and median of the predictive densities against the actual number of reported cases. These results are in figs. (d)d(c)c. In both the New York and Tennessee data, we see that the median of our predictions exhibit good agreement with the ground truth, as does the mean in the case of both SIR models.

We also note that while the credible intervals contain the ground truth observations, both the posterior means and the credible intervals of the underlying predictive posteriors are skewed to the left in comparison to the ground truth observations. This agrees with our knowledge of how the virus can spread asymptomatically

(Moghadas et al., 2020). The predictive plots show a greater uncertainty in the Tennessee predictions when juxtaposed against the predictions from the New York data. This result is seemingly counterintuitive, as the New York data exhibit a grater variability than the Tennessee data. However, there is an explanation for the higher predictive uncertainty associated with the Tennessee data; the Tennessee Department of Health changed how they defined an active case on September 3, resulting in a one-day decrease of approximately 20,000 reported cases (of Health, b). Consequently, the predictive distribution is skewed to the left in comparison to the ground truth value, as anticipated after considering the change in the definition of an active case by the Tennessee Department of Health.

(a) Weekly incidence and model based reconstruction of the New York data.
(b) Weekly incidence and model based reconstruction of the Tennessee data.
(c) NY Predictions
(d) TN Predictions
Figure 14: Top: Weekly incidence data from the 2020 COVID-19 pandemic, and our model based reconstruction of the case counts in New York and Tennessee, resp. The weekly counts are shown in blue, the black line denotes the median of our posterior, and its mean is denoted by the red triangles. Light and dark shaded grey areas represent 75% and 95% credible intervals, respectively. Bottom: Predictions on New York and Tennessee data for week starting August 31. The error bars are 25% and 75% credible intervals for the predictive densities estimated via kernel density estimation.

5 Conclusions

We carried out Bayesian state space modelling of COVID-19, combining hidden Markov modelling with SIR and SEIR epidemiological modelling and with PMCMC sampling. Several conclusions arise from our analysis.

Firstly, Bayesian model selection indicates that SIR modelling is a more probable choice than SEIR modelling when only COVID-19 incidence data are available. Secondly, a novelty in our modelling approach has been to include a varying reporting rate. Bayesian model selection suggests that a varying reporting rate leads to models that are more likely to fit and explain COVID-19 incidence data. This conclusion is intuitive, since changes in the reporting rate imply changes in the resulting data, so a model with a varying reporting rate is more likely to fit data affected by changes in reporting procedures. Our recommendation, based on Bayesian model selection, is to use and develop SIR models with varying reporting rate for COVID-19 incidence data in the face of changes in the ways that incidence data are reported.

Thirdly, we provide a Bayesian approach to quantify uncertainty in relevant epidemiological parameters and in predictions, yielding a source of important information to public health officials tasked with assessing the present state and with suggesting mitigation strategies for subsequent weeks. Our one-week ahead predictions are accurate, since relevant credible intervals contain the ground truth. Notably, SIR modelling with a varying reporting rate yields the most tight credible intervals (figs. (d)d(c)c) and therefore the smallest predictive uncertainty among the considered models.

In the future, we plan to develop a more flexible model that allows a different reporting rate at each time step or that uses change point detection to decide flexibly at which time step to change the reporting rate.

Acknowledgments

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable,world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Research was supported by the National Virtual Biotechnology Laboratory, a consortium of DOE national laboratories focused on response to COVID-19, with funding provided by the Coronavirus CARES Act. All numerical experiments were completed employing the computing resources of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory. The latter of which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

References

  • Aliyu, Barkin, Buntin, Clayton, Edwards, Griffin, Grijalva, and McPheeters (Aliyu et al.) Aliyu, M., S. Barkin, M. Buntin, E. W. Clayton, E. Edwards, Kathy, M. Griffin, C. Grijalva, and M. McPheeters. Vanderbilt covid-19 report for tennessee: Statewide hospitalizations and mobility.
  • Anderson et al. (1992) Anderson, R. M., B. Anderson, and R. M. May (1992). Infectious diseases of humans: dynamics and control. Oxford university press.
  • Andersson and Britton (2012) Andersson, H. and T. Britton (2012). Stochastic epidemic models and their statistical analysis, Volume 151. Springer Science & Business Media.
  • Andrieu et al. (2010) Andrieu, C., A. Doucet, and R. Holenstein (2010). Particle markov chain monte carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72(3), 269–342.
  • Arevalo-Rodriguez et al. (2020) Arevalo-Rodriguez, I., D. Buitrago-Garcia, D. Simancas-Racines, P. Zambrano-Achig, R. del Campo, A. Ciapponi, O. Sued, L. Martinez-Garcia, A. Rutjes, N. Low, et al. (2020). False-negative results of initial rt-pcr assays for covid-19: a systematic review. medRxiv.
  • Astolfi et al. (2012) Astolfi, R., L. Lorenzoni, and J. Oderkirk (2012). Informing policy makers about future health spending: a comparative analysis of forecasting methods in oecd countries. Health Policy 107(1), 1–10.
  • Blackwood and Childs (2018) Blackwood, J. C. and L. M. Childs (2018). An introduction to compartmental modeling for the budding infectious disease modeler. Letters in Biomathematics 5(1), 195–221.
  • Bolker (2008) Bolker, B. M. (2008). Ecological models and data in R. Princeton University Press.
  • Breslow (1984) Breslow, N. E. (1984). Extra-poisson variation in log-linear models. Journal of the Royal Statistical Society. Series C (Applied Statistics) 33(1), 38–44.
  • Brooks and Gelman (1998) Brooks, S. P. and A. Gelman (1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics 7(4), 434–455.
  • Dai and Jones (2017) Dai, N. and G. L. Jones (2017). Multivariate initial sequence estimators in markov chain monte carlo.

    Journal of Multivariate Analysis

     159, 184–199.
  • Del Moral et al. (2006) Del Moral, P., A. Doucet, and A. Jasra (2006). Sequential monte carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(3), 411–436.
  • Del Moral and Murray (2015) Del Moral, P. and L. M. Murray (2015). Sequential monte carlo with highly informative observations. SIAM/ASA Journal on Uncertainty Quantification 3(1), 969–997.
  • Dureau et al. (2013) Dureau, J., K. Kalogeropoulos, and M. Baguelin (2013). Capturing the time-varying drivers of an epidemic using stochastic dynamical systems. Biostatistics 14(3), 541–555.
  • Fisher and Wilder-Smith (2020) Fisher, D. and A. Wilder-Smith (2020). The global community needs to swiftly ramp up the response to contain covid-19. The Lancet 395(10230), 1109–1110.
  • Funk (2016) Funk, S. (2016). Rbi. helpers: Rbi helper functions.
  • Funk et al. (2018) Funk, S., A. Camacho, A. J. Kucharski, R. M. Eggo, and W. J. Edmunds (2018). Real-time forecasting of infectious disease dynamics with a stochastic semi-mechanistic model. Epidemics 22, 56–61.
  • He et al. (2020) He, X., E. H. Lau, P. Wu, X. Deng, J. Wang, X. Hao, Y. C. Lau, J. Y. Wong, Y. Guan, X. Tan, et al. (2020). Temporal dynamics in viral shedding and transmissibility of covid-19. Nature medicine 26(5), 672–675.
  • Jacob et al. (2020) Jacob, P. E., A. Lee, L. M. Murray, S. Funk, and S. Abbott (2020). Rbi: R interface to libbi.
  • Kass and Raftery (1995) Kass, R. E. and A. E. Raftery (1995). Bayes factors. Journal of the american statistical association 90(430), 773–795.
  • Kermack and McKendrick (1927) Kermack, W. O. and A. G. McKendrick (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(772), 700–721.
  • Lau et al. (2020) Lau, H., T. Khosrawipour, P. Kocbach, H. Ichii, J. Bania, and V. Khosrawipour (2020). Evaluating the massive underreporting and undertesting of covid-19 cases in multiple global epicenters. Pulmonology.
  • Lauer et al. (2020) Lauer, S. A., K. H. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. Meredith, A. S. Azman, N. G. Reich, and J. Lessler (2020). The incubation period of 2019-ncov from publicly reported confirmed cases: estimation and application. medRxiv.
  • Lloyd-Smith (2007) Lloyd-Smith, J. O. (2007). Maximum likelihood estimation of the negative binomial dispersion parameter for highly overdispersed data, with applications to infectious diseases. PloS one 2(2), e180.
  • Moghadas et al. (2020) Moghadas, S. M., M. C. Fitzpatrick, P. Sah, A. Pandey, A. Shoukat, B. H. Singer, and A. P. Galvani (2020). The implications of silent transmission for the control of covid-19 outbreaks. Proceedings of the National Academy of Sciences 117(30), 17513–17515.
  • Murray (2013) Murray, L. M. (2013). Bayesian state-space modelling on high-performance hardware using libbi. arXiv preprint arXiv:1306.3277.
  • Myers et al. (2000) Myers, M. F., D. Rogers, J. Cox, A. Flahault, and S. I. Hay (2000). Forecasting disease risk for increased epidemic preparedness in public health.  47, 309–330.
  • of Health (a) of Health, N. Y. S. D. Percentage positive results by region dashboard.
  • of Health (b) of Health, T. D. Covid-19 critical indicators.
  • Oksendal (2013) Oksendal, B. (2013). Stochastic differential equations: an introduction with applications. Springer Science & Business Media.
  • Raifman et al. (2020) Raifman, J., K. Nocka, D. Jones, J. Bor, S. Lipson, J. Jay, P. Chan, S. Galea, et al. (2020). Covid-19 us state policy database.
  • Rothe et al. (2020) Rothe, C., M. Schunk, P. Sothmann, G. Bretzel, G. Froeschl, C. Wallrauch, T. Zimmer, V. Thiel, C. Janke, W. Guggemos, et al. (2020). Transmission of 2019-ncov infection from an asymptomatic contact in germany. New England Journal of Medicine 382(10), 970–971.
  • Silverman (1986) Silverman, B. W. (1986). Density estimation for statistics and data analysis, Volume 26. CRC press.
  • Times (2020) Times, N. Y. (2020, September). New york times covid-19 data.
  • Unacast (2020) Unacast (2020, September). Social distancing scoreboard.
  • Vats et al. (2019) Vats, D., J. M. Flegal, and G. L. Jones (2019). Multivariate output analysis for Markov chain Monte Carlo. Biometrika 106(2), 321–337.
  • Vehtari et al. (2020) Vehtari, A., A. Gelman, D. Simpson, B. Carpenter, and P.-C. Bürkner (2020). Rank-normalization, folding, and localization: An improved for assessing convergence of MCMC.
  • Zhigljavsky and Zilinskas (2007) Zhigljavsky, A. and A. Zilinskas (2007). Stochastic global optimization, Volume 9. Springer Science & Business Media.