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

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 $t$, a fixed population can be split into three compartments: $S(t)$, $I(t)$, and $R(t)$, 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

$$\begin{split}&\frac{dS}{dt}=-\beta SI/N+\phi R,\\ &\frac{dI}{dt}=\beta SI/N-\gamma I,\\ &\frac{dR}{dt}=\gamma I-\phi R,\end{split}$$

(1)

where $\beta$ 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, $\gamma$ 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 $\phi$ 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

$$\frac{dS}{dt}+\frac{dI}{dt}+\frac{dR}{dt}=0.$$

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

$$S(t)+I(t)+R(t)=N,$$

(2)

where $N$ is a fixed total population. In cases with discrete time $t=1,\ldots,T$ (in units of days), we have

$$\begin{split}S(t+1)=S(t)-\beta S(t)I(t)/N+\phi R(t),\\ I(t+1)=I(t)+\beta S(t)I(t)/N-\gamma I(t),\\ R(t+1)=R(t)+\gamma I(t)-\phi R(t),\end{split}$$

(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 $t=1,\ldots,T$, the data models are

$$Z_{R}(t)|P_{R}(t)\sim\text{Possion}(N\times P_{R}(t)),$$

(4)

and

$$Z_{I}(t)|P_{I}(t)\sim\text{Possion}(N\times P_{I}(t)),$$

(5)

where $Z_{R}(t)$ and $Z_{I}(t)$ are the observed number of “recovereds” (includes deaths) and infectious individuals at time $t$, respectively; $N$ is known total number of population and $Z_{S}(t)=N-Z_{I}(t)-Z_{R}(t)$; $P_{R}(t)$ and $P_{I}(t)$ are underlying true rates of recovered and infectious individuals. Thus, our observed data are $\{(Z_{R}(t),Z_{I}(t))\}:t=1,2\ldots,T$. Based on the relationship between the number of “recovereds”, infectious, and suspects, we have

$$P_{R}(t)+P_{I}(t)+P_{S}(t)=1,$$

(6)

where $P_{S}(t)$ underlying rate of susceptible individuals.

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

$$\begin{split}P_{R}(t+1)&=P_{R}(t)+\gamma P_{I}(t)-\phi P_{R}(t),\\ P_{I}(t+1)&=P_{I}(t)+\beta P_{S}(t)P_{I}(t)-\gamma P_{I}(t),\\ P_{S}(t+1)&=P_{S}(t)-\beta P_{S}(t)P_{I}(t)+\phi P_{R}(t).\end{split}$$

(7)

In order to model the hidden uncertainties in SIRS model, we define following transformation of $P_{R}(t)$, $P_{I}(t)$ and $P_{S}(t)$ based on (6)

$$\begin{split}W_{S}(t)&\equiv\log(\frac{P_{S}(t)}{P_{R}(t)}),\\ W_{I}(t)&\equiv\log(\frac{P_{I}(t)}{P_{R}(t)}).\\ \end{split}$$

(8)

The time-varying process of $W_{R}(t)$ and $W_{I}(t)$ is defined as

$$\begin{split}W_{S}(t+1)&=\mu_{S}(t)+\epsilon_{S}(t+1),\\ W_{I}(t+1)&=\mu_{I}(t)+\epsilon_{I}(t+1),\end{split}$$

(9)

where $\epsilon_{S}(t)\sim N(0,\sigma^{2}_{S})$ and $\epsilon_{I}(t)\sim N(0,\sigma^{2}_{I})$ for $t=1,2,\ldots,T$. Based on (6) and (7), we have

$$\begin{split}\mu_{S}(t)=&W_{S}(t)\\ &+\log(1+\frac{\phi}{\exp(W_{S}(t))}-\frac{\beta\exp(W_{I}(t))}{1+\exp(W_{S}(t))+\exp(W_{I}(t))})\\ &+\log(\frac{1}{1+\gamma\exp(W_{I}(t))-\phi})\end{split}$$

(10)

and

$$\begin{split}\mu_{I}(t)=&W_{I}(t)\\ &+\log(1-\gamma+\frac{\beta\exp(W_{S}(t))}{1+\exp(W_{S}(t))+\exp(W_{I}(t))})\\ &+\log(\frac{1}{1+\gamma\exp(W_{I}(t))-\phi})\end{split}$$

(11)

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

$$\begin{split}Z_{R}(t)|W_{S}(t),W_{I}(t)\sim\text{Possion}\left(N\times\frac{1}{1+\exp(W_{S}(t))+\exp(W_{I}(t))}\right),\\ Z_{I}(t)|W_{S}(t),W_{I}(t)\sim\text{Possion}\left(N\times\frac{\exp(W_{I}(t))}{1+\exp(W_{S}(t))+\exp(W_{I}(t))}\right).\end{split}$$

(12)

For the simplicity, we refer the model from (9) to (12) as $\{(Z_{R}(t),Z_{I}(t),N),t=1,2\ldots,T\}\sim\text{SIRS}(\beta,\gamma,\phi,\sigma^{2}_{S},\sigma^{2}_{I})$. Based on the transmission rate and recover rate, the basic reproduction number, $R_{0}$, can be calculated by

$$R_{0}=\frac{\beta}{\gamma}.$$

(13)

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 $n$ different regions. In this section, we focus on the the transmission rate $\beta$ for different regions. Recovery rate and loss of immunity rate can be parameterized in the same way.

Let $z_{1},\ldots,z_{n}\in\{1,\ldots,k\}$ denote clustering labels of $n$ regions and $\beta_{1},\ldots,\beta_{n}$ denote the corresponding parameters in SIRS model for $n$ regions. Our goal is to cluster $\beta_{1},\ldots,\beta_{n}$ into $k$ clusters with distinct values $\beta^{*}_{1},\ldots,\beta^{*}_{k}$, which is usually unknown in practice. A popular solution for unknown $k$ is to introduce the Dirichlet process mixture prior models (Antoniak, 1974) as following:

$\displaystyle\beta_{i}\sim G,\quad G\sim DP(\alpha G_{0}),$

(14)

where $G_{0}$ is a base measure and $\alpha$ is a concentration parameter. If a set of values of $\beta_{1},\ldots,\beta_{n}$ are drawn from $G$, a conditional prior can be obtained by integration (Blackwell et al., 1973):

$\displaystyle p(\beta_{n+1}\mid\beta_{1},\ldots,\beta_{n})=\dfrac{1}{n+\alpha}\sum_{j=1}^{n}\delta_{\beta_{j}}(\beta_{n+1})+\dfrac{\alpha}{n+\alpha}G_{0}(\beta_{n+1}).$

(15)

Here, $\delta_{\beta_{j}}(\beta_{\ell})=I(\beta_{\ell}=\beta_{j})$ is a point mass at $\beta_{j}$. We can obtain the following equivalent models by introducing cluster membership $z_{j}$’s and letting the unknown number of clusters $k$ go to infinity (Neal, 2000).

$$\begin{split}z_{i}\mid\bm{\pi}&\sim\text{Discrete}(\pi_{1},\ldots,\pi_{k}),\\ \beta^{*}_{c}&\sim G_{0}\\ \bm{\pi}&\sim\text{Dirichlet}(\alpha/k,\ldots,\alpha/k)\end{split}$$

(16)

where $\bm{\pi}=(\pi_{1},\ldots,\pi_{k})$. In addition, the distribution of $z_{i}$ can be marginally given by a stick-breaking representation (Sethuraman, 1994) of Dirichlet process (DP) as

$$\begin{split}z_{i}&\sim\sum_{h=1}^{\infty}\pi_{h}\delta_{h},\\ \pi_{h}&=\nu_{h}\prod_{\ell\leq h}(1-\nu_{\ell}),\\ \nu_{h}&\sim\text{Beta}(1,\alpha),\end{split}$$

(17)

where $\delta_{h}$ is the Dirac function with mass at $h$.

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

$\displaystyle k\sim p(\cdot),\quad(\pi_{1},\ldots,\pi_{k})\mid k\sim\mbox{Dirichlet}(\eta,\ldots,\eta),\quad z_{i}\mid k,\bm{\pi}\sim\sum_{h=1}^{k}\pi_{h}\delta_{h},\quad i=1,\ldots,n,$

(18)

where $p(\cdot)$ is a proper probability mass function (p.m.f.) on $\{1,2,\ldots,\}$.

Like the stick-breaking representation in (17) of Dirichlet process, the MFM also has a similar construction. If we choose $k-1\sim\mbox{Poisson}(\lambda)$ and $\eta=1$ in (18), the mixture weights $\pi_{1},\cdots,\pi_{k}$ is constructed as follows:

1. 1.

   Generate $\eta_{1},\eta_{2},\cdots\overset{\text{iid}}{\sim}\text{Exp}(\lambda)$,

2. 2.

   $k=\min\{j:\sum_{i=1}^{j}\eta_{i}\geq 1\}$,

3. 3.

   $\pi_{i}=\eta_{i}$, for $i=1,\cdots,k-1$,

4. 4.

   $\pi_{k}=1-\sum_{i}^{k-1}\pi_{i}$.

For ease of exposition, we refer the stick-breaking representation of MFM above as $\text{MFM}(\lambda)$ with default choice of $p(\cdot)$ being $\mbox{Poisson}(\lambda)$ and $\eta=1$.

In order to allow for simultaneously heterogeneity learning of three parameters in SIRS model, the MFM prior is introduced for parameters $\beta$, $\gamma$ and $\phi$ in the SIRS model. Our observed data are $\{(Z_{R}(t,\bm{s}_{i}),Z_{I}(t,\bm{s}_{i}),N_{i}):t=1,2,\ldots,T,i=1,2,\ldots,n\}$, where $t$ denotes each discrete time point and $i$ denotes each state. The hierarchical SIRS model with MFM is given as

$$\begin{split}&\{(Z_{R}(t,\bm{s}_{i}),Z_{I}(t,\bm{s}_{i}),N_{i}),t=1,2\ldots,T\}\sim\text{SIRS}(\beta_{z^{\beta}_{i}},\gamma_{z^{\gamma}_{i}},\phi_{z^{\phi}_{i}},\sigma^{2}_{S,i},\sigma^{2}_{I,i}),i=1,2,\ldots,n\\ &z^{\beta}_{i}\sim\text{MFM}(\lambda_{\beta}),i=1,2,\ldots,n,\\ &z^{\gamma}_{i}\sim\text{MFM}(\lambda_{\gamma}),i=1,2,\ldots,n,\\ &z^{\phi}_{i}\sim\text{MFM}(\lambda_{\phi}),i=1,2,\ldots,n,\\ &\beta_{z^{\beta}_{i}}\sim G_{\beta},\\ &\gamma_{z^{\gamma}_{i}}\sim G_{\gamma},\\ &\phi_{z^{\phi}_{i}}\sim G_{\phi},\\ &\sigma^{2}_{S,i},\sigma^{2}_{I,i}\sim\text{IG}(0.01,0.01),i=1,2,\ldots,n,\end{split}$$

(19)

where $z^{\beta}$, $z^{\gamma}$, and $z^{\phi}$ denote the cluster assignments of parameter $\beta$,$\gamma$, and $\phi$, respectively. $G_{\beta}$, $G_{\gamma}$, and $G_{\phi}$ is the base distribution for parameter $\beta$,$\gamma$, and $\phi$, respectively. The choices of $G_{\beta}$, $G_{\gamma}$, and $G_{\phi}$ will be discussed in Section 3.4.

For the hierarchical SIRS model with MFM introduced in Section 3.3, the set of parameters is denoted as $\Theta=\{\beta_{z^{\beta}_{i}},\gamma_{z^{\gamma}_{i}},\phi_{z^{\phi}_{i}},\sigma^{2}_{S,i},\sigma^{2}_{I,i},\lambda_{\beta},\lambda_{\gamma},\lambda_{\phi}:i=1,2\ldots,n\}$. 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 $G_{\beta}$, $G_{\gamma}$ and $G_{\phi}$, respectively:

$$\begin{split}&\beta_{z_{i}^{\beta}}\sim\text{Uniform}(0,1),\\ &\gamma_{z^{\gamma}_{i}}\sim\text{Uniform}(0,1),\\ &\phi_{z^{\phi}_{i}}\sim\text{Uniform}(0,1).\end{split}$$

(20)

For the hyperparameters for three MFM processes, we assign gamma prior $\text{Gamma}(1,1)$ on $\lambda_{\beta},\lambda_{\gamma},\lambda_{\phi}$. With the joint prior distributions $\pi(\Theta)$, the posterior distribution of these parameters based on the data $D=\{(Z_{R}(t,\bm{s}_{i}),Z_{I}(t,\bm{s}_{i}),N_{i}):t=1,2,\ldots,T,i=1,2,\ldots,n\}$ is given as

$$\pi(\Theta|(Z_{R}(t,\bm{s}_{i}),Z_{I}(t,\bm{s}_{i}),N_{i}):t=1,2,\ldots,T,i=1,2,\ldots,n)\propto L(D|\Theta)\times\pi(\Theta),$$

(21)

where $L(D|\Theta)$ is the full data likelihood given from the model (9) to (12). The analytical form of the posterior distribution of $\pi(\Theta|(Z_{R}(t,\bm{s}_{i}),Z_{I}(t,\bm{s}_{i}),N_{i}):t=1,2,\ldots,T,i=1,2,\ldots,n)$ 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.

After obtaining posterior samples, an important task is do inference on posterior samples. Using posterior mean or posterior median for grouping label $\bm{z}$ 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, $B^{(1)},\ldots,B^{(M)}$, where $B^{(t)}$ for the $t$-th post-burn-in MCMC iteration is defined as:

$\displaystyle B^{(t)}=[B^{(t)}(i,j)]_{i,j\in\{1:n\}}=1(z_{i}^{(t)}=z_{j}^{(t)})_{n\times n},~{}~{}t=1,\ldots,M,$

(22)

Here $1(\cdot)$ denotes the indicator function, which means $B^{(t)}(i,j)=1$ indicates observations $i$ and $j$ are in the same cluster in the $t$-th posterior sample post burn-in. After obtaining the membership matrices of the posterior samples, a Euclidean mean for membership matrices is calculated by:

$$\overline{B}=\frac{1}{M}\sum_{t=1}^{M}B^{(t)}.$$

Based on $\overline{B}$ and $B^{(1)},\ldots,B^{(M)}$, we find the iteration with the least squares distance to $\overline{B}$ as

$\displaystyle C_{L}=\text{argmin}_{t\in(1:M)}\sum_{i=1}^{n}\sum_{j=1}^{n}\{B(i,j)^{(t)}-\overline{B}(i,j)\}^{2}.$

(23)

The estimated parameters, together with the cluster assignments $\bm{z}$, are then extracted from the $C_{L}$-th post burn-in iteration.

In order to mimic the real dataset we analyze, we choose $n=51$ and the population for each location is assigned as the real data population for $51$ states. The time length $T$ 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 $\beta_{1}=0.06,\beta_{2}=0.6$, $\phi_{1}=0.06,\phi_{2}=0.6,$, and $\gamma_{1}=0.06,\gamma_{2}=0.6$. We randomly assign the labels to $51$ locations and fix them over 100 replicates. The grouping labels for three parameters is given in Figure 2.

For each replicates, we have $15,000$ iterations MCMC samples and have first $5,000$ iterations burn-in in order to obtain samples from every 5th iteration of the last $10,000$ 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 $\beta$ as an example:

	$$\displaystyle\text{MAB}=\frac{1}{100}\sum_{r=1}^{100}\bigg{\{}\frac{1}{n}\sum_{i=1}^{n}\widehat{\beta}^{r}(\bm{s}_{i})-\beta(\bm{s}_{i})\bigg{\}},$$
	$$\displaystyle\text{MSD}=\sqrt{\frac{1}{100}\sum_{r=1}^{100}\bigg{\{}\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{\beta}^{r}(\bm{s}_{i})-\overline{\widehat{\beta}}(\bm{s}_{i})\right)^{2}}\bigg{\}},$$

where $\overline{\widehat{\beta}}(\bm{s}_{i})$ is the mean of the posterior estimate over 100 replicates.

For clustering estimation evaluation, the estimated number of clusters $\widehat{K}$ 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.

The parameter estimation performance and clustering performance results are shown in Table 1 and Table 2.

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.

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 $50,000$ MCMC iterations and burnin the first $20,000$ iterations in order to obtain samples from every 10th iteration of the last $30,000$ iterations. The group labels are obtained by Dalh’s method in Section 3.5.

For $\beta$, one group is identified. $\beta=0.079$ with $95\%$ Highest Probability Density (HPD) interval (Chen and Shao, 1999) $(0.058,0.098)$. For $\gamma$, three groups are identified with $\gamma_{1}=0.0054$ with HPD interval $(0.0021,0.0207)$, $\gamma_{2}=0.0419$ with HPD interval $(0.0022,0.0609)$ and $\gamma_{3}=0.0164$ with HPD interval$(0.0035,0.0241)$. 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 $\gamma$s are shown in Figure 3.

For $\phi$, one group is identified. $\phi=0.0015$ with HPD interval $(1.181\times 10^{-7},0.0047)$.

With the estimated values of $\beta$ and $\gamma$, the basic reproduction number, $R_{0}$, is calculated among different states. The values of $R_{0}$ among different states are shown in Figure 4.

The second time period we analyze is from May 1st, 2020. Other settings are same with previous analysis.

For $\beta$, one group is identified. $\beta=0.0042$ with $95\%$ Highest Probability Density (HPD) interval (Chen and Shao, 1999) $(3.056\times 10^{-8},0.1083)$. Compared with previous 30-day data, in this time period, the transmission rate decrese a lot. For $\gamma$, two groups are identified with $\gamma_{1}=0.0381$ with HPD interval $(0.0048,0.3713)$ and $\gamma_{2}=0.0007$ with HPD interval $(0.0003,0.0013)$. Two states, Oregon and Vermont, are identified in group 1. Other states are identified in group 2. For $\phi$, one groups is identified. $\phi=0.0006$ with HPD interval $(2.747\times 10^{-7},0.0026)$.

With the estimated values of $\beta$ and $\gamma$, the basic reproduction number, $R_{0}$, is calculated among different states. There are two different groups for the basic reproduction number. One group include Oregon and Vermont with $R_{0}=0.1102$. The other group includes other 49 states with $R_{0}=5.4619$. Comparing to the 30 days period from April 1st, we can see a decrease for $R_{0}$ in general.