Heterogeneity Learning for SIRS model: an Application to the COVID-19

by   Guanyu Hu, et al.

We propose a Bayesian Heterogeneity Learning approach for Susceptible-Infected-Removal-Susceptible (SIRS) model that allows underlying clustering patterns for transmission rate, recovery rate, and loss of immunity rate for the latest coronavirus (COVID-19) among different regions. Our proposed method provides simultaneously inference on parameter estimation and clustering information which contains both number of clusters and cluster configurations. Specifically, our key idea is to formulates the SIRS model into a hierarchical form and assign the Mixture of Finite mixtures priors for heterogeneity learning. The properties of the proposed models are examined and a Markov chain Monte Carlo sampling algorithm is used to sample from the posterior distribution. Extensive simulation studies are carried out to examine empirical performance of the proposed methods. We further apply the proposed methodology to analyze the state level COVID-19 data in U.S.



There are no comments yet.


page 14

page 17


Time Fused Coefficient SIR Model with Application to COVID-19 Epidemic in the United States

In this paper, we propose a Susceptible-Infected-Removal (SIR) model wit...

Bayesian Heterogeneity Pursuit Regression Models for Spatially Dependent Data

Most existing spatial clustering literatures discussed the cluster algor...

Graphical Assistant Grouped Network Autoregression Model: a Bayesian Nonparametric Recourse

Vector autoregression model is ubiquitous in classical time series data ...

Analysis of professional basketball field goal attempts via a Bayesian matrix clustering approach

We propose a Bayesian nonparametric matrix clustering approach to analyz...

Bayesian Group Learning for Shot Selection of Professional Basketball Players

In this paper, we develop a group learning approach to analyze the under...

Bayesian Spatial Homogeneity Pursuit of Functional Data: an Application to the U.S. Income Distribution

An income distribution describes how an entity's total wealth is distrib...

Spatial homogeneity learning for spatially correlated functional data with application to COVID-19 Growth rate curves

We study the spatial heterogeneity effect on regional COVID-19 pandemic ...
This week in AI

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

1 Introduction

The Coronavirus Disease 2019 (COVID-19) 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 COVID-19 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 COVID-19 has been a significant public-health concern and people pay more attention on precise and timely estimates and predictions of COVID-19. The Susceptible-Infectious-Recovered (SIR) model and its variation approaches, such as Susceptible-Infected-Removal-Susceptible (SIRS) (Kermack and McKendrick, 1932, 1933) and Susceptible-Exposed-Infected-Removal (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 COVID-19 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 time-varying model under SIR and its variations is also fully discussed for COVID-19 (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 COVID-19. They lack discussions on heterogeneity pattern of COVID-19 among different regions.

The aim of this paper is to propose a new hierarchical SIRS model for detecting heterogeneity pattern of COVID-19 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 pre-specify 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 over-clustering 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 two-fold. 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 COVID-19 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 COVID-19 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 COVID-19 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 COVID-19, and 607,279 people are recovered from COVID-19. The fatality rate of COVID-19 is .

Figure 1: Exploratory Analysis of COVID-19 on June 10th

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



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 that

Thus, 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


where is a fixed total population. In cases with discrete time (in units of days), we have


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




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


where underlying rate of susceptible individuals.

Similar to (Zhuang and Cressie, 2014), we have following hidden processes:


In order to model the hidden uncertainties in SIRS model, we define following transformation of , and based on (6)


The time-varying process of and is defined as


where and for . Based on (6) and (7), we have




Based on the transformation in (8), we have our data in (4) and (5) as


For the simplicity, we refer the model from (9) to (12) as . Based on the transmission rate and recover rate, the basic reproduction number, , can be calculated by


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:


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):


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).


where . In addition, the distribution of can be marginally given by a stick-breaking representation (Sethuraman, 1994) of Dirichlet process (DP) as


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):


where  is a proper probability mass function (p.m.f.) on .

Like the stick-breaking 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:

  1. Generate ,

  2. ,

  3. , for ,

  4. .

For ease of exposition, we refer the stick-breaking 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


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:


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


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 post-burn-in MCMC iteration is defined as:


Here denotes the indicator function, which means indicates observations and  are in the same cluster in the -th posterior sample post burn-in. 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


The estimated parameters, together with the cluster assignments , are then extracted from the -th post burn-in 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.

Figure 2: Grouping Labels for , , and

For each replicates, we have iterations MCMC samples and have first iterations burn-in 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 R-package 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
Table 1: Estimation Performance for Simulation Study
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
Table 2: Grouping Performance for Simulation Study

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 30-Day Analysis from April 1st

We analyze 30-Day 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 corona-virus https://www.whitehouse.gov/presidential-actions/ on March 11th, 2020. From the April 1st, we can assume that there are very limited imported cases from outside U.S.. We analyze 30-day 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 .

Figure 3: Group Labels for s of April 1st data

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.

Figure 4: for 51 States from April 1st

5.2 30-Day 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 30-day 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 COVID-19 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 COVID-19 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.


  • 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 time-dependent sir model for covid-19 with undetectable infected persons. arXiv preprint arXiv:2003.00122.
  • Dahl (2006) Dahl, D. B. (2006). Model-based 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. Anderson-Bergman, 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 state-space 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 covid-19 spread in south korea using the sir model with time-dependent parameters and deep learning.

  • 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 2019-ncov: 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 covid-19 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 2019-ncov 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 2019-ncov 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.