1 Introduction
The Coronavirus Disease 2019 (COVID19) has created a profound public health emergency around world. It has become an epidemic with more than 5,000,000 confirmed infections worldwide as on May 21 2020. The spreading speed of COVID19 which is caused by a new coronavirus is faster than severe acute respiratory syndrome (SARS) and Middle East respiratory syndrome (MERS). Recently, the risk of COVID19 has been a significant publichealth concern and people pay more attention on precise and timely estimates and predictions of COVID19. The SusceptibleInfectiousRecovered (SIR) model and its variation approaches, such as SusceptibleInfectedRemovalSusceptible (SIRS) (Kermack and McKendrick, 1932, 1933) and SusceptibleExposedInfectedRemoval (SEIR) model (Hethcote, 2000), have been widely discussed to study the dynamical evolution of an infectious disease in a certain region. There are rich literatures producing early results on COVID19 based on SIR model and its variations (Wu et al., 2020; Read et al., 2020; Tang et al., 2020). From statistician’s perspectives, building a timevarying model under SIR and its variations is also fully discussed for COVID19 (Chen et al., 2020; Sun et al., 2020; Jo et al., 2020). In most existing literature, people focus more on dynamic regimes of the SIR models for COVID19. They lack discussions on heterogeneity pattern of COVID19 among different regions.
The aim of this paper is to propose a new hierarchical SIRS model for detecting heterogeneity pattern of COVID19 among different regions under a Bayesian framework. Bayesian nonparametric methods such as Dirichlet process (DP) offer choices to do simultaneously inference on parameters’ estimation and parameters’ heterogeneity information which contains the number of clusters and clustering configurations. Compared with existing approaches such as finite mixtures models, Bayesian nonparametric approach does not need to prespecify the number of clusters, which provides probabilistic framework for simultaneous inference of the number of clusters and the clustering labels. Miller and Harrison (2013) points out that the estimation of the number of clusters under Dirichlet process mixture (DPM) model is inconsistent which will produce extremely small clusters. One remedy for overclustering problem under DPM is mixture of finite mixtures model (MFM) (Miller and Harrison, 2018). The clustering properties of MFM are fully discussed in Miller and Harrison (2018); Geng et al. (2019) and it has been widely applied in different areas such as regional economics (Hu et al., 2020), environmental science (Geng et al., 2019), and social science (Geng et al., 2019). Thus, the key idea of this paper is to assign MFM priors on different parameters of the SIRS model to capture the heterogeneity of parameters among different regions. The contribution of this paper are twofold. First, we formulate a Bayesian heterogeneity learning model for SIRS under MFM. To our best knowledge, this is the first time when MFM is applied into epidemiology models such as SIRS. Our proposed Bayesian approach successfully captures the heterogeneity of three different parameters under the SIRS model among different regions while also considering uncertainty in estimation of the number of clusters. Several interesting findings based on our proposed method are discovered for COVID19 data in US.
This paper is organized as follows. Section 2 presents the motivating data we analyze. We discuss our proposed Bayesian hierarchical model for heterogeneity learning under SIRS model framework in Section 3. The performance of our proposed method is illustrated via simulation studies in Section 4. Section 5 is devoted to the analysis of state level COVID19 data in U.S. A brief discussion is presented in Section 6.
2 Motivating Data
Our motivating data comes from the COVID tracking project https://covidtracking.com. State Level COVID19 Data are recorded for the 50 states plus Washington, DC. For simplicity, we refer to them as “51 states” in the rest of this paper. Up to June 10th, 2020, United States totally confirmed 2,043,031 cases. 114,533 people died because of COVID19, and 607,279 people are recovered from COVID19. The fatality rate of COVID19 is .
Figure 1 shows state level confirmed numbers, death numbers, incident rate, and mortality rate. We can see that New York state has the highest confirmed number, death numbers, and incident rate; Connecticut has the highest mortality rate among 51 states; Montana has the least confirmed number; Alaska has the least death number; the incident rate of Hawaii is lowest among 51 states; and Texas has the lowest mortality rate.
3 Method
3.1 SIRS Model
Compartment epidemic models are widely used to study infectious disease which spreads through human populations across a large region. SIR model (Kermack and McKendrick, 1927) has been universally discussed for analyzing the dynamical evolution of an infectious disease in a large population. SIR model is extended to SIRS model for imperfect immunity situation (Kermack and McKendrick, 1932, 1933). For a given time , a fixed population can be split into three compartments: , , and
, which denotes the number of susceptibles, the number of infectious, and the number of “recovereds” (which includes deaths), respectively. The dynamical process of SIRS model can be written as following nonlinear ordinary differential equations of three given compartments
(1) 
where
denotes the average rate of contact per unit time multiplied by the probability of disease transmission per contact between a susceptible and an infectious subject,
denotes the rate of “recovery” per unit time, which is the rate at which infectious individuals are removed from being infectious due to recovery (or death), and denotes the rate of loss of immunity of recovered individuals per unit time, which is the rate at which recovered individuals become susceptible again (Anderson et al., 1992; Zhuang and Cressie, 2014). By adding the equations in (1), we notice thatThus, the model postulates a fixed total population without entry and exits of demographic type. For example, there are no births or deaths caused by other than the disease we study in a certain time. The SIRS model assume the sum of all three compartments together is constant within a short period of time such that
(2) 
where is a fixed total population. In cases with discrete time (in units of days), we have
(3) 
with the same constraints as (2).
Based on the models in (3) and (2) and assmuptions in (Dukic et al., 2012)
, the data model of SIRS assumes conditional independent Poisson distributions evolving at discrete time points. For a given time
, the data models are(4) 
and
(5) 
where and are the observed number of “recovereds” (includes deaths) and infectious individuals at time , respectively; is known total number of population and ; and are underlying true rates of recovered and infectious individuals. Thus, our observed data are . Based on the relationship between the number of “recovereds”, infectious, and suspects, we have
(6) 
where underlying rate of susceptible individuals.
Similar to (Zhuang and Cressie, 2014), we have following hidden processes:
(7) 
3.2 Heterogeneity Learning
In section 2, our motivating data is at state level in US and we are interested in whether there are heterogeneity patterns on the parameters of interest among different states. As an assumption, we believe that different states might have different parameters, however, some states will share similar pattern in transmission rate, recovery rate, or loss of immunity rate. Next, we introduce nonparametric Bayesian methods for heterogeneity learning of SIRS parameters over different regions. In this section, we focus on the the transmission rate for different regions. Recovery rate and loss of immunity rate can be parameterized in the same way.
Let denote clustering labels of regions and denote the corresponding parameters in SIRS model for regions. Our goal is to cluster into clusters with distinct values , which is usually unknown in practice. A popular solution for unknown is to introduce the Dirichlet process mixture prior models (Antoniak, 1974) as following:
(14) 
where is a base measure and is a concentration parameter. If a set of values of are drawn from , a conditional prior can be obtained by integration (Blackwell et al., 1973):
(15) 
Here, is a point mass at . We can obtain the following equivalent models by introducing cluster membership ’s and letting the unknown number of clusters go to infinity (Neal, 2000).
(16) 
where . In addition, the distribution of can be marginally given by a stickbreaking representation (Sethuraman, 1994) of Dirichlet process (DP) as
(17) 
where is the Dirac function with mass at .
However, Miller and Harrison (2018) proved that the DP mixture model produces extraneous clusters in the posterior leading to inconsistent estimation of the number of clusters even when the sample size grows to infinity. A modification of DP mixture model called Mixture of finite mixtures (MFM) model is proposed to circumvent this issue (Miller and Harrison, 2018):
(18) 
where is a proper probability mass function (p.m.f.) on .
Like the stickbreaking representation in (17) of Dirichlet process, the MFM also has a similar construction. If we choose and in (18), the mixture weights is constructed as follows:

Generate ,

,

, for ,

.
For ease of exposition, we refer the stickbreaking representation of MFM above as with default choice of being and .
3.3 Hierarchical Model
In order to allow for simultaneously heterogeneity learning of three parameters in SIRS model, the MFM prior is introduced for parameters , and in the SIRS model. Our observed data are , where denotes each discrete time point and denotes each state. The hierarchical SIRS model with MFM is given as
(19) 
where , , and denote the cluster assignments of parameter ,, and , respectively. , , and is the base distribution for parameter ,, and , respectively. The choices of , , and will be discussed in Section 3.4.
3.4 Prior and Posterior
For the hierarchical SIRS model with MFM introduced in Section 3.3, the set of parameters is denoted as . To complete the model, we now specify the joint prior distribution for the parameters. Based on the natural constraints generated by (3), we have following distribution for bases distribution , and , respectively:
(20) 
For the hyperparameters for three MFM processes, we assign gamma prior
on . With the joint prior distributions , the posterior distribution of these parameters based on the data is given as(21) 
where is the full data likelihood given from the model (9) to (12). The analytical form of the posterior distribution of is unavailable. Therefore, we carry out the posterior inference using the MCMC sampling algorithm to sample from the posterior distribution and then obtain the posterior estimates of the unknown parameters. Computation is facilitated by the nimble package (de Valpine et al., 2017) in R which generates C++ code for faster computation.
3.5 Group Inference via MCMC Samples
After obtaining posterior samples, an important task is do inference on posterior samples. Using posterior mean or posterior median for grouping label is not suitable. Instead, inference on the clustering configurations is obtained employing the modal clustering method of (Dahl, 2006). The inference is based on the membership matrices of posterior samples, , where for the th postburnin MCMC iteration is defined as:
(22) 
Here denotes the indicator function, which means indicates observations and are in the same cluster in the th posterior sample post burnin. After obtaining the membership matrices of the posterior samples, a Euclidean mean for membership matrices is calculated by:
Based on and , we find the iteration with the least squares distance to as
(23) 
The estimated parameters, together with the cluster assignments , are then extracted from the th post burnin iteration.
4 Simulation
In this section, we investigate the performance of the hierarchical SIRS model with MFM from a variety of measures.
4.1 Simulation Settings and Evaluation Metrics
In order to mimic the real dataset we analyze, we choose and the population for each location is assigned as the real data population for states. The time length equals 30 for all the simulation replicates. The total number of replicates in our simulation study is 100. For each parameter, we have two different groups and we set the true values of the parameters , , and . We randomly assign the labels to locations and fix them over 100 replicates. The grouping labels for three parameters is given in Figure 2.
For each replicates, we have iterations MCMC samples and have first iterations burnin in order to obtain samples from every 5th iteration of the last iterations.
The performance of the posterior estimates of parameters were evaluated by the mean bias (MB) and the mean standard deviation (MSD) in the following ways, take
as an example:where is the mean of the posterior estimate over 100 replicates.
For clustering estimation evaluation, the estimated number of clusters for each replicate is summarized from the MCMC iteration picked by Dahl’s method. Rand Index (RI; Rand, 1971) is applied to evaluate cluster configuration. The RI is calculated by Rpackage fossil (Vavrek, 2011). A higher value of the RI represents higher accuracy of clustering The average RI (MRI) was calculated as the mean of RIs over the 100 replicates.
4.2 Simulation Results
The parameter estimation performance and clustering performance results are shown in Table 1 and Table 2.
Parameter  MB  MSD 

0.008  0.021  
0.072  0.152  
0.007  0.017  
0.068  0.151  
0.012  0.023  
0.069  0.149 
Parameter  MRI  S.D of RI  S.D. of  

0.854  0.058  2.12  0.33  
0.857  0.057  2.33  0.55  
0.847  0.059  2.31  0.54  
From the results shown in Table 1, we see that the MABs and MSDs of the parameters are both within a reasonable range. In general, performance of posterior estimates achieve a good target.
And we see that our proposed methods successfully recover the number of groups and grouping labels within a reasonable range for all three parameters from Table 2. Average rand index for all parameters around 0.85 indicate our proposed method truly recovers the group labels for all three parameters. The mean of the estimated number of groups is close to true number of groups over 100 replicates.
5 Real Data Analysis
5.1 30Day Analysis from April 1st
We analyze 30Day data from April 1st, 2020. The reason why we analyze this time period data is that U.S. Government announced the suspension of entry as immigrants and nonimmigrants of certain additional persons who pose a risk of transmitting coronavirus https://www.whitehouse.gov/presidentialactions/ on March 11th, 2020. From the April 1st, we can assume that there are very limited imported cases from outside U.S.. We analyze 30day data based on the model in (19) and use the priors discussed in Section 3.4. We run MCMC iterations and burnin the first iterations in order to obtain samples from every 10th iteration of the last iterations. The group labels are obtained by Dalh’s method in Section 3.5.
For , one group is identified. with Highest Probability Density (HPD) interval (Chen and Shao, 1999) . For , three groups are identified with with HPD interval , with HPD interval and with HPD interval. Thirty three states including Alabama, Arizona, California, Colorado, Connecticut, District of Columbia, Florida, Georgia, Idaho, Illinois, Indiana, Kansas, Louisiana, Massachusetts, Michigan, Mississippi, Missouri, Nebraska, Nevada, New Jersey, North Carolina, Ohio, Oregon, Pennsylvania, Rhode Island, South Carolina, Texas, Utah, Vermont, Virginia, Washington, West Virginia, Wisconsin, belong to group one. Eleven states including Arkansas, Hawaii, Iowa, Maine, Minnesota, New Hampshire, North Dakota, Oklahoma, South Dakota, Tennessee, Wyoming , belong to group 2. And seven states including Alaska, Delaware, Kentucky, Maryland, Montana, New Mexico, New York, belong to group 3. The grouping labels for s are shown in Figure 3.
For , one group is identified. with HPD interval .
With the estimated values of and , the basic reproduction number, , is calculated among different states. The values of among different states are shown in Figure 4.
5.2 30Day Analysis from May 1st
The second time period we analyze is from May 1st, 2020. Other settings are same with previous analysis.
For , one group is identified. with Highest Probability Density (HPD) interval (Chen and Shao, 1999) . Compared with previous 30day data, in this time period, the transmission rate decrese a lot. For , two groups are identified with with HPD interval and with HPD interval . Two states, Oregon and Vermont, are identified in group 1. Other states are identified in group 2. For , one groups is identified. with HPD interval .
With the estimated values of and , the basic reproduction number, , is calculated among different states. There are two different groups for the basic reproduction number. One group include Oregon and Vermont with . The other group includes other 49 states with . Comparing to the 30 days period from April 1st, we can see a decrease for in general.
6 Discussion
In this paper, we develop a nonparametric Bayesian heterogeneity learning method for SIRS model based on Mixture of Finite Mixtures model. This statistical framework was motivated by the heterogeneity of COVID19 pattern among different regions. Our simulation results indicate that the proposed method can recover the heterogeneity pattern of parameters among different regions. Illustrated by the analysis of COVID19 data in U.S., our proposed methods reveal the heterogeneity pattern among different states.
In addition, three topics beyond the scope of this paper are worth further investigation. First, we can add spatially dependent structure (Hu et al., 2020; Zhao et al., 2020) on the heterogeneity of different states. Second, our model assumes parameters are constant over a certain time period. Building heterogeneity learning model with time varying coefficients is an interesting future work. Finally, proposing a measurement error correction for SIRS devotes another interesting future work.
References
 Anderson et al. (1992) Anderson, R. M., B. Anderson, and R. M. May (1992). Infectious diseases of humans: dynamics and control. Oxford university press.
 Antoniak (1974) Antoniak, C. E. (1974). Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. The Annals of Statistics 2(6), 1152–1174.
 Blackwell et al. (1973) Blackwell, D., J. B. MacQueen, et al. (1973). Ferguson distributions via Pólya urn schemes. The Annals of Statistics 1(2), 353–355.
 Chen and Shao (1999) Chen, M.H. and Q.M. Shao (1999). Monte carlo estimation of bayesian credible and hpd intervals. Journal of Computational and Graphical Statistics 8(1), 69–92.
 Chen et al. (2020) Chen, Y.C., P.E. Lu, C.S. Chang, and T.H. Liu (2020). A timedependent sir model for covid19 with undetectable infected persons. arXiv preprint arXiv:2003.00122.
 Dahl (2006) Dahl, D. B. (2006). Modelbased clustering for expression data via a dirichlet process mixture model. Bayesian inference for gene expression and proteomics 4, 201–218.
 de Valpine et al. (2017) de Valpine, P., D. Turek, C. J. Paciorek, C. AndersonBergman, D. T. Lang, and R. Bodik (2017). Programming with models: writing statistical algorithms for general model structures with nimble. Journal of Computational and Graphical Statistics 26(2), 403–413.
 Dukic et al. (2012) Dukic, V., H. F. Lopes, and N. G. Polson (2012). Tracking epidemics with google flu trends data and a statespace seir model. Journal of the American Statistical Association 107(500), 1410–1426.
 Geng et al. (2019) Geng, J., A. Bhattacharya, and D. Pati (2019). Probabilistic community detection with unknown number of communities. Journal of the American Statistical Association 114(526), 893–905.
 Geng et al. (2019) Geng, J., W. Shi, and G. Hu (2019). Bayesian nonparametric nonhomogeneous poisson process with applications to usgs earthquake data. arXiv preprint arXiv:1907.03186.
 Hethcote (2000) Hethcote, H. W. (2000). The mathematics of infectious diseases. SIAM review 42(4), 599–653.
 Hu et al. (2020) Hu, G., J. Geng, Y. Xue, and H. Sang (2020). Bayesian spatial homogeneity pursuit of functional data: an application to the us income distribution. arXiv preprint arXiv:2002.06663.

Jo
et al. (2020)
Jo, H., H. Son, S. Y. Jung, and H. J. Hwang (2020).
Analysis of covid19 spread in south korea using the sir model with timedependent parameters and deep learning.
medRxiv.  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.
 Kermack and McKendrick (1932) Kermack, W. O. and A. G. McKendrick (1932). Contributions to the mathematical theory of epidemics. ii.—the problem of endemicity. Proceedings of the Royal Society of London. Series A, containing papers of a mathematical and physical character 138(834), 55–83.
 Kermack and McKendrick (1933) Kermack, W. O. and A. G. McKendrick (1933). Contributions to the mathematical theory of epidemics. iii.—further studies of the problem of endemicity. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 141(843), 94–122.
 Miller and Harrison (2013) Miller, J. W. and M. T. Harrison (2013). A simple example of Dirichlet process mixture inconsistency for the number of components. In Advances in Neural Information Processing Systems, pp. 199–206.
 Miller and Harrison (2018) Miller, J. W. and M. T. Harrison (2018). Mixture models with a prior on the number of components. Journal of the American Statistical Association 113(521), 340–356.
 Neal (2000) Neal, R. M. (2000). Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics 9(2), 249–265.
 Rand (1971) Rand, W. M. (1971). Objective criteria for the evaluation of clustering methods. Journal of the American Statistical association 66(336), 846–850.
 Read et al. (2020) Read, J. M., J. R. Bridgen, D. A. Cummings, A. Ho, and C. P. Jewell (2020). Novel coronavirus 2019ncov: early estimation of epidemiological parameters and epidemic predictions. MedRxiv.
 Sethuraman (1994) Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica, 639–650.

Sun
et al. (2020)
Sun, H., Y. Qiu, H. Yan, Y. Huang, Y. Zhu, J. Gu, and S. X. Chen (2020).
[discussion paper] tracking reproductivity of covid19 epidemic in
china with varying coefficient sir model.
Journal of Data Science
, 2.  Tang et al. (2020) Tang, B., X. Wang, Q. Li, N. L. Bragazzi, S. Tang, Y. Xiao, and J. Wu (2020). Estimation of the transmission risk of the 2019ncov and its implication for public health interventions. Journal of clinical medicine 9(2), 462.
 Vavrek (2011) Vavrek, M. J. (2011). Fossil: palaeoecological and palaeogeographical analysis tools. Palaeontologia Electronica 14(1), 16.
 Wu et al. (2020) Wu, J. T., K. Leung, and G. M. Leung (2020). Nowcasting and forecasting the potential domestic and international spread of the 2019ncov outbreak originating in wuhan, china: a modelling study. The Lancet 395(10225), 689–697.
 Zhao et al. (2020) Zhao, P., H.C. Yang, D. K. Dey, and G. Hu (2020). Bayesian spatial homogeneity pursuit regression for count value data. arXiv preprint arXiv:2002.06678.
 Zhuang and Cressie (2014) Zhuang, L. and N. Cressie (2014). Bayesian hierarchical statistical sirs models. Statistical Methods & Applications 23(4), 601–646.
Comments
There are no comments yet.