# The Basic Reproduction Number for the Markovian SIR-Type Epidemic Models: Comparison and Consistency

Authors: Muteb Faraj Alharthi
Journal: Journal of Mathematics
 (2022)

---

## Abstract

This paper is concerned with a well-known epidemiological concept to measure the spread of infectious disease, that is, the basic reproduction number. This paper has two major objectives. The first is to examine Bayesian sensitivity and consistency of this measure for the case of Markov epidemic models. The second is to assess the Martingale method by comparing its performance to that of Markov Chain Monte Carlo (MCMC) methods in terms of estimating this parameter and the infection and removal rate parameters given only removal data. We specifically consider the Markovian SIR (Susceptible-Infective-Removed) epidemic model routinely employed in the literature with exponentially distributed infectious periods. For illustration, numerical simulation studies are performed. Abakaliki smallpox data are examined as a real data application.


---

## Body

## 1. Introduction and Related Work

The basic reproduction number (or ratio) is a fundamental concept in epidemiology, commonly denoted asR0. For homogeneously mixing stochastic epidemic models, it can be roughly specified as the expected number of infections generated by one infectious individual in a large susceptible population [1, p. 11]. When R0>1, each generation of infectious individuals is more likely to infect the next with a larger number until the susceptible population is significantly depleted. Interventions in the infection or removal process, such as vaccination policies, are frequently carried out with the purpose of decreasing R0 below unity.One of the approaches to estimateR0 via fitting stochastic epidemic models to data, that offers a high level of modeling flexibility, is the data augmentation methods [2] such as MCMC method. Fitting stochastic epidemic models to data is a nonstandard problem [3], as data on the infection processes defined in such models is rarely observed directly. As a result, the likelihood of the observed data is intractable in the sense that calculating it is computationally expensive. Although the data-augmented approaches often weaken in large populations due to poor mixing and increasing computation time, they give a solution to this issue. MCMC methods have already been used successfully for model-based analyses of partially temporal data [4–7].For certain types of epidemic models (models that attain Markov property), the Martingale method [8, 9] can be utilized for inference. The essential concept behind this technique is to build Martingales for counting processes that are included in the epidemic mode from which estimating equations for the parameters of interest may be derived. This technique may be used to temporal and final size data and could be utilized to undertake nonparametric inference in specific scenarios. While this approach is very elegant, it is somewhat specialized and not as widely applicable as most other methods for fitting epidemic models to data. This paper aims to explore the efficiency of this approach in estimating model parameters as well as comparing its estimation performance with the MCMC method outcomes.When analysing temporal data, the likelihood of the observed data can become very difficult to evaluate, and so is the posterior distribution, as calculating the likelihood involves integration over all possible unobserved infection times, which is rarely analytically possible. However, simulating realizations from a stochastic epidemic model can be relatively straightforward. Implicit methods such as Approximate Bayesian Computation (ABC) may be suited to make inference for the parameters of epidemic models based on partially observed data, and this has been illustrated when temporal data are available in [10–13]. However, the choice of the tolerance parameter and the summary statistics are nontrivial tasks and may require expensive pilot runs to select the appropriate ones. Moreover, when working with large amounts of data, it often becomes necessary to use simplistic models due to computational restraints, especially when performing time-sensitive analyses. Approximation methods [14, 15] could be a natural solution under which likelihood can be computed. However, the simplifying assumptions used in approximation methods are not always appropriate in epidemic context.Our main focus in this paper is on the basic reproduction number, a quantity of central importance in mathematical epidemic theory, whose value essentially dictates whether or not a large epidemic outbreak can occur. As the emphasis is on the Martingale methods which need models to be Markovian, we specifically consider the SIR model that is routinely employed in the literature, namely, the SIR model with exponentially distributed infectious period. The motivations for this study are that Bayesian analysis of this important quantityR0 is of interest in its own right. In addition, the Markovian SIR model offers guidance that is valuable in analysing more complex and realistic models as SIR models themselves are frequently used as components of more complicated epidemic models, for instance those featuring populations divided into households [16] or epidemics on networks [1].A relative asymptotic efficiency ofR0 was derived, using Martingale method in [17] by comparing the variance of R0 with that of the maximum likelihood estimator based on complete observation of the epidemic. The asymptotic performance of a Martingale estimator for the infection rate was compared with the maximum likelihood estimator in [18]. However, these comparisons are of limited value since complete observation of epidemic is not possible. Large and finite sample properties for the Martingale estimator of R0 were studied in [19], where the performance of Martingale estimate with the maximum likelihood estimate based on complete observation and the maximum likelihood estimate based on the eventual size of the epidemic was compared.Initial values and continuous observation of the removal process are a plausible data set that may be observed during the duration of an epidemic. The time of diagnosis or commencement of symptoms, which is often accessible, might be used to approximate the time of removal for many diseases. Using such type of data, Martingale technique was applied in [9] to obtain estimating equations that enable the model parameters to be estimated separately. Data augmentation technique was implemented to deal with such data by [6] in which the missing data (the infection times) were treated as extra (unknown) parameters to be inferred from the data.Despite the extensive literature aboutR0, as far as the author knows, Bayesian sensitivity and consistency are not covered in much detail. Furthermore, given only the removal times, the performance of the Martingale and MCMC methods in terms of estimating R0 and the infection and removal rates is not assessed yet. This gap in the literature needs to be filled.The main purpose of this paper is to study Bayesian sensitivity and consistency for the posterior distribution ofR0 when Markovian SIR model is considered. Another important aim of this work is to explore the performance of the Martingale method in terms of estimating the basic reproduction number and the infection and removal rates, and compare it with the performance of the MCMC method performance as the most widely used method in epidemic inference.The rest of this paper is organized as follows: In Section2, relevant preliminary information will be presented which includes the definition of our stochastic epidemic model. Then, MCMC methods for the Markovian SIR model based on removal data will be introduced in Section 3. In Section 4, Martingale methods for the Markovian SIR model will be discussed. Bayesian consistency and sensitivity for R0 will be investigated in Section 5. A comparison study of Bayesian and Martingale estimation for R0 and the infection and removal rates will be conducted using both simulated and real data sets in Section 6. Some concluding remarks will be given in Section 7.


## 2. Preliminaries

### 2.1. The Stochastic SIR Epidemic Model

The stochastic SIR epidemic model is defined as follows [1, 6]. Consider a closed population of size N individuals divided to n and k initial susceptible and infected individuals, respectively. At any given time during the epidemic, each individual is either susceptible, infected, or removed. Susceptible individuals have not yet developed the disease, but they are at risk. Infected individuals carry the disease and can spread it to others. Individuals who have been removed are no longer infectious and have no role in disease transmission. The removed state is context-specific in applications, and it might equate to immunity, isolation, death, or other comparable outcomes. Each infective individual remains so for a period of time, called the infectious period, which is drawn from a specified nonnegative probability distribution TI. At the end of its infectious period, an individual is immediately removed. During its infectious period, a given infectious has contacts with a given susceptible in the population at times given by the points of a Poisson process of rate βn−1, where β is the infection rate parameter. The first such contact, if it occurs, results in the susceptible immediately becoming infected, and later contacts have no effect. All infectious periods and all Poisson processes are mutually independent. The epidemic ends as soon as there are no more infections remaining in the population.For any timet, let Xt,Yt, and Zt denote, respectively, the numbers of susceptible, infected, and removed individuals in the population. The above model is known as the Markovian SIR model when the infectious periods follow an exponential distribution with rate γ, say. In this case, the epidemic process Xt,Yt:t≥0 is a bivariate continuous time Markov chain with the following transition rates:(1)i,j⟶i−1,j+1:βn−1XtYt,i,j⟶i,j−1:γYt.The corresponding transition probabilities to an infection and removal are(2)PXt+δt−Xt=−1,Yt+δt−Yt=1|ℋt=βn−1XtYt+oδt,PXt+δt−Xt=0,Yt+δt−Yt=−1|ℋt=γYt+oδt,where all other transitions having probability oδt, ℋt is the sigma-algebra generated by the history of the process up to time t, and γ denotes the removal rate parameter.


### 2.2. The Basic Reproduction NumberR0

For the Markovian SIR model, the basic reproduction number is defined asR0=βETI=βγ−1 (see for example [1]). The importance of this measure comes from its relation to the spread of epidemic. In other words, in an infinitely large population, if R0≤1, then with probability one, only a finite number of individuals will become infected, i.e., the epidemic will die out, whereas when R0>1 there is a positive probability that the epidemic will take off. Consequently, and more importantly, knowing the value of R0 enables us to compute the proportion of a population that should be vaccinated in order to prevent an epidemic from spreading [1]. That is to say, vaccinating a fraction vc=1−1/R0 (critical vaccination coverage with perfect vaccine) of the community chosen uniformly at random is sufficient to prevent a large outbreak.


### 2.3. Epidemic Outbreak Data

Although the details of outbreak data vary considerably from study to study, broadly speaking, there are three types of epidemic data, namely, complete temporal data, partial (incomplete) temporal data, and final size data [20]. The complete temporal data consist of the time of infection and time of removal of all infected individuals in the population. The final size data do not contain explicit temporal information regarding the disease propagation throughout the population but instead consist of the initial number of susceptible individuals and which of these individuals were ultimately infected during the course of the epidemic. The partial temporal data are less detailed than complete temporal data, but more detailed than final size data. This type of data is more common in practice and typically consists of removal times.


## 2.1. The Stochastic SIR Epidemic Model

The stochastic SIR epidemic model is defined as follows [1, 6]. Consider a closed population of size N individuals divided to n and k initial susceptible and infected individuals, respectively. At any given time during the epidemic, each individual is either susceptible, infected, or removed. Susceptible individuals have not yet developed the disease, but they are at risk. Infected individuals carry the disease and can spread it to others. Individuals who have been removed are no longer infectious and have no role in disease transmission. The removed state is context-specific in applications, and it might equate to immunity, isolation, death, or other comparable outcomes. Each infective individual remains so for a period of time, called the infectious period, which is drawn from a specified nonnegative probability distribution TI. At the end of its infectious period, an individual is immediately removed. During its infectious period, a given infectious has contacts with a given susceptible in the population at times given by the points of a Poisson process of rate βn−1, where β is the infection rate parameter. The first such contact, if it occurs, results in the susceptible immediately becoming infected, and later contacts have no effect. All infectious periods and all Poisson processes are mutually independent. The epidemic ends as soon as there are no more infections remaining in the population.For any timet, let Xt,Yt, and Zt denote, respectively, the numbers of susceptible, infected, and removed individuals in the population. The above model is known as the Markovian SIR model when the infectious periods follow an exponential distribution with rate γ, say. In this case, the epidemic process Xt,Yt:t≥0 is a bivariate continuous time Markov chain with the following transition rates:(1)i,j⟶i−1,j+1:βn−1XtYt,i,j⟶i,j−1:γYt.The corresponding transition probabilities to an infection and removal are(2)PXt+δt−Xt=−1,Yt+δt−Yt=1|ℋt=βn−1XtYt+oδt,PXt+δt−Xt=0,Yt+δt−Yt=−1|ℋt=γYt+oδt,where all other transitions having probability oδt, ℋt is the sigma-algebra generated by the history of the process up to time t, and γ denotes the removal rate parameter.


## 2.2. The Basic Reproduction NumberR0

For the Markovian SIR model, the basic reproduction number is defined asR0=βETI=βγ−1 (see for example [1]). The importance of this measure comes from its relation to the spread of epidemic. In other words, in an infinitely large population, if R0≤1, then with probability one, only a finite number of individuals will become infected, i.e., the epidemic will die out, whereas when R0>1 there is a positive probability that the epidemic will take off. Consequently, and more importantly, knowing the value of R0 enables us to compute the proportion of a population that should be vaccinated in order to prevent an epidemic from spreading [1]. That is to say, vaccinating a fraction vc=1−1/R0 (critical vaccination coverage with perfect vaccine) of the community chosen uniformly at random is sufficient to prevent a large outbreak.


## 2.3. Epidemic Outbreak Data

Although the details of outbreak data vary considerably from study to study, broadly speaking, there are three types of epidemic data, namely, complete temporal data, partial (incomplete) temporal data, and final size data [20]. The complete temporal data consist of the time of infection and time of removal of all infected individuals in the population. The final size data do not contain explicit temporal information regarding the disease propagation throughout the population but instead consist of the initial number of susceptible individuals and which of these individuals were ultimately infected during the course of the epidemic. The partial temporal data are less detailed than complete temporal data, but more detailed than final size data. This type of data is more common in practice and typically consists of removal times.


## 3. MCMC Methods for the Markovian SIR Model Based on Removal Data

In the context of epidemic modeling, MCMC techniques are appealing for two reasons. First, they provide a great deal of modeling flexibility. Second, when used in combination with the Bayesian framework, they allow for the analysis of all model parameters or any function of them.Suppose that an epidemic outbreak which results in a total ofm removals in a closed population of size N, n initial susceptibles and one initial infectious is observed on t1I,tmR, so that all removals are assumed to be observed. Let 0=t1R,t2R,…,tmR be the observed ordered removal times and write tR=t1R,…,tmR. Also, let t1I<0 be the first unobserved infection time and tI=t2I,…,tmI be the rest of unobserved ordered infection times, where thI<th+1I<thR for h=1,2,…,m−1. Then, the joint likelihood of tI and tR given the model parameters and t1I is [6, 21].(3)LtI,tR|β,γ,t1I=∏j=2mβn−1XtjI−YIj−×∏j=1mγYtjR−×exp−∫t1ItmRβn−1XtYt+γYtdt,where tjI− denotes the time just prior to tjI and tjR−is defined similarly.In practice, the infection times are hard to observe which in turn makes the likelihood of observing only the removal timestR given the model parameters intractable. One possible solution to make the likelihood tractable is to use the data augmentation technique [2, 22] by treating the missing data as extra (unknown) parameters to be estimated from the data [6]. By doing so, we can sample from the conditional posterior distribution of the infection times in order to produce samples from the conditional posterior distributions of the model parameters.Now, we assign an independent gamma prior distribution for each rate parameter, namely,Gammaλζ,νζ, where ζ=β,γ, and assume a priori that −t1I∼Expθ, where Gammaa,b represents a gamma distribution with mean a/b and variance a/b2 and Expa denotes an exponential distribution with mean a−1. Then, by multiplying the likelihood and the priors, we obtain the following full conditional posterior densities for β,γ, and t1I, that are(4)β|γ,t1I,tI,tR∼Gammaλβ+m−1,νβ+n−1∫t1ItmRXtYtdt,γ|β,t1I,tI,tR∼Gammaλγ+m,νγ+∫t1ItmRYtdt,and(5)t1I|β,γ,tI,tR∼t2I−Expβ+γ+θ.Also, we have the following proportional full posterior density fortI:(6)πtI|β,γ,t1I,tR∝∏j=2mXtjI−YtjI−×∏j=1mYtjR−×exp−∫t1ItmRβn−1XtYt+γYtdt.The model parametersβ,γ, and the initial infection time t1I, all can be updated using Gibbs sampling steps as they have closed forms of the posterior density. However, the infection times vector, tI, needs to be updated using Metropolis-Hastings steps as its posterior density is not available explicitly. Any function of these parameters, such as the basic reproduction number R0=β/γ, can be obtained easily by sampling from the parameters posterior distributions. Hence, we can produce an MCMC estimate for the basic reproduction number, say R01^.


## 4. Martingale Methods for the Markovian SIR Model Based on Removal Data

Maximum likelihood estimation of parameters of the Markovian epidemic model is tedious when only partial information is available, as is often the case in practice. This motivates the development of easier methods of estimation such as the one established in [8] which can lead to simple but nevertheless efficient estimators.LetNt=Yt+Zt count the number of infections that occur by time t. Under the assumptions of the Markovian SIR epidemic model, there are two zero-mean Martingale processes [8] defined by(7)MβtmR−NtmR+βn−1∫0tmRXtYtdt=0,MγtmR−ZtmR+γ∫0tmRYtdt=0.Besides the Martingales considered above, the method of moments and the Martingale processXt1+R0/nZtmR−Zt:t≥0 were used in [9] to derive explicit estimators for the infection rate β and the mean infectious period γ−1 based on observing only the removal process. These estimators can be written as(8)β^=ZtmR∑j=1ZtmR−1n/n+1−jZtmRNtmR−∫0tmRZtdt−XtmR∫0tmR1+∑j=1ZtmR−1n/n+1−jZtmR/nZtmR−Ztdt,γ−1^=NtmR−∫0tmRZtdt−XtmR∫0tmR1+∑j=1ZtmR−1n+1−jZtmR/nZtmR−ZtdtZtmR,provided at least one individual remains susceptible at time tmR. Note that these estimators depend only on the initial number of susceptible individuals n, the final number of removals ZtmR, and the removal times. A Martingale estimate for the basic reproduction number is given by(9)R02^=∑j=1ZtmR−1nn+1−jZtmR,with standard deviation given as(10)s.dR02^=1ZtmRZtmRR02^2+n2∑j=1ZtmR−11n+1−j2.


## 5. Bayesian Consistency and Sensitivity forR01

An estimatorθn of a quantity of interest θ is said to be consistent if θn⟶θ as n⟶∞, in other words, if it is guaranteed to converge to the true value. However, the posterior is consistent for θ if the posterior distribution concentrates around a mass point which is the true value of the parameter [23].Given complete observation of the epidemic, the following explicit posterior density forR01 was derived in [24] using prior settings as in Section 3.(11)πR01|tI,tR,t1I=Γ2m−1+λβ+λγΓm+λγΓm−1+λβ×n−1A+νβ/B+νγm−1+λβR01m−2+λβ1+n−1A+νβ/B+νγR012m−1+λβ+λγ,given that R01>0, where A=∫t1ItmRXtYtdt and B=∫t1ItmRYtdt.The posterior mean and variance are given as follows:(12)ER01|tI,tR,t1I=m−1+λβm−1+λγ×B+νγn−1A+νβ,VarR01|tI,tR,t1I=2m−1+λβ+λγm−1+λβm−2+λγm−1+λγ2×B+νγn−1A+νβ2,given that m+λγ>2.Common choices of the prior parameters in the literature areλβ=λγ=λ and νβ=νγ=ν, where λ≥1 and ν is a small positive number. This prior setting yields(13)πR01|tI,tR,t1I=Γ2m−1+2λΓm+λΓm−1+λ×n−1A+ν/B+νm−1+λR01m−2+λ1+n−1A+ν/B+νR012m−1+2λ,with(14)ER01|tI,tR,t1I=B+νn−1A+ν,(15)VarR01|tI,tR,t1I=2m−2+λ×B+νn−1A+ν2.It is of our interest to study the behaviour of posterior density, the posterior mean, and the posterior variance in two cases. The first case is when the initial number of susceptiblesn tends to infinity; and the second case is when the prior distributions of the model parameters (β and γ) get diffuse. We restrict our results on major outbreak situation [18] as it may not be of importance to model and investigate minor outbreak in practical situations.


### 5.1. Consistency forR01

Suppose that we have a sequence of processesXnt,Ynt indexed by X0=n, the initial number of susceptibles. Without the loss of generality, we set t1I=0 and let tmR=τ be the time when the epidemic ends.We need to examine the behaviour of the following random variables asn⟶∞(16)mn,An=∫0τnXntYntdt,andBn=∫0τnYntdt.It was shown in [18, 25] that as n⟶∞,mn/n⟶Pρ, where ρ is a constant satisfies 1−ρ=e−ρR0.Also, it was proven that(17)n−1∫0τnYntdt⟶Pρ/γandn−2∫0τnXntYntdt⟶Pρ/β.To indicate that the posterior mean is specific ton, we rewrite (14) as(18)EnR01|tI,tR,t1I=Bn+νn−1An+ν=Bn/n+ν/nn−2An+ν/n⟶ρ/γρ/β=βγ=R01,asn⟶∞.Similarly, to indicate that the posterior variance is specific ton, we rewrite (15) as(19)VarnR01|tI,tR,t1I=2mn−2+λ×Bn+νn−1An+ν2=2/nmn/n−2/n+λ/n×Bn/n+ν/nn−2An+ν/n2⟶0asn⟶∞.Therefore, we conclude that as we have more and more data, the posterior mean tends to the true value and the posterior variance tends to zero and so the posterior distribution concentrates around a mass point which is the true value of the parameter. Mathematically, asn⟶∞ we have(20)πR01|tI,tR,t1I⟶δR01˜R01,where δx˜x is a Dirac delta function and x˜ is the true value of x.That is to say, the posterior is consistent atR01˜, and the Bayes estimate (the posterior mean) is a consistent estimator for R01 as n⟶∞.


### 5.2. Sensitivity forR01

To see the effect of the prior’s diffuseness on both posterior mean and variance ofR01, we reformulate them in terms of the mean and variance of the prior distribution, that are(21)EGammaλ,ν=λ/ν=μand VarGammaλ,ν=λ/ν2=σ2,which implies λ=μ2/σ2 and ν=μ/σ2. Substituting these values into (14) and (15) gives(22)ER01|tI,tR,t1I=B+μ/σ2n−1A+μ/σ2⟶nBA,asσ2⟶∞.This is in accord with the maximum likelihood estimator ofR0 given by the ratio of the maximum likelihood estimators of β and γ, that is, R0^=nm−1B/mA (see for example [1]).Also, we have(23)VarR01|tI,tR,t1I=2m−2+μ2/σ2×B+μ/σ2n−1A+μ/σ22⟶2n2B2m−2A2,asσ2⟶∞.


## 5.1. Consistency forR01

Suppose that we have a sequence of processesXnt,Ynt indexed by X0=n, the initial number of susceptibles. Without the loss of generality, we set t1I=0 and let tmR=τ be the time when the epidemic ends.We need to examine the behaviour of the following random variables asn⟶∞(16)mn,An=∫0τnXntYntdt,andBn=∫0τnYntdt.It was shown in [18, 25] that as n⟶∞,mn/n⟶Pρ, where ρ is a constant satisfies 1−ρ=e−ρR0.Also, it was proven that(17)n−1∫0τnYntdt⟶Pρ/γandn−2∫0τnXntYntdt⟶Pρ/β.To indicate that the posterior mean is specific ton, we rewrite (14) as(18)EnR01|tI,tR,t1I=Bn+νn−1An+ν=Bn/n+ν/nn−2An+ν/n⟶ρ/γρ/β=βγ=R01,asn⟶∞.Similarly, to indicate that the posterior variance is specific ton, we rewrite (15) as(19)VarnR01|tI,tR,t1I=2mn−2+λ×Bn+νn−1An+ν2=2/nmn/n−2/n+λ/n×Bn/n+ν/nn−2An+ν/n2⟶0asn⟶∞.Therefore, we conclude that as we have more and more data, the posterior mean tends to the true value and the posterior variance tends to zero and so the posterior distribution concentrates around a mass point which is the true value of the parameter. Mathematically, asn⟶∞ we have(20)πR01|tI,tR,t1I⟶δR01˜R01,where δx˜x is a Dirac delta function and x˜ is the true value of x.That is to say, the posterior is consistent atR01˜, and the Bayes estimate (the posterior mean) is a consistent estimator for R01 as n⟶∞.


## 5.2. Sensitivity forR01

To see the effect of the prior’s diffuseness on both posterior mean and variance ofR01, we reformulate them in terms of the mean and variance of the prior distribution, that are(21)EGammaλ,ν=λ/ν=μand VarGammaλ,ν=λ/ν2=σ2,which implies λ=μ2/σ2 and ν=μ/σ2. Substituting these values into (14) and (15) gives(22)ER01|tI,tR,t1I=B+μ/σ2n−1A+μ/σ2⟶nBA,asσ2⟶∞.This is in accord with the maximum likelihood estimator ofR0 given by the ratio of the maximum likelihood estimators of β and γ, that is, R0^=nm−1B/mA (see for example [1]).Also, we have(23)VarR01|tI,tR,t1I=2m−2+μ2/σ2×B+μ/σ2n−1A+μ/σ22⟶2n2B2m−2A2,asσ2⟶∞.


## 6. Comparison of MCMC and Martingale Estimation forR0, β, and γ

In this section, we explored, via simulation study, the accuracy of Martingale method for estimatingR02, β, and γ given only a set of removal times. Moreover, a comparison between MCMC and Martingale estimation methods was performed. Finally, a widely studied removal data set obtained from a smallpox outbreak [26] was considered.We simulated epidemics under a Markovian SIR model. All outbreaks were in a closed population consisting of different initially susceptible individualsn and a single initially infectious case k=1. In all simulations, the value of γ was taken (without loss of generality) to be equal to 1, corresponding to a mean infectious period of 1 day, say.To make our comparisons more useful, we wish to eliminate realizations that resulted in early extinction, thereby presenting results based on major epidemics only. A working algorithm that excluded all simulations for whichRτ<0.1n+k was implemented in [27], where τ is the extinction time. Here, we followed [4, 21] and excluded simulated epidemics in which the final size was not consistent with the theoretical probability of a major epidemic. In other words, we simulated the model until we have the desired number of realizations that have final sizes fit well within the major outbreak part of the simulation-based final size distribution.


### 6.1. Simulation Study I

There have been several applications of Martingale results to issues of epidemiological parameter inference, but the investigation of their efficacy has been unaddressed so far. To see how efficient the Martingale method is in estimatingR02, β, and γ separately given a set of removal times, we conducted this simulation study using different population and final epidemic sizes. We set β=2.5 and γ=1, so that R02=2.5. 1000 major outbreak realizations were simulated from the model. In Table 1, the column labelled “Avge ZtmR” gives the average number of removals over the realizations. The other columns represent the average of the Martingale method outcomes for estimating R02 and its standard deviation, β, and γ. It can be seen that the method is very efficient in recovering true values and best results are achieved as both population and final epidemic size get large.Table 1
Results of simulation study I withR02=2.5, β=2.5, and γ=1. Each row is obtained based on the average of 1000 major outbreak realizations.
NAvgeZtmRAvge (Martingale outcomes)R02^s.d^R02^β^γ−1^50422.18380.552.33481.0087100892.54140.402.29141.16083002682.54010.212.36851.09277006252.50800.162.43501.038210008922.50820.132.44331.0322150013392.50750.102.46461.0214


### 6.2. Simulation Study II

To make a reasonable and meaningful comparison between MCMC and Martingale methods for estimating the basic reproduction number and infection and removal rates given only removal data, we simulated 10 typical (in the sense that the final size fits within the high density of the major outbreak component of the simulation-based final size distribution) epidemic realizations from the model for different population sizes as seen in Table2 using various values of the threshold parameter R0 (1.5,2, and 3). Each outbreak simulation ended when all the infections were removed. For each outbreak, we computed both MCMC and Martingale estimates and then average over the 10 realizations. In Table 2, the column labelled “Avge ZtmR” gives the average number of removals over the 10 realizations. Martingale results were obtained using 8, 9, and 10, whereas MCMC results were produced using standard algorithms (see for example [6, 21]).Table 2
Results of simulation study II for different values ofR0, where γ value is set to 1. Each row is obtained based on the average of 10 major outbreak realizations.
NTrueR0AvgeZtmRAvge (MCMC outcomes)Avge (Martingale outcomes)R01^s.d^β^γ−1^R02^s.d^β^γ−1^501.5331.7048 (0.43)1.85461.10421.5695 (0.40)2.09700.92072412.2245 (0.52)2.80480.96582.0556 (0.49)2.32631.08763463.0847 (0.74)3.39960.95052.7089 (0.66)2.48421.10671001.5671.6877 (0.30)1.77771.02191.6297 (0.29)1.69720.98582812.0881 (0.34)2.37130.96732.0321 (0.34)1.91761.13313912.8392 (0.48)3.69820.80262.7349 (0.47)2.87290.96353001.51821.5567 (0.17)1.64601.00221.5313 (0.16)1.70520.905722382.0055 (0.19)2.16830.98031.9856 (0.19)2.00971.025932802.9007 (0.28)3.29390.89492.9288 (0.29)2.88861.0316The findings of this simulation study can be summarized as follows: For both MCMC and Martingale outcomes, the estimation improved as both the population and epidemic final size increase. When the population size gets large, the Martingale method seems to be useful in terms of computational times and precision. When looking to standard deviation as a measure of discrepancy, MCMC results were as good as the Martingale ones. However, as pointed out in [28], in large-scale settings, MCMC methods struggle due to the increased computational burden as a result of high dimensional parameter space.


### 6.3. Abakaliki Smallpox Data

We now consider a widely studied temporal data set obtained from a smallpox outbreak in a closed population ofN=120 individual (n=119, k=1) in Abakaliki, Nigeria [26], p. 125. The data consist of times between removals measured in days, which are(24)13,7,2,3,0,0,1,4,5,3,2,0,2,0,5,3,1,4,0,1,1,1,2,0,1,5,0,5,5,where zero indicates that the case occurred at the same day as the previous one.Therefore, in our notation, the observed removal times are given by(25)tR=0,13,20,22,25,25,25,26,30,35,38,40,40,42,42,47,50,51,55,55,56,57,58,60,60,61,66,66,71,76.The use of the SIR stochastic epidemic to model these data is not entirely appropriate, since smallpox has an appreciable latent period. However, it can be used for illustrative purpose [6]. In the literature, the basic reproductive number for smallpox is about 3−6 [29]. However, using the Martingale method, the authors of [8, 9] gave an estimate of 1.1 for R02, while MCMC method produces an estimate of R01=1.13 in [6]. Using the ABC method, an estimate of 1.1 for the basic reproduction number parameter was obtained in [10].Surprisingly, no study in the literature has provided Martingale estimates for the infection rateβ and the mean infectious period γ−1 based on these removal times. Using (8), we obtained the following Martingale estimates, β^=0.335 and γ−1^=3.292. It is of concern why these estimates are not somehow similar to those obtained using different methods. However, it was pointed out in [14] that assuming that all the members are susceptible and are mixed homogeneously may cause a significant underestimation of the basic reproduction number and hence related parameters. Table 3 shows estimates of β^, γ−1^ and R0^ for Abakaliki smallpox data, using different approaches.Table 3
Estimates ofβ^, γ−1^ and R0^ for Abakaliki smallpox data, using different approaches.
Methodβ^γ−1^R0^Martingale0.3353.2921.1MCMC0.09911.4161.13ABC0.11101.1


### 6.4. Comparison Discussion

The simulation studies above revealed the accuracy and the efficiency of the Martingale method when estimating the basic reproduction number and the infection and removal rates of the Markovian SIR model. Although the Martingale method seems promising in estimating all the parameters (see Table1) and preforms well with large population sizes, it is only applicable when epidemic models are Markov. Unfortunately, most sophisticated epidemic models need to be non-Markov to mimic the mechanism of infectious diseases.The MCMC method, on the other hand, estimated all the parameters with very good precision. It is known that this method struggles when the population size gets large [28] and for this reason, we considered only N=50,100, and 300, to make a fair comparison.In combination with the Bayesian approach, MCMC method has no restriction that models be Markov. Therefore, it is an attractive choice in the epidemic context as it permits a huge amount of modeling flexibility, and it enables analysis of all of the model parameters, or any function of them.


## 6.1. Simulation Study I

There have been several applications of Martingale results to issues of epidemiological parameter inference, but the investigation of their efficacy has been unaddressed so far. To see how efficient the Martingale method is in estimatingR02, β, and γ separately given a set of removal times, we conducted this simulation study using different population and final epidemic sizes. We set β=2.5 and γ=1, so that R02=2.5. 1000 major outbreak realizations were simulated from the model. In Table 1, the column labelled “Avge ZtmR” gives the average number of removals over the realizations. The other columns represent the average of the Martingale method outcomes for estimating R02 and its standard deviation, β, and γ. It can be seen that the method is very efficient in recovering true values and best results are achieved as both population and final epidemic size get large.Table 1
Results of simulation study I withR02=2.5, β=2.5, and γ=1. Each row is obtained based on the average of 1000 major outbreak realizations.
NAvgeZtmRAvge (Martingale outcomes)R02^s.d^R02^β^γ−1^50422.18380.552.33481.0087100892.54140.402.29141.16083002682.54010.212.36851.09277006252.50800.162.43501.038210008922.50820.132.44331.0322150013392.50750.102.46461.0214


## 6.2. Simulation Study II

To make a reasonable and meaningful comparison between MCMC and Martingale methods for estimating the basic reproduction number and infection and removal rates given only removal data, we simulated 10 typical (in the sense that the final size fits within the high density of the major outbreak component of the simulation-based final size distribution) epidemic realizations from the model for different population sizes as seen in Table2 using various values of the threshold parameter R0 (1.5,2, and 3). Each outbreak simulation ended when all the infections were removed. For each outbreak, we computed both MCMC and Martingale estimates and then average over the 10 realizations. In Table 2, the column labelled “Avge ZtmR” gives the average number of removals over the 10 realizations. Martingale results were obtained using 8, 9, and 10, whereas MCMC results were produced using standard algorithms (see for example [6, 21]).Table 2
Results of simulation study II for different values ofR0, where γ value is set to 1. Each row is obtained based on the average of 10 major outbreak realizations.
NTrueR0AvgeZtmRAvge (MCMC outcomes)Avge (Martingale outcomes)R01^s.d^β^γ−1^R02^s.d^β^γ−1^501.5331.7048 (0.43)1.85461.10421.5695 (0.40)2.09700.92072412.2245 (0.52)2.80480.96582.0556 (0.49)2.32631.08763463.0847 (0.74)3.39960.95052.7089 (0.66)2.48421.10671001.5671.6877 (0.30)1.77771.02191.6297 (0.29)1.69720.98582812.0881 (0.34)2.37130.96732.0321 (0.34)1.91761.13313912.8392 (0.48)3.69820.80262.7349 (0.47)2.87290.96353001.51821.5567 (0.17)1.64601.00221.5313 (0.16)1.70520.905722382.0055 (0.19)2.16830.98031.9856 (0.19)2.00971.025932802.9007 (0.28)3.29390.89492.9288 (0.29)2.88861.0316The findings of this simulation study can be summarized as follows: For both MCMC and Martingale outcomes, the estimation improved as both the population and epidemic final size increase. When the population size gets large, the Martingale method seems to be useful in terms of computational times and precision. When looking to standard deviation as a measure of discrepancy, MCMC results were as good as the Martingale ones. However, as pointed out in [28], in large-scale settings, MCMC methods struggle due to the increased computational burden as a result of high dimensional parameter space.


## 6.3. Abakaliki Smallpox Data

We now consider a widely studied temporal data set obtained from a smallpox outbreak in a closed population ofN=120 individual (n=119, k=1) in Abakaliki, Nigeria [26], p. 125. The data consist of times between removals measured in days, which are(24)13,7,2,3,0,0,1,4,5,3,2,0,2,0,5,3,1,4,0,1,1,1,2,0,1,5,0,5,5,where zero indicates that the case occurred at the same day as the previous one.Therefore, in our notation, the observed removal times are given by(25)tR=0,13,20,22,25,25,25,26,30,35,38,40,40,42,42,47,50,51,55,55,56,57,58,60,60,61,66,66,71,76.The use of the SIR stochastic epidemic to model these data is not entirely appropriate, since smallpox has an appreciable latent period. However, it can be used for illustrative purpose [6]. In the literature, the basic reproductive number for smallpox is about 3−6 [29]. However, using the Martingale method, the authors of [8, 9] gave an estimate of 1.1 for R02, while MCMC method produces an estimate of R01=1.13 in [6]. Using the ABC method, an estimate of 1.1 for the basic reproduction number parameter was obtained in [10].Surprisingly, no study in the literature has provided Martingale estimates for the infection rateβ and the mean infectious period γ−1 based on these removal times. Using (8), we obtained the following Martingale estimates, β^=0.335 and γ−1^=3.292. It is of concern why these estimates are not somehow similar to those obtained using different methods. However, it was pointed out in [14] that assuming that all the members are susceptible and are mixed homogeneously may cause a significant underestimation of the basic reproduction number and hence related parameters. Table 3 shows estimates of β^, γ−1^ and R0^ for Abakaliki smallpox data, using different approaches.Table 3
Estimates ofβ^, γ−1^ and R0^ for Abakaliki smallpox data, using different approaches.
Methodβ^γ−1^R0^Martingale0.3353.2921.1MCMC0.09911.4161.13ABC0.11101.1


## 6.4. Comparison Discussion

The simulation studies above revealed the accuracy and the efficiency of the Martingale method when estimating the basic reproduction number and the infection and removal rates of the Markovian SIR model. Although the Martingale method seems promising in estimating all the parameters (see Table1) and preforms well with large population sizes, it is only applicable when epidemic models are Markov. Unfortunately, most sophisticated epidemic models need to be non-Markov to mimic the mechanism of infectious diseases.The MCMC method, on the other hand, estimated all the parameters with very good precision. It is known that this method struggles when the population size gets large [28] and for this reason, we considered only N=50,100, and 300, to make a fair comparison.In combination with the Bayesian approach, MCMC method has no restriction that models be Markov. Therefore, it is an attractive choice in the epidemic context as it permits a huge amount of modeling flexibility, and it enables analysis of all of the model parameters, or any function of them.


## 7. Concluding Remarks

In this paper, we have considered a quantity of central importance in mathematical epidemic theory, that is, the basic reproduction number. Bayesian sensitivity and consistency of its posterior distribution were discussed, particularly in the case where the infectious period is exponentially distributed as this case is commonly used by modellers, partly for mathematical convenience. The performance of the Martingale method in terms of estimating the basic reproduction number and the infection and removal rates was explored and compared to the performance of the most widely used method, that is, the MCMC method.The Martingale method’s key advantages were that it reduced the stochastic epidemic model’s unreasonable assumptions and it simplified the computations. However, this approach seems to be undesirable as it requires models to be Markov. The MCMC method, on the other hand, allows for a great deal of modeling adaptability. In most cases, no such modeling constraint on infectious periods is required. It also allows analysis of all model parameters, or any function of them, when used in conjunction with the Bayesian technique.This work is significant in the sense that it explores Bayesian sensitivity and consistency for the posterior distribution of the basic reproduction number in the case of the Markovian SIR model. In addition, it showed the interesting performance of the Martingale method in terms of estimating the basic reproduction number and the infection and removal rates.It is natural to study Bayesian sensitivity and consistency for other models with different choices of infectious period distribution. However, the posterior distribution of the basic reproduction number given a complete observation may not have a closed form, which in turn makes the analysis harder.



---

*Source: 1925202-2022-09-24.xml*