# Parameter Estimation on a Stochastic SIR Model with Media Coverage

Authors: Changguo Li; Yongzhen Pei; Meixia Zhu; Yue Deng
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

Media coverage reduces the transmission rate from infective to susceptible individuals and is reflected by suitable nonlinear functions in mathematical modeling of the disease. We here focus on estimating the parameters in the transmission rate based on a stochastic SIR epidemic model with media coverage. In order to reduce the computational load, the Newton-Raphson algorithm and Markov Chain Monte Carlo (MCMC) technique are incorporated with maximum likelihood estimation. Simulations validate our estimation results and the necessity of a model with media coverage when modeling the contagious diseases.


---

## Body

## 1. Introduction

The spread of a contagious disease can trigger responses of people so as to minimize the effect of the disease onto them, and that prevents themselves from contracting the disease [1, 2]. Many scholars have explored the influence of media awareness from mathematical models. Cui et al. [3] modeled the transmission rate involving the media effect by the function βexp(-αI(t)), in which the parameters α and β represent the media impact (hence we call these parameters media parameters) and the transmission rate before media alert, respectively, and I(t) denotes the number of infected individuals at time t. Liu et al. [4] proposed another transmission rate β-β1I(t)/(m+I(t)) to capture the impact of media on disease spread. Here the media parameter β1(<β) refers to the reduced maximum value of the transmission rate when the number of infected individuals I(t) approaches infinite and m reflects the reactive velocity of media coverage to the epidemic disease. Literatures concerning with media coverage mainly focus on the dynamic behaviors of epidemic models with media coverage by way of surveys, qualitative approaches, or numerical simulations [2–7].Comparatively the problem of parameter estimation concerning transmission rate involving media effect has not been fully discussed. One main challenge in these inferential tasks is attributed to the burden of computational load of minimizing/maximizing corresponding object functions, including likelihood functions or the squared differences [8–11]. In this paper, we consider a stochastic Susceptible-Infected-Removed (SIR) model [12] with media coverage. The advantage of the model [12] is that an explicit likelihood functions can be formulated, which enables us to estimate parameters of the model. Newton-Raphson algorithms and Bayesian inference techniques are adopted to alleviate the computational load.The paper is organized as follows. The stochastic SIR model with media coverage effect is introduced in Section2. Then the likelihood function of the model is derived, from which unknown parameters can be estimated via Newton-Raphson algorithm or Markov Chain Monte Carlo (MCMC) technique in Section 3. Finally, some simulations are included to help illustrate the necessity to take media coverage into account when modeling the transmission dynamics of infectious disease in Section 4.


## 2. Model Formulation

Katriel [12] proposed a stochastic discrete time model with infection age. Suppose that the population of size N is partitioned into three classes, the susceptible, infectious, and recovered, with the numbers denoted by S(t), I(t), and R(t) at day t, respectively. Furthermore, newly infected number i(t) is introduced on day t. Suppose the infected individuals can remain “infective” for d days. And the period from the moment an individual became infected to the present is named “age-of-infection” of the individuals which is denoted by τ. Hence the number of the infective with infection age τ(1≤τ≤d) on day t is i(t-τ).Faced with growing number of the infective, people will reduce their chances of contacting with others for fear of being infected. For this reason, we adopt the nonlinear contact rateβexp(-λi^(t-1)). Here parameter β is constant and i^(t-1)=i(t-1)/N is the intensity of infection occurring on day t-1, and parameter λ reflects the extent to which media coverage affects society. If λ=0, the contact rate is constant β which has been used by many classical models. As λ increases the alertness of the public to the disease, the public will be more aware of the diseases which reduces the contact rate as such.In order to deduce the likelihood, we divide(0,t^) into discrete “days” and each day is divided into m small intervals of length △t=1/m. In view of the impacts of media, a susceptible encounters an individual with infection age τ in each small interval with probability (1)βexp⁡-λi^t-1it-τN-11mand gets infected with probability (2)βexp⁡-λi^t-1N-11m∑τ=1dPτit-τ. Here we assume, when a susceptible contacts with an infective individual whose infection age is τ(1≤τ≤d), the susceptible becomes infective with probability Pτ. Therefore the possibility that this susceptible escapes the infection in each day is (3)1-βexp⁡-λi^t-1N-11m∑τ=1dPτit-τm⟶exp⁡-βexp⁡-λi^t-1N-1∑τ=1dPτit-τ as m→∞. As a result, the probability that at least one susceptible becomes infected during day t is given by (4)pt=1-e-βexp⁡-λi^t-1/N-1∑τ=1dPτit-τ. Given S(t-1), the number of newly infected individuals i(t) at day t is binomially distributed with parameters (S(t-1) and p(t)). That is, (5)it~BSt-1,pt.Since the number of infected individuals each day is equal to the reduction of susceptible of that day, we readily have (6)St=St-1-it.Based on (5) and (6) not only can we construct the likelihood so as to estimate the parameters, they are also crucial for simulations afterward. In order to realize this process, we can produce a binomial random number by applying (5) and (6) iteratively for t≥1 when starting values i(t)(-(d-1)≤t≤0) and S(0) are known.When the population is large, a deterministic model is obtained from the stochastic one in the sense of “thermodynamic limit”. Leti^(t)=i(t)/N and S^(t)=S(t)/N be the proportions of the individuals which turn into infective and susceptible by day t, respectively. Given information about the infection before day t, the conditional expected size of infected individuals on day t is (7)St-11-e-βexp⁡-λi^t-1/N-1∑τ=1dPτit-τ. Hence we derive the infected fraction at day t by the formula (8)i^t=1NSt-11-e-βexp⁡-λi^t-1/N-1∑τ=1dPτit-τ=S~t-11-e-βexp⁡-λi^t-1N/N-1∑τ=1dPτi^t-τ. Then a deterministic model is yielded as N→∞(9)i^t=S~t-11-e-βexp⁡-λi^t-1∑τ=1dPτi^t-τ,(10)S~t=S~t-1-i^t.In order to check the relationship between the stochastic model and its deterministic counterpart, we simulate the stochastic models with different values of population sizes. Together with solutions of model (9) and (10), the result (Figure 1) shows the model (9) is indeed the limit of its stochastic model as N increases.Figure 1
Three realizations are plotted for three values ofN (103,104,105 from left to right) with respect to the stochastic and its deterministic model. The rest parameters are β=2,λ=12,d=5,Pτ=1/5(1≤τ≤5), including the starting values, i(t)=0(-4≤τ≤-1),i(0)=15,S(-5)=N.


## 3. Estimation of Parameters

In the following, the parameters in the model (5) will be estimated. Assume that the data i(t)(1≤t≤t^) during t^ days are available. We will first deal with the case that S(0) is known and next consider its estimation when it is unknown.


### 3.1. Estimatingβ,λ

Before the formal estimation, we consider the approximation of the original model (5) and (6). For relatively large N(11)βexp⁡-λi^t-1N-1∑τ=1dPτit-τ≪1,and thus (12)1-e-βexp⁡-λi^t-1/N-1∑τ=1dPτit-τ≈βexp⁡-λi^t-1N-1∑τ=1dPτit-τ.Due to this approximation, the random variable i(t) is approximately Poisson distributed. That is, (13)it~Poissonβexp⁡-λi^t-1St-1N-1∑τ=1dPτit-τ.Notice that the likelihood is the probability of the following event: staring from given numbersi(1),…,i(d), the data i(d+1),…,i(T) are observed. Hence the likelihood of model of (5) and (6) is given by (14)Lβ,λ=Pid+1,…,iT∣i1,…,it=∏t=d+1TPit∣i1,…,it-1=∏t=d+1Tβexp⁡-λi^t-1St-1/N-1∑τ=1dPτit-τitit!×exp⁡-βexp⁡-λi^t-1St-1N-1∑τ=1dPτit-τ=Ae-β∑t=d+1Texp⁡-λi^t-1Gtβne-λhT, where (15)A=∏t=d+1T1it!St-1N-1∑τ=1dPτit-τit,Gt=St-1N-1∑τ=1dPτit-τ,n=∑t=d+1Tit,hT=1N∑t=d+1Tit-1it. The log-likelihood is thus (16)LLβ,λ=log⁡A-β∑t=d+1Texp⁡-λi^t-1Gt+nlog⁡β-λhT. Differentiating with respect to β and λ, we obtain (17)dLLβ,λdβ=-∑t=d+1Texp-λi^t-1Gt+n1β,dLLβ,λdλ=β∑t=d+1Texp⁡-λi^t-1i^t-1Gt-hT. Therefore the maximum likelihood estimators for β~ and λ~ are the solutions to equations: (18)β^=n∑t=d+1Te-λ^i^t-1Gt,(19)0=n∑t=d+1Te-λ^i^t-1Gt∑t=d+1Te-λ^i^t-1i^t-1Gt-hT.By the Newton-Raphson algorithm [13], it is easy to evaluate λ^ from (19). Substituting λ^ into (18), we obtain the estimator β^.To check the validity of the estimators, 100 realizations of (5) and (6) are generated with β=2.6,λ=15. In each case, two estimators β^ and λ^ are evaluated and plotted in Figure 2, which shows that β^ and λ^ give suitable estimates of the original value of β and λ, respectively.Figure 2
Histograms of 100 estimators ofβ=2.6 and λ=15 with d=5,Pτ=1/5(1≤τ≤5) and N=104. The starting values are i(t)=10(1≤t≤5),S(0)=N.


### 3.2. EstimatingS0, β, and λ

In real applications, the initial susceptible numberS0 is usually unknown either. Hence we also need to estimate S0 from data it, t=1,2,…,T. Noticing that St-1=S0-rt with rt=∑s=1t-1is, we have(20)Pit∣i1,…,it-1=exp⁡-S0-rtβe-λi^gt×S0-rtβe-λi^gtit1it!,where g(t)=1/N-1∑τ=1d‍Pτit-τ. Accordingly the likelihood function of the model is given by (21)Lβ,λ,S0∣x=∏t=d+1TPit∣i1,…,it-1=∏t=d+1Texp⁡-S0-rtβe-λi^gt×S0-rtβe-λi^gtit1it!∝βN×exp⁡-β∑t=d+1TS0-rtexp-λi^t-1gt×∏t=d+1TS0-rtit×e-λhT.Owing to the difficulty in maximizingL(β,λ,S0∣x), we pose it as a Bayes parameter inference problem. The prior for β and λ is chosen to be noninformative, i.e., π(β)=1/β and π(λ)=1/λ. In addition S(0) is uniformly distributed on (r(t),N). By assumption of independence, the prior distribution of (β,λ,S0) is defined as (22)πβ,λ,S0∝1β1λIrt,NS0. In this way we get the posterior (23)πβ,λ,S0∣x∝πβ,λ,S0Lβ,λ,S0∣x. By drawing samples from π(β,λ,S(0)∣x) using a MCMC scheme, the Bayesian estimates λ^B, β^B, and S^B(0) of β,λ, and S(0) are obtained by the sample means of these simulating parameters, respectively.To examine the quality of these Bayesian estimators, a population of sizeN=20,000 with S0=12,000 initially susceptible is considered. The parameters are set as β=2.6, λ=28, d=5, Pτ=1/5(1≤τ≤5) and i(t)=10 for 1≤t≤5. The histograms of the estimates are presented in Figure 3 which shows that β^B,λ^B, and S^B(0) are suitable estimates of the original value of β,λ, and S0, respectively.Figure 3
Histograms of50,000 Bayesian estimators of β,λ, and S(0). The data are simulated from the model with parameters β=2.6,λ=28,S(0)=12,000, d=5,Pτ=1/5(1≤τ≤5), N=20,000, and i(t)=10 for 1≤t≤5.


## 3.1. Estimatingβ,λ

Before the formal estimation, we consider the approximation of the original model (5) and (6). For relatively large N(11)βexp⁡-λi^t-1N-1∑τ=1dPτit-τ≪1,and thus (12)1-e-βexp⁡-λi^t-1/N-1∑τ=1dPτit-τ≈βexp⁡-λi^t-1N-1∑τ=1dPτit-τ.Due to this approximation, the random variable i(t) is approximately Poisson distributed. That is, (13)it~Poissonβexp⁡-λi^t-1St-1N-1∑τ=1dPτit-τ.Notice that the likelihood is the probability of the following event: staring from given numbersi(1),…,i(d), the data i(d+1),…,i(T) are observed. Hence the likelihood of model of (5) and (6) is given by (14)Lβ,λ=Pid+1,…,iT∣i1,…,it=∏t=d+1TPit∣i1,…,it-1=∏t=d+1Tβexp⁡-λi^t-1St-1/N-1∑τ=1dPτit-τitit!×exp⁡-βexp⁡-λi^t-1St-1N-1∑τ=1dPτit-τ=Ae-β∑t=d+1Texp⁡-λi^t-1Gtβne-λhT, where (15)A=∏t=d+1T1it!St-1N-1∑τ=1dPτit-τit,Gt=St-1N-1∑τ=1dPτit-τ,n=∑t=d+1Tit,hT=1N∑t=d+1Tit-1it. The log-likelihood is thus (16)LLβ,λ=log⁡A-β∑t=d+1Texp⁡-λi^t-1Gt+nlog⁡β-λhT. Differentiating with respect to β and λ, we obtain (17)dLLβ,λdβ=-∑t=d+1Texp-λi^t-1Gt+n1β,dLLβ,λdλ=β∑t=d+1Texp⁡-λi^t-1i^t-1Gt-hT. Therefore the maximum likelihood estimators for β~ and λ~ are the solutions to equations: (18)β^=n∑t=d+1Te-λ^i^t-1Gt,(19)0=n∑t=d+1Te-λ^i^t-1Gt∑t=d+1Te-λ^i^t-1i^t-1Gt-hT.By the Newton-Raphson algorithm [13], it is easy to evaluate λ^ from (19). Substituting λ^ into (18), we obtain the estimator β^.To check the validity of the estimators, 100 realizations of (5) and (6) are generated with β=2.6,λ=15. In each case, two estimators β^ and λ^ are evaluated and plotted in Figure 2, which shows that β^ and λ^ give suitable estimates of the original value of β and λ, respectively.Figure 2
Histograms of 100 estimators ofβ=2.6 and λ=15 with d=5,Pτ=1/5(1≤τ≤5) and N=104. The starting values are i(t)=10(1≤t≤5),S(0)=N.


## 3.2. EstimatingS0, β, and λ

In real applications, the initial susceptible numberS0 is usually unknown either. Hence we also need to estimate S0 from data it, t=1,2,…,T. Noticing that St-1=S0-rt with rt=∑s=1t-1is, we have(20)Pit∣i1,…,it-1=exp⁡-S0-rtβe-λi^gt×S0-rtβe-λi^gtit1it!,where g(t)=1/N-1∑τ=1d‍Pτit-τ. Accordingly the likelihood function of the model is given by (21)Lβ,λ,S0∣x=∏t=d+1TPit∣i1,…,it-1=∏t=d+1Texp⁡-S0-rtβe-λi^gt×S0-rtβe-λi^gtit1it!∝βN×exp⁡-β∑t=d+1TS0-rtexp-λi^t-1gt×∏t=d+1TS0-rtit×e-λhT.Owing to the difficulty in maximizingL(β,λ,S0∣x), we pose it as a Bayes parameter inference problem. The prior for β and λ is chosen to be noninformative, i.e., π(β)=1/β and π(λ)=1/λ. In addition S(0) is uniformly distributed on (r(t),N). By assumption of independence, the prior distribution of (β,λ,S0) is defined as (22)πβ,λ,S0∝1β1λIrt,NS0. In this way we get the posterior (23)πβ,λ,S0∣x∝πβ,λ,S0Lβ,λ,S0∣x. By drawing samples from π(β,λ,S(0)∣x) using a MCMC scheme, the Bayesian estimates λ^B, β^B, and S^B(0) of β,λ, and S(0) are obtained by the sample means of these simulating parameters, respectively.To examine the quality of these Bayesian estimators, a population of sizeN=20,000 with S0=12,000 initially susceptible is considered. The parameters are set as β=2.6, λ=28, d=5, Pτ=1/5(1≤τ≤5) and i(t)=10 for 1≤t≤5. The histograms of the estimates are presented in Figure 3 which shows that β^B,λ^B, and S^B(0) are suitable estimates of the original value of β,λ, and S0, respectively.Figure 3
Histograms of50,000 Bayesian estimators of β,λ, and S(0). The data are simulated from the model with parameters β=2.6,λ=28,S(0)=12,000, d=5,Pτ=1/5(1≤τ≤5), N=20,000, and i(t)=10 for 1≤t≤5.


## 4. Numerical Simulation

In this section, we make some contrasts on the models with and without media coverage.


### 4.1. Contrast on Fitting Simulated Data (All Data Used to Support Our Findings of This Study Are Simulated by Gillespie Algorithm)

We make a contrast on fitting models with and without media coverage to the data simulated from a model with media coverage. First a set of data are simulated from model (5) and (6) with parameters β=1.5,λ=30 and initial parameters Pτ=1/5(1≤τ≤5),i(t)=0(-4≤τ≤-1),i(0)=15,S(-5)=105. The scatter plot is illustrated in Figure 4. Next we estimated the parameters β,λ of model (5) and (6) and β in the model [12] by MLE, yielding β^=1.5265,λ^=27 in our model and β^=1.2847 in model [12]. Finally the corresponding deterministic models are obtained from the estimated parameters and the average infected levels are plotted in Figure 4. This simulation shows that model without consideration of media effect fails to capture the characteristics of data effected by media coverage.Figure 4
Different estimations on the same data produced by a model with media coverage whereβ=1.5 and λ=30,Pτ=1/5(1≤τ≤5),i(t)=0(-4≤τ≤-1),i(0)=15,S(-5)=105. The solid line corresponds to the estimated model with media coverage (the estimates are β^=1.5265 and λ^=27). The dash line corresponds to the estimated model without media coverage (the estimate is β^=1.2847).


### 4.2. Influence of Media Coverage on the Size of Newly Infected People

We now compare the number of newly infected individuals every day under differentλ′s. The solutions of the deterministic model (9) and (10) are plotted in Figure 5 with λ=0, λ=6, and λ=12. The simulation indicates that high involvement of media will significantly delay the epidemics peak and thus decreases the severity of the outbreak. This contrast shows that media coverage can decrease the prevalence of epidemic effectively so that we can take further measures to enhance the media coverage against the disease.Figure 5
Effect of media coverage on the epidemics outbreak. Comparison of epidemics withλ=0 (dash line), λ=6 (solid line), and λ=12 (dot solid line). In all these cases β=2,d=5,Pτ=1/5(1≤τ≤5), and the starting values are i(t)=0(-4≤τ≤-1),i(0)=15,S(-5)=105.


## 4.1. Contrast on Fitting Simulated Data (All Data Used to Support Our Findings of This Study Are Simulated by Gillespie Algorithm)

We make a contrast on fitting models with and without media coverage to the data simulated from a model with media coverage. First a set of data are simulated from model (5) and (6) with parameters β=1.5,λ=30 and initial parameters Pτ=1/5(1≤τ≤5),i(t)=0(-4≤τ≤-1),i(0)=15,S(-5)=105. The scatter plot is illustrated in Figure 4. Next we estimated the parameters β,λ of model (5) and (6) and β in the model [12] by MLE, yielding β^=1.5265,λ^=27 in our model and β^=1.2847 in model [12]. Finally the corresponding deterministic models are obtained from the estimated parameters and the average infected levels are plotted in Figure 4. This simulation shows that model without consideration of media effect fails to capture the characteristics of data effected by media coverage.Figure 4
Different estimations on the same data produced by a model with media coverage whereβ=1.5 and λ=30,Pτ=1/5(1≤τ≤5),i(t)=0(-4≤τ≤-1),i(0)=15,S(-5)=105. The solid line corresponds to the estimated model with media coverage (the estimates are β^=1.5265 and λ^=27). The dash line corresponds to the estimated model without media coverage (the estimate is β^=1.2847).


## 4.2. Influence of Media Coverage on the Size of Newly Infected People

We now compare the number of newly infected individuals every day under differentλ′s. The solutions of the deterministic model (9) and (10) are plotted in Figure 5 with λ=0, λ=6, and λ=12. The simulation indicates that high involvement of media will significantly delay the epidemics peak and thus decreases the severity of the outbreak. This contrast shows that media coverage can decrease the prevalence of epidemic effectively so that we can take further measures to enhance the media coverage against the disease.Figure 5
Effect of media coverage on the epidemics outbreak. Comparison of epidemics withλ=0 (dash line), λ=6 (solid line), and λ=12 (dot solid line). In all these cases β=2,d=5,Pτ=1/5(1≤τ≤5), and the starting values are i(t)=0(-4≤τ≤-1),i(0)=15,S(-5)=105.


## 5. Discussion

In this paper, we study the parameter estimation in a stochastic SIR age-of-infection epidemic model with media impact. The work can be generalized to other models with media coverage. We replace the age of infection by a stochastic recovery and suppose the probability of recovery to be proportional to the length of the subinterval of time. That is, an infected individual is supposed to recover from infection with probabilityγ1/m during any time interval Δt=1/m. Thus an infected individual does not recover within one day with probability (1-γ/m)m→e-γ as m→∞. Then the probability that an infected individual recovers within one day is qγ=1-e-γ. So, given the number of the infected I(t-1) at day t-1, the recovered number r(t) on day t is binomially distributed. That is, (24)rt~BIt-1,qγ.Then a more general SIR model with media coverage turns into (25)St=St-1-it,It=It-1-it+rt,Rt=Rt-1+rt,it~bSt-1,pt-1,rt~bIt-1,qγ, where pt-1=1-exp-(βexp-λIt-1/N/NI(t-1)),qγ=1-e-γ.Given infected data,x=(i(0),i(1),i(2),…,i(T)) and the likelihood is (26)Lβ,λ=Pid+1,…,iT∣i1,…,id=∏t=d+1TPit∣i1,…,it-1 where (27)Pit∣i1,…,it-1=Pit∣it-1∝pt-1it1-pt-1St-1-itFrom this formula, we know that the new added parameter γ does not affect the likelihood; thus parameters β and λ can be estimated in the similar way.We can also consider the case that the recovered individuals are not immune afterd days and immediately become susceptible to reinfection again. In such case, the model is adjusted to the following case: (28)St=St-1-it+it-d,it~BN-∑τ=1dPτit-τ,pt, where (29)pt=1-exp-βexp⁡-λi^t-1N-1∑τ=1dPτit-τ.This is the SI model concerning media coverage and age of infection and the parameters can be estimated likewise.



---

*Source: 3187807-2018-06-12.xml*