key: cord-0623310-3k9is4h6
authors: Zhou, Tianjian; Ji, Yuan
title: Semiparametric Bayesian Inference for the Transmission Dynamics of COVID-19 with a State-Space Model
date: 2020-06-10
journal: nan
DOI: nan
sha: 8ccec96e210f679762f90f0e289640cb4084a695
doc_id: 623310
cord_uid: 3k9is4h6

The outbreak of Coronavirus Disease 2019 (COVID-19) is an ongoing pandemic affecting over 200 countries and regions. Inference about the transmission dynamics of COVID-19 can provide important insights into the speed of disease spread and the effects of mitigation policies. We develop a novel Bayesian approach to such inference based on a probabilistic compartmental model using data of daily confirmed COVID-19 cases. In particular, we consider a probabilistic extension of the classical susceptible-infectious-recovered model, which takes into account undocumented infections and allows the epidemiological parameters to vary over time. We estimate the disease transmission rate via a Gaussian process prior, which captures nonlinear changes over time without the need of specific parametric assumptions. We utilize a parallel-tempering Markov chain Monte Carlo algorithm to efficiently sample from the highly correlated posterior space. Predictions for future observations are done by sampling from their posterior predictive distributions. Performance of the proposed approach is assessed using simulated datasets. Finally, our approach is applied to COVID-19 data from four states of the United States: Washington, New York, California, and Illinois. An R package BaySIR is made available at https://github.com/tianjianzhou/BaySIR for the public to conduct independent analysis or reproduce the results in this paper.

Researchers have made substantial efforts to study the transmission dynamics of COVID-19, evaluate the effects of government interventions, and forecast infection and death counts. These works include Aguilar et al. (2020) ; Chen and Qiu (2020) ; Flaxman et al. (2020) ; Giordano et al. (2020) ; Gomez et al. (2020) ; Gu et al. (2020) ; IHME COVID-19 health service utilization forecasting team and Murray (2020) ; Li et al. (2020a,b) ; Pan et al. (2020) ; Sun et al. (2020) ; Wang et al. (2020a,b); Woody et al. (2020) ; Wu et al. (2020) ; Zhang et al. (2020) , among many others. The modeling approaches taken by these works can be broadly categorized into three groups: (i) curve fitting, (ii) compartmental modeling, and (iii) agent-based modeling. Curve fitting approaches fit a curve to the observed number of confirmed cases or deaths. For example, IHME COVID-19 health service utilization forecasting team and Murray (2020) use a Gaussian error function to model the cumulative death rate at a specific location. Compartmental modeling approaches (e.g., Li et al., 2020b) consider a partition of the population into compartments corresponding to different stages of the disease, and characterize the transmission dynamics of the disease by the flow of individuals through compartments. Finally, agent-based modeling approaches (e.g., Gomez et al., 2020) use computer simulations to study the dynamic interactions among the agents (e.g., people in epidemiology) and between an agent and the environment.

In this paper, we develop a novel semiparametric Bayesian approach to modeling the transmission dynamics of COVID-19, which is critical for characterizing disease spread. We aim to address a few issues related to the COVID-19 pandemic. First, we provide estimation of key epidemiological parameters, such as the effective reproduction number of COVID-19. The Bayesian framework allows us to elicit informative priors for parameters that are difficult to estimate due to lack of data based on clinical characteristics of COVID-19, and also offers coherent uncertainty quantification for the parameter estimates. Our second goal is to make predictions about the future trends of the spread of COVID-19 (e.g., future case counts), which will be done by calculating the posterior predictive distributions for the future observations. Although such predictions are technically straightforward, we avoid overinterpretation of the predictions because they rely on extrapolation of highly unpredictable human behaviors and the number of diagnostic tests that will be deployed. Nevertheless, such predictions may be useful for the public and decision makers to understand the trends and future impacts of COVID-19 based on current rates of transmission. We shall see this in our case studies later. Our analysis will be based on a probabilistic compartmental model motivated by the classical susceptible-infectious-recovered model (Kermack and McKendrick, 1927) . We will use data of daily confirmed cases reported by the Center for Systems Science and Engineering at Johns Hopkins University (JHU CSSE) (Dong et al., 2020) . We provide an R package BaySIR, available at https://github.com/tianjianzhou/BaySIR, that can be used to conduct independent analysis of COVID-19 data or reproduce the results in this paper.

The proposed Bayesian approach attempts to improve COVID-19 modeling in at least four aspects. First, we explicitly model the number of undocumented infections, which is only considered by some, but not all, existing works. Due to the potentially limited testing capacity and the existence of pre-symptomatic and asymptomatic COVID-19 cases He et al., 2020) , many infected individuals may not have been detected as having the disease. Therefore, modeling of undocumented infections is essential for accurate inference. Second, we estimate the disease transmission rate via Gaussian process regression (GPR), a semiparametric regression method.

The GPR approach is highly flexible and captures nonlinear and non-monotonic relationships without the need of specific parametric assumptions. Third, we develop a parallel-tempering Markov chain Monte Carlo (PTMCMC) algorithm to efficiently sample from the posterior distribution of the epidemiological parameters, which leads to improvements in convergence and mixing compared to a standard MCMC procedure. We find that standard MCMC cannot produce reliable inference due to poor mixing. Lastly, we rigorously assess our approach through simulation studies, sensitivity analyses, cross-validation and goodness-of-fit tests. Such validations provide insights into the modeling of COVID-19 data, not only for our approach, but also for others based on similar assumptions such as the popular compartmental modeling approaches.

The remainder of the paper is organized as follows. In Section 2, we provide a brief review of the susceptible-infectious-recovered (SIR) compartmental model. In Section 3, we develop a probabilistic state-space model for COVID-19 motivated by the classical SIR model. In Section 4, we present strategies for posterior inference. In Section 5, we carry out simulation studies to assess the performance of our method in estimating the epidemiological parameters. In Section 6, we apply our method to COVID-19 data from six states of the United States (U.S.): Washington, New York, California, Florida, Texas, and Illinois. We conclude with a discussion in Section 7.

We start with a review of the susceptible-infectious-recovered (SIR) model (Kermack and McKendrick, 1927; Weiss, 2013) , a simple type of compartmental model.

The purpose of this review is to introduce the reader to the basics of epidemic modeling and motivate our proposed approach.

Consider a closed population of size N . Here, "closed" means that N does not vary over time. It is a good approximation for a fast-spreading and less fatal pandemic like COVID-19. The SIR model divides the population into the following three compartments:

(S) Susceptible individuals: those who do not have the disease but may be infected;

(I) Infectious individuals: those who have the disease and are able to infect the susceptible individuals;

(R) Recovered/removed individuals: those who had the disease but are then removed from the possibility of being infected again or spreading the disease. Here, the removal can be due to several possible reasons, including death, recovery with immunity against reinfection, and quarantine and isolation from the rest of the population.

At time t (t ≥ 0), denote by S t , I t and R t the numbers of individuals in the S, I and R compartments, respectively, and write V t = (S t , I t , R t ). We have S t + I t + R t ≡ N .

2.1. Deterministic SIR Models. The classical SIR model (Kermack and McKendrick, 1927) describes the flow of people from S to I to R via the following system of differential equations:

Here, β is the disease transmission rate, and α is the removal rate. The rationale behind the first equation in (1) is as follows: suppose each infectious individual makes effective contacts (sufficient for disease transmission) with β others per unit time; therefore, βS/N of these contacts are with susceptible individuals per unit time, and as a result, I infectious individuals lead to a rate of new infections (βS/N ) · I. The third equation in (1) describes that the infectious individuals leave the infective class at a rate of αI. The second equation in (1) follows immediately from the first and third equations. The parameters β and α are determined according to the natural history of the disease. The quantities R 0 = β/α and R e = (βS 0 )/(αN ) are referred to as the basic reproduction number and effective reproduction number, respectively, where S 0 is the initial number of susceptibles at time t = 0.

In some applications, it may be convenient to consider a discrete-time approximation of the differential equations in Equation (1), which can be expressed as follows:

for t = 1, 2, . . .. This discretization replaces the derivatives in Equation (1) by the differences per unit time.

The SIR models given by Equation (1) and (2) are both deterministic models, meaning that their behaviors are completely determined by their initial conditions and parameter values.

2.2. Stochastic SIR Models. The deterministic SIR models are appealing due to their simplicity. However, the spread of disease is naturally stochastic. The disease transmission between two individuals is random rather than deterministic. Therefore, a stochastic formulation of the SIR model may be preferred for epidemic modeling, because it allows one to more readily capture the randomness of the epidemic process.

In a stochastic SIR model, {V t : t ≥ 0} is treated as a stochastic process. A commonly used formulation is as follows (Gibson and Renshaw, 1998; ONeill and Roberts, 1999; Andersson and Britton, 2000) . Suppose that an infectious individual makes effective contacts with any given individual in the population at times given by a Poisson process of rate β/N , and assume all these Poisson processes are independent of each other. Therefore, the expected number of effective contacts made by each infectious individual is β per unit time. Furthermore, suppose each infectious individual remains so (before being removed) for a period of time, known as the infectious period. Lastly, assume that the length of the infectious period for each individual is independent and follows an exponential distribution with mean α −1 . It can be shown that {V t : t ≥ 0} is a Markov process with transition probabilities:

Here, δ is a small increment in time.

2.3. State-space SIR Models. There are, of course, other ways to model the uncertainty of the epidemic process. Probabilistic state-space modeling approaches that build on deterministic models have recently been popular in the statistics literature (Dukic et al., 2012; Osthus et al., 2017 Osthus et al., , 2019 . A state-space SIR model typically consists of two components: an evolution model for the epidemic process, and an observation model for the data. As an example, the model in Osthus et al. (2017) has the form

for t = 1, 2, . . .. In the evolution model, f (V t−1 , β, α) is the solution to Equation (1) at time t with a initial value of V t−1 at time (t − 1) and parameters β and α, and V t is assumed to be centered at f (·) with its variance characterized by κ. In other words, κ measures the derivation of V t from the solution given by the deterministic model. In the observation model,Ĩ t is the number of patients seen with the disease reported by healthcare providers, which can be thought of as a proxy to the true number of infectious individuals I t . The observationĨ t is assumed to be centered at I t with variance characterized by λ. State-space epidemic models are quite flexible and are in general more computationally manageable compared to stochastic epidemic models as in Equation (3).

The SIR model can be extended in many different ways, such as by considering vital dynamics (births and deaths) and demographics, adding more compartments to the model, and allowing more possible transitions across compartments. For example, the susceptible-exposed-infectious-recovered (SEIR) model includes an additional compartment for exposed individuals who are exposed to the disease but are not yet infectious, and the susceptible-infectious-recovered-infectious (SIRS) model allows recovered individuals to return to a susceptible state. These extensions may better capture the characteristics of the disease under consideration. For a comprehensive review of deterministic epidemic models, see, for example, Anderson and May (1991) , Hethcote (2000) or Brauer (2008) . For a comprehensive review of stochastic epidemic models, see, for example, Becker and Britton (1999) , Andersson and Britton (2000) or Allen (2008) .

We now turn to our proposed model for the COVID-19 data, which belongs to the state-space model category (Section 2.3). Our approach integrates the discrete-time deterministic SIR model (Equation 2) and semiparametric Bayesian inference. To capture some unique features of COVID-19, we consider the following extensions of the classical SIR model. First, we split the infectious individuals into two subgroups:

undocumented infectious individuals and documented infectious individuals. The reason is that many people infected with SARS-CoV-2 have not been tested for the virus thus are not detected or reported as having the infection (Li et al., 2020b) .

Second, we allow some epidemiological parameters (such as the disease transmission rate β) to be time-varying to reflect the impact of mitigation policies such as stay-athome orders and the change of public awareness of the disease over time. We discuss details next.

3.1. Model for the Epidemic Process. Consider the transmission dynamics of COVID-19 in a specific country or region (e.g., a state, province or county). For simplicity, we consider a closed population (with no immigration and emigration) and also ignore nature births and deaths. Let N denote the population size. At any time point, we assume that each individual in the population precisely belongs to one of the following four compartments:

(S) Susceptible individuals who do not have the disease but are susceptible to it;

(UI) Undocumented infectious individuals who have the disease and may infect the susceptible individuals. However, they have not been detected as having the disease for several possible reasons. For example, they may have limited symptoms and are thus not tested for the disease;

(DI) Documented infectious individuals who have been confirmed as having the disease and are capable of infecting the susceptible individuals;

(R) Removed individuals who had the disease but are then removed from the possibility of being infected again or spreading the disease.

We further assume that the infectious individuals (including both the UI-and DIindividuals) infect the S-individuals with a transmission rate of β. After being infected, a S-individual first becomes an UI-individual before being detected as a DIindividual. All the infectious (UI-and DI-) individuals recover or die with a removal rate of α. Those UI-individuals who have not been removed are diagnosed with the disease with a diagnosis rate of γ. In total, there are four possible transitions across compartments: S to UI, UI to R, UI to DI, and DI to R. See Figure 1 . Note that it is possible to assume different transmission rates for the UI-and DI-individuals, or to further split the UI and DI compartments into smaller subgroups (e.g., quarantined, hospitalized, etc.) with each subgroup having its distinct transmission rate. It is also possible to consider an extra compartment for the exposed (but not yet infectious) individuals as in the SEIR model. Here, we use a more parsimonious model without the exposed compartment for simplicity and characterize the average transmission rate for all infectious individuals with a single parameter β (β depends on time, which will be clear later). Finally, we assume recovery from COVID-19 confers immunity to reinfection, although there is only limited evidence for this assumption (Long et al., 2020; Kirkcaldy et al., 2020) .

Undocumented We define day t = 0 as the date when the 100th case is confirmed in the country/region under consideration, and index subsequent dates by t = 1, 2, . . . , T , where T is the current date. The reason for choosing day 0 in this way is because we believe the transmission dynamics of the disease is more trackable after a sufficient number of infectious individuals are reported in the country/region, although the choice of "the 100th case" is arbitrary and can be modified. Denote by S t , I U t , I D t and R t the numbers of individuals belonging to compartments S, UI, DI and R on day t, respectively. We have S t + I U t + I D t + R t ≡ N . The transmission rate and diagnosis rate are allowed to vary over time and are hereafter denoted by β t and γ t , respectively.

The number of individuals diagnosed with the disease between day (t − 1) and day t is observed and is denoted by B t−1 . This is our data. We propose modeling the transmission dynamics of COVID-19 over time by the following equations:

The epidemic process, {V t , t = 0, 1, . . . , T }, is determined by its initial value V 0 , the parameters {β t , α}, and the ob-

Rigorously speaking, V t should be a vector of non-negative integers, but for computational convenience, we relax this restriction and only require it to be a vector of non-negative real numbers. Model (4) is a simple extension of (2) by adding a component of I U , the undocumented infections, and by incorporating the observed daily new cases B t−1 into the equations. Later, we introduce a model for the observation B t−1 to complete the state-space model.

With time-varying disease transmission rates, the basic reproduction number and effective reproduction number are also functions of time. That is, R 0 (t) = β t /α and

Here, R e (t) is interpreted as the rate of secondary infections generated by each infectious case at time t, scaled by the length of the infectious period (α −1 ). If R e (t) < 1 for t ≥ t * , then the number of infectious individuals (I U t + I D t ) will monotonically decrease after time t * , because each infectious individual will only be able to infect less than 1 other during the course of his/her infectious period. In other words, an R e (t) < 1 indicates containment of the disease. Due to the important role of R e (t) in characterizing disease spread, we consider the estimation of R e (t) as our main interest.

3.2. Model for the Observed Data. Our observations only consist of the daily new confirmed COVID-19 cases, B t . Assume that on day t, the UI-individuals who have not been removed are diagnosed with the disease with a diagnosis rate of γ t .

where γ t is between 0 and 1. We consider the logit transformation of γ t ,γ

Other transformations, such as the probit and complementary log-log transformations, can also be specified in the BaySIR package. Empirically we find the proposed model to be robust to different specifications of the link function (See Appendix C). We assume a prior transformation,γ

where y t is a vector of covariates that are thought to be related to the diagnosis rate.

In other words, the sampling model for B t can be written as

In the simulation studies and real data analyses, we use a simple choice of y t = 1, assuming the mean diagnosis rate is a constant. It is possible to include other covariates in y t , such as the number of tests, but empirically we find it hard to detect the effects of these covariates. In the BaySIR package, the user has the option to include any covariates. The parameters η and σ 2 γ are the regression coefficients and variance term, respectively, where σ 2 γ captures random fluctuations of confirmed case counts and report errors.

For some countries and regions, the numbers of recoveries and deaths are also available, and one may think of using them as the observed number of removed individuals. We choose not to use these data for two reasons. First, many infected individuals, even with confirmed disease, are not hospitalized, and their recoveries are not recorded. In other words, the reported number of recoveries and deaths is a significant underestimate of the size of the removed population. Second, according to Wölfel et al. (2020) and He et al. (2020) , the ability of a COVID-19 patient to infect others becomes negligible several days before the patient recovers or dies, suggesting that "removal" in our application is not equivalent to "recovery or death".

In what follows, we discuss prior specification for the initial condition and parameters. Due to the limited amount of observable information, many latent variables and parameters in the proposed model are unidentifiable. See Appendix A for a detailed discussion with an example showing that two epidemic processes with distinct parameters lead to exactly the same observed data. We note that this problem is pervasive in most existing methods, and a typical solution to the problem is to prespecify some parameter values based on prior knowledge. Here, we elicit informative priors for some parameters based on the clinical characteristics of COVID-19, which favor more clinically plausible estimates.

Initial condition. The initial condition of the epidemic process refers to the vector

. We assume that there are no removed individuals on day 0, i.e., R 0 = 0. As a result, the number of DI-individuals on day 0, I D 0 , equals to the cumulative number of confirmed cases on that day and is observed. We further assume

where Ga(ν 1 , ν 2 ) refers to a gamma distribution with shape and rate parameters ν 1 and ν 2 , respectively. We set ν 1 = 5 and ν 2 = 1, such that E(I U 0 /I D 0 ) = 5. This choice is based on the findings in Li et al. (2020b) that 86% of all infections were undocumented at the beginning of the epidemic in China. Lastly, note that

The disease transmission rate β t must be non-negative. We considerβ(t) = log(β t ) and assumẽ

where GP[m(t), C(t, t )] refers to a Gaussian process (GP) with mean function m(t) and covariance function C(t, t ). The GP (Rasmussen and Williams, 2006 ) is a very flexible prior model for a stochastic process. It enables one to capture potential nonlinear relationships between t andβ(t) without the need to impose any parametric assumptions. Specifically, for any t 1 , . . . , t n ≥ 0, the vector (β(t 1 ), . . . ,β(t n )) follows a multivariate Gaussian distribution with mean (m(t 1 ), . . . , m(t n )) and covariance matrix C with the (i, j)-th entry being C(t i , t j ). For applications of GP to epidemic modeling, see, for example, Xu et al. (2016) and Kypraios and ONeill (2018) .

We specify m(t) and C(t, t ) as below:

Here, x t is a vector of covariates that are thought to be related to the transmission rate, and µ is a vector of regression coefficients. In the simulation studies and real data analyses, we use x t = (1, t) , which contains an intercept term and the time.

Other covariates, such as indicators for mitigation policies at time t, may also be included in x t . Nevertheless, in practice, we find our GP model with a time trend is sufficient to capture the change ofβ(t) over time and the potential effects of mitigation policies and public awareness. Users of our software may include other covariates using the R package BaySIR. The variance parameter σ 2 β characterizes the amplitude of the difference betweenβ(t) and m(t), and the correlation parameter ρ characterizes the correlation betweenβ(t) andβ(t ) for any t and t . We note that based on our specification of the covariance function, our GP model is equivalent to a first-order autoregressive model. Indeed, autoregressive models of any orders are discrete-time equivalents of GP models with Matérn covariance functions (Roberts et al., 2013) .

We place the following priors on µ, σ β and ρ:

such that E(σ 2 β ) = 0.1 and E(ρ) = 0.8. Here, Inv-Ga(·, ·) refers to an inverse gamma distribution, and Beta(·, ·) refers to a beta distribution. The prior choices for σ 2 β and ρ shrinkβ(t) toward its mean function (i.e., a linear regression model) and

impose a strong prior correlation between the transmission rates for two consecutive days. For the prior of µ, we use µ * = (−1.31, 0) and Σ µ = diag(0.3 2 , 1 2 ), where diag(·) represents a diagonal matrix. In this way, the prior median of the basic reproduction number on day 0 is 2.5 (with 95% credible interval 1.4 to 4.5), assuming the infectious period is 9.3 days. This is based on the findings in Li et al. (2020a) and Wu et al. (2020) . The prior also induces a mild shrinkage (towards 0) for the regression coefficient of the time trend.

Removal rate. The removal rate is between 0 and 1. The inverse of the removal rate, α −1 , corresponds to the average time to removal after infection. We assume

We take ν α 1 = 325.5 and ν α 1 = 35, such that E(α −1 ) = 9.3 with prior 95% credible interval between 8.3 and 10.3 days. The mean infectious period of 9.3 days is chosen based on the findings in He et al. (2020) , who estimated that the infectiousness of COVID-19 starts from around 2.3 days before symptom onset and declines quickly within 7 days after symptom onset.

Diagnosis rate. We place the following standard weakly informative priors on η and σ 2 γ , the regression coefficients and variance term in the diagnosis rate model (Equation 5):

When y t only has an intercept term, we use η ∼ N(0, 1 2 ).

4.1. Posterior Sampling. Let θ = {I U 0 , β, α, µ, σ β , ρ, η, σ 2 γ } denote all model parameters and hyperparameters, where β = (β 0 , β 1 , . . . , β T ), and let B = (B 0 , B 1 , . . . , B T ) be the vector of daily increments in confirmed cases. The joint posterior distribution of θ is given by

where φ(· | µ, σ 2 ) denotes the density function of a normal distribution with mean µ and standard deviation σ 2 , and π * (θ) represents the prior density of θ. Recall that

We use a Markov chain Monte Carlo (MCMC) algorithm (see, e.g., Liu, 2008) We therefore use parallel tempering (PT) to improve the convergence and mixing of the Markov chains (Geyer, 1991) . Consider J parallel Markov chains with a target distribution of

for the j-th chain, where ∆ j is the temperature. The temperatures {∆ 1 , ∆ 2 , . . . , ∆ J } are decreasing with ∆ J = 1. Thus the target distribution of the J-th chain is the original posterior π(θ | B, I D 0 ). At each MCMC iteration, we first independently update all J chains based on Gibbs transition probabilities. Then, for j = 1, 2, . . . , J− 1, we propose a swap between θ j and θ j+1 and accept the proposal with probability

The draws from the J-th chain are kept. A chain with a higher temperature can more freely explore the posterior space, and the swap proposal allows interchange of states between adjacent chains. Therefore, the PT scheme helps the Markov chain avoid getting stuck at local optima. In Appendix B, we demonstrate the advantage of the PT scheme with an example.

In the simulation studies and real data analyses, we run J = 10 parallel Markov chains with a temperature of ∆ j = 1.5 10−j for the j-th chain. We run MCMC simulation for 50,000 iterations, discard the first 20,000 draws as initial burn-in, and keep one sample every 30 iterations. This leaves us a total of 1,000 posterior samples.

Inference. In addition to the estimation of epidemiological parameters, one may be interested in the prediction of a future observation, which can be achieved by sampling from its posterior predictive distribution. As an example, let B * = (B T +1 , . . . , B T +T * ) denote the vector of daily confirmed cases for future days t = T + 1, . . . , T + T * . The posterior predictive distribution of B * is given by

Sampling from (8) involves computing π(β * | B, I D 0 ) = π(β * |β) · π(β | B, I D 0 ) dβ forβ * = (β T +1 , . . . ,β T +T * ). We havẽ

the (i, j)-th entry being C(T + i, j − 1), and C * * is a T * × T * matrix with the (i, j)-th entry being C(T + i, T + j). This is based on a GP prediction rule (Rasmussen and Williams, 2006) .

We assess the performance of the proposed method in estimating the epidemiological parameters by applying it to simulated epidemic time series. Consider a closed population of size N = 20, 000, 000. We assume the initial condition on day 0 is I D 0 = 100, I U 0 = 800, R 0 = 0, and S 0 = N − I D 0 − I U 0 . We set the removal rate α = 9.3 −1 . For the transmission rate, we consider the following three scenarios: Recall that R 0 (t) = β t /α. In all the scenarios, R 0 (t) → 0 + as t → ∞. For scenario 2, R 0 (t) is non-monotonic, and for scenario 3, R 0 (t) is discontinuous. Next, we generate γ t ∼ N [logit(0.2), 0.25 2 ] and γ t = 1 − exp(− exp(γ t )). Finally, for each scenario, we generate a hypothetical epidemic process for 80 days according to Equation (4) with B t = γ t (1 − α)I U t . We keep B = (B 0 , . . . , B T ) and I D 0 as our observations (T = 79). The simulated datasets, shown in Figure 2 (upper panel) , are similar to a real COVID-19 dataset (e.g., Figure 4) .

We fit the proposed model to the simulated datasets using the PTMCMC algorithm. Figure 2 ( 

To illustrate the practical application of the proposed method, we carry out data analysis based on daily counts of confirmed COVID-19 cases reported by JHU CSSE.

This is the B t in our model. We limit our analysis to six U.S. states (Washington, New York, California, Florida, Texas, and Illinois) to keep the paper in reasonable length. The reader can carry out independent analysis for other states, countries or regions using the R package BaySIR. The populations of these states are obtained from U.S. Census Bureau.

6.1. Estimation of the Effective Reproduction Number. Figure 3 shows the estimated R e (t) for the six states. The start dates of statewide stay-at-home orders and state reopening plans are also displayed in the figure for reference (data source: Mervosh et al., 2020 and Washington Post Staff, 2020) . The estimated initial R e ranges from 2.5 to 4.0. Specifically, R e (0) = 2.8, 3.7, 2.6, 3.3, 2.7 and 3.1 for Washington, New York, California, Florida, Texas and Illinois, respectively. During the early stage of the outbreak, the R e generally has a decreasing trend. We suspect that the decline in R e may be associated with the implementation of mitigation policies (e.g., statewide stay-at-home orders, shown in Figure 3 ) and the increase of public awareness. Starting from April, the R e for these states is maintained around or below 1, indicating (partial) containment of the disease. However, with the gradual lift of stay-at-home orders and reopening of businesses, we can clearly observe rebounds of R e for some states (e.g., Florida) since May. For all the states, we can observe local fluctuations of R e over time, which may potentially be attributed to some unobserved factors such as social distancing fatigue. Our analysis is preliminary and does not lead to definitive conclusions about whether a specific intervention is effective in controlling disease spread. Due to the issue of (potentially unmeasured) confounding, it is very challenging to draw causal inference about the effectiveness of an intervention.

and may be used as a reference for decision-makers.

6.2. Test of Fit. We carry out the Bayesian χ 2 test (Johnson, 2004) to assess the goodness-of-fit of our model using Illinois data as an example. First, we choose quantiles 0 ≡ a 0 < a 1 < · · · < a G−1 < a G ≡ 1, with p g = a g − a g−1 , g = 1, . . . , G.

As suggested by Johnson (2004) , we use (a 0 , . . . , a 5 ) = (0, 0.2, 0.4, 0.6, 0.8, 1), so ] falls between the a g−1 and a g quantiles of the distribution

Then, under the null hypothesis of a good model fit, the statistic ω should follow a χ 2 -distribution with G − 1 = 4 degrees of freedom. A quantile-quantile plot of the posterior samples of ω against the expected order statistics from a χ 2 4 distribution (Appendix Figure D) shows that ω plausibly comes from a χ 2 4 distribution. In addition, we find the proportion of posterior samples of ω exceeding the 95% quantile of a χ 2 4 distribution to be 0.053. There is no evidence of a lack of fit.

6.3. Forecasts.

Retrospective forecasts. As described in Section 4.2, the proposed method can be used to predict a future observation based on its posterior predictive distribution.

To evaluate the forecasting performance of the proposed model, we conduct withinsample forecasts using Illinois as an example. Specifically, we split the observations B into a training set B tr and a testing set B te , where B tr = (B 0 , B 1 , . . . , B t * ) and

B te = (B t * +1 , B t * +2 , . . . , B T ). We consider three different scenarios, t * ∈ {19, 39, 59}, so that the training set consists of observations for 20, 40 and 60 days, respectively.

We first sample from the posterior distribution of the parameters evaluated on the training set, π(θ | B tr , I D 0 ), and then sample from the posterior predictive distribution of the testing observations, π(B te | B tr , I D 0 ). Figure 4 shows the forecasting results for the three scenarios. The 97.5% percentile of π(B te | B tr , I D 0 ) (i.e., the upper bound of the 95% credible interval) is truncated in the figure for better display, because it becomes huge with exponential growth.

To better understand the forecasting behavior of the proposed model, the predictions of future R e (t)'s are also displayed. Using 20-day training data, the median of π(B te | B tr , I D 0 ) underestimates the actual observations, although the 95% credible interval covers the observed values. In general, prediction of an epidemic process is challenging, especially when the epidemiological parameters vary over time. To see this, notice that there is a rebound of R e (t) around April 21, which cannot be captured by the GP prediction rule with 20-day training data. Since the stay-at-home order is still in effect on April 21, this rebound can neither be captured by policyrelated covariates. To summarize, future predictions are made based on extrapolation of the current trend, and if the trend changes unexpectedly, the predictions will be inaccurate.

With more training data, the prediction accuracy improves, as seen in Figure 4 (b, c). Using 60-day training data, the median of π(B te | B tr , I D 0 ) matches well with the actual observations. Lastly, the short-term predictions (within, say, the next 10 days) are reasonably accurate in all the scenarios.

Prospective forecasts. To make future predictions, we first sample from π(θ | B, I D 0 ) and then sample from π(B * | B, I D 0 ); recall that B * = (B T +1 , . . . , B T +T * ). Figure  5 shows the projected daily confirmed cases and R e (t)'s for Illinois in the next 30 days (i.e., T * = 30). The projections are based on the assumption that the decreasing trend of R e (t) continues. With the lift of the stay-at-home order and the reopening of businesses, it is possible that R e (t) will rebound, thus caution is needed in interpreting the forecasting results.

We developed a Bayesian approach to statistical inference about the transmission dynamics of COVID-19. We proposed to estimate the disease transmission rate using GPR, which captures nonlinear and non-monotonic trends without the need of specific parametric assumptions. A PTMCMC algorithm was used to efficiently sample from the posterior distribution of the epidemiological parameters. Case studies based on the proposed method revealed the overall decreasing trend of R e in six U.S. states (Washington, New York, California, Florida, Texas, and Illinois), which may be associated with the implementation of mitigation policies and the increasing public awareness of the disease. Projections for future case counts can be made based on extrapolation, although caution is needed in interpreting the forecasting results.

Extensions of the proposed compartmental model can be made in a number of ways.

As described in Section 3.1, it is possible to further split the UI and DI compartments and to incorporate an exposed compartment. We may also split the removed Zhang et al. (2020) and Aguilar et al. (2020) . Considering more compartments will make the model more realistic. However, by adding complexity to the current parsimonious model, sampling, estimation and model unidentifiability problems are likely exacerbated (Capaldi et al., 2012; Osthus et al., 2017) . A possible way out could be to utilize more observable information, such as numbers of recoveries and hospitalizations. Nevertheless, not every country/region has these data (or accurate measurements of these data) available, and we chose to model only daily confirmed cases to keep our method general enough and applicable to most countries/regions. (4) is a state-space model motivated by a deterministic SIR model. A future direction is to consider a stochastic epidemic model. For example, a model similar to Lekone and Finkenstädt (2006) may be used,

Compared to Equation (4), this model may better reflect the stochastic nature of the epidemic process. The cost is increased computational complexity.

In our models for the diagnosis rate (Equation 5) and transmission rate (Equation 7), we allow incorporation of covariates. Currently, we only considered an intercept term and a time trend, because empirically we found it hard to identify the effects of other covariates. Due to Ockham's razor (Jefferys and Berger, 1992) , we preferred the simpler model. More efficient ways to incorporate covariates, potentially based on model selection or variable selection techniques, are worth further investigation.

Our data analysis was carried out separately for each country/region. A nature extension is to model multiple countries/regions jointly using a hierarchical model to achieve borrowing of information, which usually leads to improvements in parameter estimations. We assumed that the population in each country/region is closed, ignoring immigration and emigration. Arguably, a more realistic model should take into account spatial spread of the disease, as seen in Li et al. (2020b) . Again, the main drawbacks to these extensions would be increased computation time.

As discussed in Appendix A, the parameters in model (4) are unidentifiable with only daily confirmed cases (B t ) observed, thus parameter estimates are sensitive to prior choices and modeling assumptions. Many existing models for COVID-19 share the same situation, which could partially explain why different studies may lead to substantially different estimates. For example, some consider the infectious period as the time from infection to recovery or death, which is around 20-30 days (Verity et al., 2020) . Under this definition, the estimated effective reproduction numbers would be higher (e.g., Aguilar et al., 2020) . Therefore, when interpreting the results, it is important to recognize their reliance on underlying assumptions.

Lastly, since the proposed model (4) is a state-space model, it is of interest to further explore online and sequential algorithms for posterior sampling, such as sequential Monte Carlo (Doucet et al., 2001; Dukic et al., 2012) . In that way, when data at more time points become available, one can update the posterior in an efficient way rather than re-fitting the model to the complete data.

With only daily confirmed cases observed, the parameters in model (4) are unidentifiable. To see this, consider the following two epidemic processes (indexed by j = 1

and 2),

for t = 1, . . . , T . The observation is the daily increment in confirmed cases, B j,t = γ j,t (1 − α j )I U j,t . These two processes give rise to identical observations B 1,t and B 2,t for all t, if

and

for t = 1, . . . , T . In other words, different sets of parameters can lead to exactly the same observed data. Even if we restrict that (S 1,0 , I U 1,0 , I D 1,0 , R 1,0 ) = (S 2,0 , I U 2,0 , I D 2,0 , R 2,0 ) (same initial conditions), γ 1,t ≡ γ 1 , and γ 2,t ≡ γ 2 (constant diagnosis rate), for any α 1 = α 2 we can still solve Equations (9) and (10) and get distinct {γ 1 , β 1,t } and {γ 2 , β 2,t } that lead to the same observed data.

A specific example is given below. Consider a population size of N = 20, 000, 000.

Suppose there are two epidemic processes with the same initial conditions, I U 1,0 = I U 2,0 = 800, I D 1,0 = I D 2,0 = 100, R 1,0 = R 2,0 = 0, and S 1,0 = S 2,0 = N − 900. Suppose further α 1 = 0.3, α 2 = 0.05, γ 1,t ≡ γ 1 = 0.2, and γ 2,t ≡ γ 2 = γ 1 · (1 − α 1 )/(1 − α 2 ) = 0.147. Then, the parameters β 1,t and β 2,t can be chosen ( Figure A(a) ) such that {B 1,t } and {B 2,t } are identical ( Figure A(b) ). The resulting effective reproduction numbers for the two epidemic processes, R j e (t) = (β j,t S j,t )/(α j N ), are shown in Figure  A (c) and are quite different. This example highlights that the parameters in (4) are unidentifiable in the absence of strong prior information. (a) β 1,t and β 2,t (b) B 1,t and B 2,t (identical) (c) R 1 e (t) and R 2 e (t) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 

Investigating the impact of asymptomatic carriers on COVID-19 transmission

An introduction to stochastic epidemic models

Infectious Diseases of Humans: Dynamics and Control

Stochastic Epidemic Models and Their Statistical Analysis

Statistical studies of infectious disease incidence

Compartmental Models in Epidemiology

Parameter estimation and uncertainty quantication for an epidemic model

Scenario analysis of non-pharmaceutical interventions on global COVID-19 transmissions

An interactive web-based dashboard to track COVID-19 in real time

An introduction to sequential

Monte Carlo methods

Tracking epidemics with Google flu trends data and a state-space SEIR model

Estimating the number of infections and the impact of non-pharmaceutical interventions on COVID-19

Evaluating the Accuracy of Sampling-based Approaches to the Calculation of Posterior Moments

Markov chain Monte Carlo maximum likelihood

Estimating parameters in stochastic compartmental models using Markov chain methods

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

INFEKTA: A general agent-based model for transmission of infectious diseases: Studying the COVID-19 propagation in

Better strategies for containing COVID-19 epidemics-a study of 25 countries via an extended varying coefficient SEIR model

Temporal dynamics in viral shedding and transmissibility of COVID-19

The mathematics of infectious diseases

IHME COVID-19 health service utilization forecasting team and Murray

Forecasting COVID-19 impact on hospital bed-days, ICU-days, ventilator-days and deaths by US state in the next 4 months

Ockham's razor and Bayesian analysis

A Bayesian χ2 test for goodness-of-fit. The Annals of Statistics

A contribution to the mathematical theory of epidemics

COVID-19 and postinfection immunity: Limited evidence, many remaining questions

Bayesian nonparametrics for stochastic epidemic models

Statistical inference in a stochastic epidemic SEIR model with control intervention: Ebola as a case study

Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

Monte Carlo Strategies in Scientific Computing

Antibody responses to SARS-CoV-2 in patients with COVID-19

See which states and cities have told residents to stay at home. The New York Times

Dynamic Bayesian influenza forecasting in the United States with hierarchical discrepancy (with discussion)

Forecasting seasonal influenza with a state-space SIR model

Bayesian inference for partially observed stochastic epidemics

Association of public health interventions with the epidemiology of the COVID-19 outbreak in Wuhan

Gaussian Process for Machine Learning

Gaussian processes for time-series modelling

Transmission of 2019-nCoV infection from an asymptomatic contact in Germany

Tracking reproductivity of COVID-19 epidemic in China with varying coefficient SIR model

Estimates of the severity of coronavirus disease 2019: a model-based analysis. The Lancet Infectious Diseases

Spatiotemporal dynamics, nowcasting and forecasting of COVID-19 in the United States

An epidemiological forecast model and software assessing interventions on COVID-19 epidemic in China

Where states are reopening after the U.S. shutdown. The Washington Post

The SIR model and the foundations of public health

Virological assessment of hospitalized patients with COVID-2019

Projections for first-wave COVID-19 deaths across the US using social-distancing measures derived from mobile phones

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: A modelling study

Bayesian non-parametric inference for stochastic epidemic models using Gaussian processes

Prediction of the COVID-19 outbreak based on a realistic stochastic model

To demonstrate the advantage of the PT scheme, we show in Figure B the Markov chains for I U 0 and η generated using or not using PT based on a simulated dataset. We evaluate the convergence of the chains using Geweke's diagnostic (Geweke, 1991) .Under the null hypothesis of chain convergence, Geweke's z-score should follow a standard normal distribution. The z-score indicates lack of convergence for the chains generated without PT. Appendix C. Simulation Studies: Sensitivity AnalysisWe carry out sensitivity analyses to explore how the choice of the link function (Equation 6 ) and priors can affect the performance of the proposed method. We consider the following four settings:(Set. 1) Replacing the default logit link for γ t by the probit link;(Set. 2) Replacing the default logit link for γ t by the complementary log-log link;(Set. 3) Replacing the default prior on α −1 by α −1 ∼ Ga(46.5, 5)1(α −1 ≥ 1). This leads to a larger prior variance for α −1 compared to the default. Recall that the default prior is α −1 ∼ Ga(325.5, 35)1(α −1 ≥ 1);(Set. 4) Replacing the default prior on α −1 by α −1 ∼ Ga(700, 35)1(α −1 ≥ 1). This leads to a different prior mean for α −1 compared to the default. We fit our model to the Simulation Scenario 1 dataset. Figure Figure C(a, b) ). Also, increasing the prior variance for α −1 does not lead to much change in the estimates ( Figure C(c) ). Lastly, altering the prior mean for α −1 can lead to substantially different estimates ( Figure   C(d) ). This is due to parameter unidentifiability issues (Appendix A). Multiple solutions may explain the observed data equally well, thus the solutions that are more consistent with the prior would be preferred. Under Setting 4, the prior for α −1 is centered around 20, while the true α −1 = 9.3. As a result, the parameter estimates deviate from the simulation truth.Appendix D. Case Studies: Test of FitWe carry out the Bayesian χ 2 test (Johnson, 2004) to assess the goodness-of-fit of our model using Illinois data as an example. Under the null hypothesis of a good model fit, the statistic ω should follow a χ 2 4 distribution. Figure D shows a quantilequantile plot of posterior samples of ω against expected order statistics from a χ 2 4 distribution. There is no evidence that ω deviates from a χ 2 4 distribution.