1. Introduction
The outbreak of Coronavirus Disease 2019 (COVID19), caused by Severe Acute Respiratory Syndrome Coronavirus 2 (SARS‑CoV‑2), was declared a pandemic on March 11, 2020 by the World Health Organization. As of July 2, 2020, the number of confirmed COVID19 cases worldwide has exceeded 10.6 million, and the death toll has surpassed 516,000. In order to control the spread of the virus, countries around the world have implemented unprecedented nonpharmaceutical interventions, such as case isolation, closure of schools, stayathome orders, banning of mass gatherings, and local and national lockdowns. At the same time, social distancing and mask wearing by the public also contribute to the containment of COVID19.
Researchers have made substantial efforts to study the transmission dynamics of COVID19, evaluate the effects of government interventions, and forecast infection and death counts. These works include aguilar2020investigating; chen2020scenario; flaxman2020estimating; giordano2020modelling; gomez2020infekta; gu2020better; ihme2020forecasting; li2020early; li2020substantial; pan2020association; sun2020tracking; wang2020spatiotemporal; wang2020epidemiological; woody2020projections; wu2020nowcasting; zhang2020prediction, 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) agentbased modeling. Curve fitting approaches fit a curve to the observed number of confirmed cases or deaths. For example, ihme2020forecasting use a Gaussian error function to model the cumulative death rate at a specific location. Compartmental modeling approaches (e.g., li2020substantial) 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, agentbased modeling approaches (e.g., gomez2020infekta) 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 COVID19, which is critical for characterizing disease spread. We aim to address a few issues related to the COVID19 pandemic. First, we provide estimation of key epidemiological parameters, such as the effective reproduction number of COVID19. 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 COVID19, 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 COVID19 (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 COVID19 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 susceptibleinfectiousrecovered model (kermack1927contribution). We will use data of daily confirmed COVID19 cases reported by the Center for Systems Science and Engineering at Johns Hopkins University (JHU CSSE) (dong2020interactive). We provide an R package BaySIR, available at https://github.com/tianjianzhou/BaySIR, that can be used to conduct independent analysis of COVID19 data or reproduce the results in this paper.
The proposed Bayesian approach attempts to improve COVID19 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 presymptomatic and asymptomatic COVID19 cases (rothe2020transmission; he2020temporal), 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 nonmonotonic relationships without the need of specific parametric assumptions. Third, we develop a paralleltempering 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, crossvalidation and goodnessoffit tests. Such validations provide insights into the modeling of COVID19 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 susceptibleinfectiousrecovered (SIR) compartmental model. In Section 3, we develop a probabilistic statespace model for COVID19 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 COVID19 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.
2. Review of the SusceptibleInfectiousRecovered Model
We start with a review of the susceptibleinfectiousrecovered (SIR) model (kermack1927contribution; weiss2013sir), 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 . Here, “closed” means that does not vary over time. It is a good approximation for a fastspreading and less fatal pandemic like COVID19. The SIR model divides the population into the following three compartments:

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

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

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 (), denote by , and the numbers of individuals in the S, I and R compartments, respectively, and write . We have .
2.1. Deterministic SIR Models
The classical SIR model (kermack1927contribution) describes the flow of people from S to I to R via the following system of differential equations:
(1) 
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, of these contacts are with susceptible individuals per unit time, and as a result, infectious individuals lead to a rate of new infections . The third equation in (1) describes that the infectious individuals leave the infective class at a rate of . 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 and are referred to as the basic reproduction number and effective reproduction number, respectively, where is the initial number of susceptibles at time .
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, is treated as a stochastic process. A commonly used formulation is as follows (gibson1998estimating; o1999bayesian; andersson2000stochastic). Suppose that an infectious individual makes effective contacts with any given individual in the population at times given by a Poisson process of rate , 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
. It can be shown thatis a Markov process with transition probabilities:
(3) 
Here, is a small increment in time.
2.3. Statespace SIR Models
There are, of course, other ways to model the uncertainty of the epidemic process. Probabilistic statespace modeling approaches that build on deterministic models have recently been popular in the statistics literature (dukic2012tracking; osthus2017forecasting; osthus2019dynamic). A statespace 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 osthus2017forecasting has the form
Evolution:  
Observation: 
for . In the evolution model, is the solution to Equation (1) at time with a initial value of at time and parameters and , and is assumed to be centered at
with its variance characterized by
. In other words, measures the derivation of from the solution given by the deterministic model. In the observation model, 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 . The observation is assumed to be centered at with variance characterized by . Statespace 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 susceptibleexposedinfectiousrecovered (SEIR) model includes an additional compartment for exposed individuals who are exposed to the disease but are not yet infectious, and the susceptibleinfectiousrecoveredinfectious (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, anderson1991infectious, hethcote2000mathematics or Brauer2008. For a comprehensive review of stochastic epidemic models, see, for example, becker1999statistical, andersson2000stochastic or allen2008introduction.
3. Proposed Model for COVID19
We now turn to our proposed model for the COVID19 data, which belongs to the statespace model category (Section 2.3). Our approach integrates the discretetime deterministic SIR model (Equation 2) and semiparametric Bayesian inference. To capture some unique features of COVID19, 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 SARSCoV2 have not been tested for the virus thus are not detected or reported as having the infection (li2020substantial). Second, we allow some epidemiological parameters (such as the disease transmission rate ) to be timevarying to reflect the impact of mitigation policies such as stayathome 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 COVID19 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 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:

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

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;

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

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 Sindividuals with a transmission rate of . After being infected, a Sindividual first becomes an UIindividual before being detected as a DIindividual. All the infectious (UI and DI) individuals recover or die with a removal rate of . Those UIindividuals 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 DIindividuals, 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 COVID19 confers immunity to reinfection, although there is only limited evidence for this assumption (long2020antibody; kirkcaldy2020covid).
We define day as the date when the 100th case is confirmed in the country/region under consideration, and index subsequent dates by , where 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 , , and the numbers of individuals belonging to compartments S, UI, DI and R on day , respectively. We have . The transmission rate and diagnosis rate are allowed to vary over time and are hereafter denoted by and , respectively. The number of individuals diagnosed with the disease between day and day is observed and is denoted by . This is our data. We propose modeling the transmission dynamics of COVID19 over time by the following equations:
(4) 
for . Denote by . The epidemic process, , is determined by its initial value , the parameters , and the observations . Rigorously speaking,
should be a vector of nonnegative integers, but for computational convenience, we relax this restriction and only require it to be a vector of nonnegative real numbers. Model (
4) is a simple extension of (2) by adding a component of , the undocumented infections, and by incorporating the observed daily new cases into the equations. Later, we introduce a model for the observation to complete the statespace model.With timevarying disease transmission rates, the basic reproduction number and effective reproduction number are also functions of time. That is, and
Here, is interpreted as the rate of secondary infections generated by each infectious case at time , scaled by the length of the infectious period (). If for , then the number of infectious individuals will monotonically decrease after time , 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 indicates containment of the disease. Due to the important role of in characterizing disease spread, we consider the estimation of as our main interest.
3.2. Model for the Observed Data
Our observations only consist of the daily new confirmed COVID19 cases, . Assume that on day , the UIindividuals who have not been removed are diagnosed with the disease with a diagnosis rate of . Mathematically, this means , where
is between 0 and 1. We consider the logit transformation of
, . Other transformations, such as the probit and complementary loglog 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,(5) 
where is a vector of covariates that are thought to be related to the diagnosis rate. In other words, the sampling model for can be written as
(6) 
In the simulation studies and real data analyses, we use a simple choice of , assuming the mean diagnosis rate is a constant. It is possible to include other covariates in , 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 are the regression coefficients and variance term, respectively, where 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 wolfel2020virological and he2020temporal, the ability of a COVID19 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”.
3.3. Prior Specification
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 COVID19, 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 , i.e., . As a result, the number of DIindividuals on day , , equals to the cumulative number of confirmed cases on that day and is observed. We further assume
where
refers to a gamma distribution with shape and rate parameters
and , respectively. We set and , such that . This choice is based on the findings in li2020substantial that 86% of all infections were undocumented at the beginning of the epidemic in China. Lastly, note that .Transmission rate
The disease transmission rate must be nonnegative. We consider and assume
where refers to a Gaussian process (GP) with mean function and covariance function . The GP (rasmussen2006gaussian) is a very flexible prior model for a stochastic process. It enables one to capture potential nonlinear relationships between and without the need to impose any parametric assumptions. Specifically, for any , the vector
follows a multivariate Gaussian distribution with mean
and covariance matrix with the th entry being . For applications of GP to epidemic modeling, see, for example, xu2016bayesian and kypraios2018bayesian.We specify and as below:
(7) 
Here, 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 , which contains an intercept term and the time. Other covariates, such as indicators for mitigation policies at time , may also be included in . Nevertheless, in practice, we find our GP model with a time trend is sufficient to capture the change of 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 characterizes the amplitude of the difference between and , and the correlation parameter characterizes the correlation between and for any and
. We note that based on our specification of the covariance function, our GP model is equivalent to a firstorder autoregressive model. Indeed, autoregressive models of any orders are discretetime equivalents of GP models with Matérn covariance functions
(roberts2013gaussian).We place the following priors on , and :
such that and . Here, refers to an inverse gamma distribution, and
refers to a beta distribution. The prior choices for
and shrinktoward 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 and , whererepresents 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
li2020early and wu2020nowcasting. 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, , corresponds to the average time to removal after infection. We assume
We take and , such that with prior 95% credible interval between and days. The mean infectious period of 9.3 days is chosen based on the findings in he2020temporal, who estimated that the infectiousness of COVID19 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 , the regression coefficients and variance term in the diagnosis rate model (Equation 5):
When only has an intercept term, we use .
4. Inference
4.1. Posterior Sampling
Let
denote all model parameters and hyperparameters, where
, and let be the vector of daily increments in confirmed cases. The joint posterior distribution of is given bywhere
denotes the density function of a normal distribution with mean
, and represents the prior density of . Recall that .We use a Markov chain Monte Carlo (MCMC) algorithm (see, e.g., liu2008monte), in particular the Gibbs sampler, to simulate from the posterior distribution and implement posterior inference. MetropolisHastings steps are used when the conditional posterior distribution of a parameter is not available in closed form. The regular Gibbs sampler is not very efficient in our application because of the strong correlations among the model parameters. This issue was also noted by osthus2017forecasting. We therefore use parallel tempering (PT) to improve the convergence and mixing of the Markov chains (geyer1991markov). Consider parallel Markov chains with a target distribution of
for the th chain, where is the temperature. The temperatures are decreasing with . Thus the target distribution of the th chain is the original posterior . At each MCMC iteration, we first independently update all chains based on Gibbs transition probabilities. Then, for , we propose a swap between and and accept the proposal with probability
The draws from the 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 parallel Markov chains with a temperature of for the th chain. We run MCMC simulation for 50,000 iterations, discard the first 20,000 draws as initial burnin, and keep one sample every 30 iterations. This leaves us a total of 1,000 posterior samples.
4.2. Predictive 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 denote the vector of daily confirmed cases for future days . The posterior predictive distribution of is given by
(8) 
Sampling from (8) involves computing for . We have
where , , is a matrix with the th entry being , and is a matrix with the th entry being . This is based on a GP prediction rule (rasmussen2006gaussian).
5. Simulation Studies
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 . We assume the initial condition on day is , , , and . We set the removal rate . For the transmission rate, we consider the following three scenarios:

, where and are chosen such that , and ;

, where and are chosen such that , and ;

, where represents the largest integer that is smaller than .
Recall that . In all the scenarios, as . For scenario 2, is nonmonotonic, and for scenario 3, is discontinuous. Next, we generate and . Finally, for each scenario, we generate a hypothetical epidemic process for 80 days according to Equation (4) with . We keep and as our observations (). The simulated datasets, shown in Figure 2 (upper panel), are similar to a real COVID19 dataset (e.g., Figure 4).
We fit the proposed model to the simulated datasets using the PTMCMC algorithm. Figure 2 (lower panel) shows a comparison of the estimated timevarying effective reproduction numbers with the simulation truth. The simulation truth is nicely recovered, and the 95% credible intervals of ’s always cover the true values.
(a) Scenario 1  (b) Scenario 2  (c) Scenario 3 
We also 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. Details of the sensitivity analyses are reported in Appendix C. In general, our method is robust to different specifications of the link function. The choice of the priors, on the other hand, may have an impact on the parameter estimates, because of parameter unidentifiability issues (see Appendix A).
6. Case Studies
To illustrate the practical application of the proposed method, we carry out data analysis based on daily counts of confirmed COVID19 cases reported by JHU CSSE. This is the 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 for the six states. The start dates of statewide stayathome orders and state reopening plans are also displayed in the figure for reference (data source: mervosh2020see and washington2020where). The estimated initial ranges from 2.5 to 4.0. Specifically, , , , , and for Washington, New York, California, Florida, Texas and Illinois, respectively. During the early stage of the outbreak, the generally has a decreasing trend. We suspect that the decline in may be associated with the implementation of mitigation policies (e.g., statewide stayathome orders, shown in Figure 3) and the increase of public awareness. Starting from April, the for these states is maintained around or below 1, indicating (partial) containment of the disease. However, with the gradual lift of stayathome orders and reopening of businesses, we can clearly observe rebounds of for some states (e.g., Florida) since May. For all the states, we can observe local fluctuations of 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. Nevertheless, our analysis can shed light on the transmission dynamics of COVID19 and may be used as a reference for decisionmakers.
6.2. Test of Fit
We carry out the Bayesian test (johnson2004bayesian)
to assess the goodnessoffit of our model using Illinois data as an example. First, we choose quantiles
, with , . As suggested by johnson2004bayesian, we use , so and . Next, let be a posterior sample of the model parameters , and let denote the number of observations (i.e., ’s) such that falls between the and quantiles of the distribution . LetThen, under the null hypothesis of a good model fit, the statistic
should follow a distribution with degrees of freedom. A quantilequantile plot of the posterior samples of against the expected order statistics from a distribution (Appendix Figure D) shows that plausibly comes from a distribution. In addition, we find the proportion of posterior samples of exceeding the 95% quantile of a 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 into a training set and a testing set , where and . We consider three different scenarios, , 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, , and then sample from the posterior predictive distribution of the testing observations, .
(a) 20day training data  (b) 40day training data  (c) 60day training data 
Figure 4 shows the forecasting results for the three scenarios. The 97.5% percentile of (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 ’s are also displayed. Using 20day training data, the median of 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 around April 21, which cannot be captured by the GP prediction rule with 20day training data. Since the stayathome 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 60day training data, the median of matches well with the actual observations. Lastly, the shortterm predictions (within, say, the next 10 days) are reasonably accurate in all the scenarios.
Prospective forecasts
To make future predictions, we first sample from and then sample from ; recall that . Figure 5 shows the projected daily confirmed cases and ’s for Illinois in the next 30 days (i.e., ). The projections are based on the assumption that the decreasing trend of continues. With the lift of the stayathome order and the reopening of businesses, it is possible that will rebound, thus caution is needed in interpreting the forecasting results.
(a) Projected daily confirmed cases  (b) Projected 
7. Discussion
We developed a Bayesian approach to statistical inference about the transmission dynamics of COVID19. We proposed to estimate the disease transmission rate using GPR, which captures nonlinear and nonmonotonic 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 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 compartment into recovered and deceased compartments. See, for example, giordano2020modelling, zhang2020prediction and aguilar2020investigating. 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 (capaldi2012parameter; osthus2017forecasting). 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.
The proposed model as in Equation (4) is a statespace model motivated by a deterministic SIR model. A future direction is to consider a stochastic epidemic model. For example, a model similar to lekone2006statistical may be used,
where
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 (jefferys1992ockham), 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 li2020substantial. 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 () observed, thus parameter estimates are sensitive to prior choices and modeling assumptions. Many existing models for COVID19 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 (verity2020estimates). Under this definition, the estimated effective reproduction numbers would be higher (e.g., aguilar2020investigating). Therefore, when interpreting the results, it is important to recognize their reliance on underlying assumptions.
Lastly, since the proposed model (4) is a statespace model, it is of interest to further explore online and sequential algorithms for posterior sampling, such as sequential Monte Carlo (doucet2001introduction; dukic2012tracking). In that way, when data at more time points become available, one can update the posterior in an efficient way rather than refitting the model to the complete data.
References
Appendix A Parameter Identifiability
With only daily confirmed cases observed, the parameters in model (4) are unidentifiable. To see this, consider the following two epidemic processes (indexed by and ),
for . The observation is the daily increment in confirmed cases, . These two processes give rise to identical observations and for all , if
(9) 
and
(10) 
for . In other words, different sets of parameters can lead to exactly the same observed data. Even if we restrict that (same initial conditions), , and (constant diagnosis rate), for any we can still solve Equations (9) and (10) and get distinct and that lead to the same observed data.
A specific example is given below. Consider a population size of . Suppose there are two epidemic processes with the same initial conditions, , , , and . Suppose further , , , and . Then, the parameters and can be chosen (Figure A(a)) such that and are identical (Figure A(b)). The resulting effective reproduction numbers for the two epidemic processes, , 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) and  (b) and (identical)  (c) and 
Appendix B Posterior Sampling: Parallel Tempering
To demonstrate the advantage of the PT scheme, we show in Figure B the Markov chains for and generated using or not using PT based on a simulated dataset. We evaluate the convergence of the chains using Geweke’s diagnostic (geweke1991evaluating). Under the null hypothesis of chain convergence, Geweke’s score should follow a standard normal distribution. The score indicates lack of convergence for the chains generated without PT.
(a) using PT,  (b) not using PT, 
(c) using PT,  (d) not using PT, 
Appendix C Simulation Studies: Sensitivity Analysis
We 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:

Replacing the default logit link for by the probit link;

Replacing the default logit link for by the complementary loglog link;

Replacing the default prior on by . This leads to a larger prior variance for compared to the default. Recall that the default prior is ;

Replacing the default prior on by . This leads to a different prior mean for compared to the default.
We fit our model to the Simulation Scenario 1 dataset. Figure C shows the estimated timevarying effective reproduction numbers under the four settings. The estimates are robust to the choice of the link function (Figure C(a, b)). Also, increasing the prior variance for does not lead to much change in the estimates (Figure C(c)). Lastly, altering the prior mean for 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 is centered around 20, while the true . As a result, the parameter estimates deviate from the simulation truth.
(a) Probit link  (b) Complementary loglog link 
(c) Larger prior variance for  (d) Different prior mean for 
Appendix D Case Studies: Test of Fit
We carry out the Bayesian test (johnson2004bayesian) to assess the goodnessoffit of our model using Illinois data as an example. Under the null hypothesis of a good model fit, the statistic should follow a distribution. Figure D shows a quantilequantile plot of posterior samples of against expected order statistics from a distribution. There is no evidence that deviates from a distribution.
Comments
There are no comments yet.