Monitoring SEIRD model parameters using MEWMA for the COVID-19 pandemic with application to the State of Qatar
Edward L. Boone, Abdel-Salam G. Abdel-Salam, Indranil Sahoo, Ryad Ghanam, Xi Chen, Aiman Hanif
aa r X i v : . [ s t a t . M E ] J a n ARTICLE
Monitoring SEIRD model parameters using MEWMA for theCOVID-19 pandemic with application to the State of Qatar.
E. L. Boone a , Abdel-Salam G. Abdel-Salam b , Indranil Sahoo a , Ryad Ghanam c , XiChen d , and Aiman Hanif a a Department of Statistical Sciences and Operations Research, Virginia CommonwealthUniversity, Richmond, Virginia, USA; b Department of Mathematics, Statistics and Physics,College of Arts and Sciences, Qatar University, Doha, Qatar. [email protected]; c Departmentof Liberal Arts and Science, Virginia Commonwealth University in Qatar, Doha,[email protected]; d Grado Department of Industrial and Systems Engineering,Virginia Tech, Blacksburg, Virginia, USA.
ARTICLE HISTORY
Compiled February 1, 2021
ABSTRACT
During the current COVID-19 pandemic decision makers are tasked with implement-ing and evaluating strategies for both treatment and disease prevention. In orderto make effective decisions they need to simultaneously monitor various attributesof the pandemic such as transmission rate and infection rate for disease prevention,recovery rate which indicates treatment effectiveness as well as the mortality rateand others. This work presents a technique for monitoring the pandemic by employ-ing an Susceptible, Exposed, Infected, Recovered Death model regularly estimatedby an augmented particle Markov chain Monte Carlo scheme in which the posteriordistribution samples are monitored via Multivariate Exponentially Weighted Aver-age process monitoring. This is illustrated on the COVID-19 data for the State ofQatar.
KEYWORDS
Epidemiology; Augmented particle Markov chain Monte Carlo; MultivariateExponentially Weighted Moving Average; process monitoring; COVID-19
1. Introduction
Coronavirus Disease 2019 (COVID-19) [22, 27] is a severe pandemic affecting thewhole world with a fast spreading regime, requiring to perform strict precautions tokeep it under control. As there are limited cure and target treatment at the moment,establishing those precautions become inevitable. These limitations [10] can be listedas social distancing, closure of businesses / schools and travel prohibitions [4].Corona Virus is a new human Betacoronavirus that uses densely glycosylated spikeprotein to penetrate host cells. The COVID-19 belongs to the same family classificationwith Nidovirales, viruses that use a nested set of mRNAs to replicate and it further fallsunder the subfamily of alpha, beta, gamma and delta Co-Vis. The virus that causesCOVID-19 belongs to the Betacoronavirus 2B lineage and has a close relationship with
CONTACT E. L. Boone. Email: [email protected]
ARS species. It is a novel virus since the monoclonal antibodies do not exhibit a highdegree of binding to SARS-CoV-2. Replication of the viral RNA occurs when RNApolymerase binds and re-attaches to multiple locations [8, 18].Cases of COVID-19 started in December 2019 when a strange condition was reportedamong a cluster of patients in Wuhan, China. Within a few weeks of this, the COVID-19 virus had spread to different parts of the world. On January 20, 2020, the first case ofCOVID-19 was recorded in the United States; Italy reported its first confirmed case onJanuary 31, 2020. With COVID-19 cases rising across the world, the governments weresoon seen intervening in financial and healthcare sectors. In late January, 2020, thefirst U.S. travel restrictions were imposed on travel from China. Weeks later, additionaltravel bans were imposed on countries in Europe and the United Kingdom. The WorldHealth Organization (WHO) declared COVID-19 a pandemic on March 11, 2020, witha total of more than 100,000 cases globally. As of January 26, 2021, the worldwidetotal number of confirmed COVID-19 cases eclipses 100 million with over 2.15 milliondeaths. This virus has a global mortality rate of 3.4%, which makes it more severethan the flu. The elderly who have other pre-existing illnesses are succumbing moreto the COVID-19. People with only mild symptoms recover within 3 to 7 days, whilethose with conditions such as pneumonia or severe diseases take weeks to recover. Asof January 26, 2021, the recovery percentage of patients, for example, in China standsat 95%. The global recovery percentage rate of COVID-19 is somewhere around 97%[26].The main efforts in the literature is focused on model estimation and forecastingthe dynamic nature of the COVID-19 pandemic versus monitoring the process. TheSusceptible, Exposed, Infected, Recovered, Death (SEIRD) model is a common com-partmental model used for modeling disease through a population and variants of ithave been used in modeling the COVID-19 pandemic, such as [15] for Italy, [2] forIndia and [9] for the State of Qatar. For more on disease modeling in general, see[6, 13, 17, 24].Other modeling approaches include a time-series model to analyze the outbreak ofthe pandemic [7], a time-varying Bayesian semi-parametric model to look at short-term projections of the pandemic [23]. [12] studies the dynamic pattern of COVID-19deaths over time. The impact of government intervention such as a lockdown hasbeen studied for China in [25] and for India in [2, 3]. Monitoring of the mathematicalprocess has been done for Ukraine’s COVID-19 outbreak [14]; however the parameterestimates are not generated directly from the data and no clear monitoring schemeis presented. As the vaccine is beginning to be distributed to the public, monitoringthe pandemic is more important to decision makers as they can determine how thepandemic is progressing as well as any shocks to the system that may be problematic.The monitoring approach needs to be able to react quickly to any large shifts in thesystem.The goal of this work is to develop a method of monitoring a pandemic using a basemathematical model, such as SEIRD, that can be quickly updated as soon as newinformation comes in and can “signal” if there is a change in the parameters of themathematical model. The literature does not seem to provide any approaches that meetthis goal. The big challenge in this problem is updating the parameters in an automatedway and converting those parameters to something that can be monitored. Here aBayesian approach is taken for the parameter estimation via a sampling algorithm thatwill allow for quick updating, avoid particle depletion and from which the samples canbe monitored using a standard process control regime.This work is organized in the following manner. Section 2 introduces the SEIRD2odel specific to the State of Qatar, the mean model, and the likelihood used. Sec-tion 3 introduces the Reproduction number and illustrates its deficiency in this case.This is followed by the Markov chain Monte Carlo sampling algorithm with particleaugmentation to update the parameters at each time step in Section 4. The Multi-variate Exponentially Weighted Moving Average (MEWMA) monitoring approach ispresented in Section 5. The method is illustrated on real data from the State of Qatarin Section 6, which shows how the monitoring can be employed to identify criticalchanges in the pandemic. Finally, Section 7 provides a discussion of the method andprovides some suggestions for implementation as well as possible areas for improve-ment.
2. SEIRD Model
Let t be a time index that is the number of days from the first recorded case ofCOVID-19 in the population of interest, S ( t ) be the number of subjects Susceptibleat time t , E ( t ) be the number of subjects Exposed at time t, I ( t ) be the number ofInfected (symptomatic) subjects at time t , R ( t ) be the cumulative number of Recoveredsubjects at time t and D ( t ) be the cumulative number of subject Deaths at time t .This can be modeled with the following system of ordinary differential equations: dS ( t ) dt = − αS ( t ) E ( t ) dE ( t ) dt = αS ( t ) E ( t ) − βE ( t ) − γE ( t ) dI ( t ) dt = βE ( t ) − γI ( t ) − ηI ( t ) dR ( t ) dt = γI ( t ) dD ( t ) dt = ηI ( t ) , (1)where the parameters are explained as follows: α is the transmission rate (per day × individual ) from Susceptible to Exposed, the rate (per day) at which Exposed becomeInfected (symptomatic) is denoted by β , γ is the rate (per day) at which Infectedbecome recovered and the mortality rate (per day) for those Infected is denoted by η . Notice that, this model formulation makes several key assumptions which are asfollows. Immigration, emigration, natural mortality and births are negligible over thetime frame and hence are not in the model. Once a person is in the Infected group,they are quarantined and hence they do not mix with the Susceptible population. TheRecovered and Deaths compartments are for those who first are infected. Those whoare exposed and are asymptomatic recover at the same rate γ as those who become sickand recover. The SEIRD model presented here is the same as the one presented in [9]and matches the assumptions needed for the example provided in Section 6, however,the estimation and monitoring method is not specific to this particular model.Due to the dynamic nature of how the pandemic has developed, assuming the systemis in “steady state” is invalid as governments have intervened into the system in aneffort to influence various parameters as well as medical treatment of the disease haschanged across the time frame. Hence the parameters are also functions of time denoted3s α ( t ) , β ( t ) , γ ( t ) and η ( t ). Specifically, dλ S ( t ) dt = − α ( t ) λ S ( t ) λ E ( t ) dλ E ( t ) dt = α ( t ) λ S ( t ) λ E ( t ) − β ( t ) λ E ( t ) − γ ( t ) λ E ( t ) dλ I ( t ) dt = β ( t ) λ E ( t ) − γ ( t ) λ I ( t ) − η ( t ) λ I dλ R ( t ) dt = γ ( t ) λ I ( t ) dλ D ( t ) dt = η ( t ) λ I ( t ) , (2)where λ s ( t ), λ E ( t ), λ s ( I ), λ R ( t ) and λ D ( t ) denote the respective mean parameters.At each time point t the parameters must have a prior distribution. For this work theprior distribution specification will be the same for all t , however, this is not necessaryif one has information that needs to be included at a specific time.Since a Bayesian methodology is being employed the likelihood is specified to be: I ( t ) ∼ P oisson ( λ I ( t )) R ( t ) ∼ P oisson ( λ R ( t )) D ( t ) ∼ P oisson ( λ D ( t )) . (3)Notice that S ( t ) and E ( t ) are not in the likelihood as they are latent states in thatthey are not directly observed. The true likelihood for { S ( t ) , E ( t ) , I ( t ) , R ( t ) , D ( t ) } should be Multinomial. However, with two latent states, one of which is the largeststate, the Multinomial approach is challenging to apply, thus this work uses adoptsPoisson likelihood as an approximation.
3. The Basic Reproduction Number R The basic reproduction number, R , is defined as the expected number of secondarycases produced by a single infection in a completely susceptible population. It is di-mensionless and can be calculated as a product of the transmissibility, the averagerate of contact between susceptible and infected individuals, and the duration of theinfectiousness [5]. In model (1), the last two equations do not contribute to R and so R = αβ + γ . Since our model is time varying, it follows that R ( t ) = α ( t ) β ( t ) + γ ( t )In this work the Reproduction number is not considered since it does not accountfor all the parameters in the model. One of the goals of this work is to ensure themonitoring process accounts for all of the parameters.4 . Sequential Sampling with Particle Augmentation From a monitoring perspective, we need to essentially estimate the parameters, θ t = { α ( t ) , β ( t ) , γ ( t ) , η ( t ) } ), at each time step dependent primarily on the latest data. If attime t , we want to estimate the change in the parameters from t − t , we need anestimation algorithm that can update the parameters so that changes can be detected.In a dynamic system where the parameters are changing due to interventions into thesystem such as vaccines or quarantines, the estimation algorithm needs to also takeinto account the previous changes that have taken place in the system. Traditionalapproaches to modeling varying coefficient models from a Bayesian perspective havedifficulties as the sampling methods required often suffer from particle depletion as theprocess proceeds through time.The proposed Algorithm 1 is a variant of the Sampling Importance Resamplingalgorithm at each time step. To ensure there are enough particles to work through thesampling process, the basic idea is to augment each accepted sample by some randomperturbations to generate new particles to move through to the next step.Let us fix the notation first. Denote D ( k ) = { S ( k ) , E ( k ) , I ( k ) , R ( k ) , D ( k ) } as theactual state vector of the system at time step k , where the first two components arelatent, ∀ k ∈ { , , . . . , T } and T denotes the total number of time steps considered.Denote g ( θ ) and g ( θ k | θ k − ), respectively, as the candidate distributions to sample θ from at time step 0 and at time step k ( ∀ k = 1 , , . . . , T ). Let ˜ θ k denote the setof accepted samples of θ at time k that will be passed on to the next time step( ∀ k ∈ { , , . . . , T − } ). Denote the set that contains all the accepted values of θ upto time k by ˜Θ( k ).We elaborate on the key steps of Algorithm 1 below. At time step 0, we draw n c candidate samples of θ from g ( θ ), denoted by { θ ∗ ,j } n c j =1 . We then evaluate theposterior distribution using the data at time 0 (i.e., D (0)) and the candidate particlesto obtain the unnormalized weights { w ∗ j , j = 1 , , . . . , n c } at time step 0, which aresubsequently normalized to the b w ∗ j . We then obtain n p posterior samples by selectingfrom { θ ∗ ,j , j = 1 , , . . . , n c } with the corresponding probabilities given by { b w ∗ j } n c j =1 ;denote this set of n p samples by ˜ θ = { ˜ θ , , ˜ θ , , . . . , ˜ θ ,n p } . This set of accepted valuesof θ will be passed on to time step 1. The sampling for time step k is similar to thatfor time step 0 but enhanced with sample augmentation. Specifically, at time step k ( ∀ k ∈ { , , . . . , T } ), using an appropriate candidate density g ( θ k | ˜ θ k − ,ℓ ) which isconditioned on each accepted sample from time step k −
1, we generate a batch of n b candidate samples in the neighborhood of ˜ θ k − ,ℓ ; denote the batch by { θ ∗ k,ℓ,j } n b j =1 for ℓ = 1 , , . . . , n p . We then evaluate the posterior distribution using the data up to time k (i.e., D ( k )), the accepted samples up to time step k − k − { w ∗ ℓ,j , j = 1 , , . . . , n b , ℓ = 1 , , . . . , n p } at time step k , which are subsequently normalized to the b w ∗ ℓ,j . We then obtain n p posterior samples by selecting from { θ ∗ k,ℓ,j , j = 1 , , . . . , n b , ℓ = 1 , , . . . , n p } with thecorresponding probabilities given by { b w ∗ ℓ,j , j = 1 , , . . . , n b , ℓ = 1 , , . . . , n p } ; denotethis set of n p accepted samples by ˜ θ k = { ˜ θ k, , ˜ θ k, , . . . , ˜ θ k,n p } . Then by the end of timestep k , we update the set of accepted values of θ to ˜Θ( k ) = ˜Θ( k − ∪ ˜ θ k .5 lgorithm 1 Enhanced sampling importance resampling algorithm Specify a prior distribution p ( θ ), and define D (0) = { S (0) , E (0) , I (0) , R (0) , D (0) } . Let D (0) = D (0) . Draw n c candidate samples from a candidate distribution g ( θ ), i.e., θ ∗ ,j ∼ g ( θ ), j = 1 , , . . . , n c , then D (1) is observed. Evaluate w ∗ j = p ( D (1) | D (0) , θ ∗ ,j ) p ( θ ∗ ,j ) /g ( θ ∗ ,j ), j = 1 , , . . . , n c . Obtain normalized weights b w ∗ j = w ∗ j / ( P n c ℓ =1 w ∗ ℓ ), j = 1 , , . . . , n c . Resample from { θ ∗ ,j } n c j =1 using the set of normalized weights, { b w ∗ j } n c j =1 , and obtaina set of n p samples, ˜ θ = { ˜ θ , , ˜ θ , , . . . , ˜ θ ,n p } . Let ˜Θ(0) = ˜ θ . for k = 1 , , . . . , T do Update the total dataset with new observations, D ( k ) = D ( k ) ∪ D ( k − D ( k ) = { S ( k ) , E ( k ) , I ( k ) , R ( k ) , D ( k ) } . Specify a prior distribution p k ( θ k ), then D ( k + 1) is observed. For the ℓ th sample in the collection ˜ θ k − from step k −
1, ˜ θ k − ,ℓ , draw n b samples from the candidate distribution g ( θ k | ˜ θ k − ,ℓ ), i.e., θ ∗ k,ℓ,j ∼ g ( θ k | ˜ θ k − ,ℓ ), j = 1 , , . . . , n b . Evaluate w ∗ ℓ,j = p ( D ( k + 1) | D ( k ) , θ ∗ k,ℓ,j , ˜Θ( k − p k ( θ ∗ k,ℓ,j ) /g ( θ ∗ k,ℓ,j | ˜ θ k − ,ℓ ), ℓ =1 , , . . . , n p , j = 1 , , . . . , n b . Normalize b w ∗ ℓ,j = w ∗ j,ℓ / (cid:16)P n p ℓ =1 P n b h =1 w ∗ h,ℓ (cid:17) , ℓ = 1 , , . . . , n p , j = 1 , , . . . , n b . Resample from { θ ∗ k,ℓ,j , ℓ = 1 , , . . . , n p ; j = 1 , , . . . , n b } using the set of nor-malized weights, { b w ∗ ℓ,j , ℓ = 1 , , . . . , n p ; j = 1 , , . . . , n b } , and obtain a set of n p samples, ˜ θ k = { ˜ θ k, , ˜ θ k, , . . . , ˜ θ k,n p } . Let ˜Θ( k ) = ˜Θ( k − ∪ ˜ θ k . end for5. Monitoring The presented method for estimating the SEIRD model parameters through time isquite responsive to changes in the system. Hence, these parameters can be monitoredto look for changes in the parameters which will be manifested in the data. Since thereare four parameters in the SEIRD model that need to be simultaneously monitored,the Multivariate Exponentially Weighted Moving Average (MEWMA) approach waschosen as an appropriate monitoring method. [16] developed a multivariate EWMA(MEWMA) control chart, which is an extension to the univariate EWMA.First the parameter samples were differenced using a single backward lag of one,∆ α ( t ) i = α ( t ) i − α ( t − i , ∆ γ ( t ) i = γ ( t ) i − γ ( t − i , ∆ β ( t ) i = β ( t ) i − β ( t − i and∆ η ( t ) i = η ( t ) i − η ( t − i . These form a vector (∆ α ( t ) i , ∆ γ ( t ) i , ∆ β ( t ) i , ∆ η ( t ) i ) T to bemonitored for significant deviations from zero, which would correspond to a significantchange in the set of parameters. The multivariate parameters are given by the mean:∆¯ θ ( t ) = 1 n p Tn p (∆ α ( t ) , ∆ γ ( t ) , ∆ β ( t ) , ∆ η ( t ))and variance: Cov (∆ θ ( t )) = 1 n p − ∆ α ( t ) T ∆ γ ( t ) T ∆ β ( t ) T ∆ η ( t ) T (∆ α ( t ) , ∆ γ ( t ) , ∆ β ( t ) , ∆ η ( t )) .
6n order to have a monitoring process that is not too sensitive, the MEWMA is em-ployed which is given by
M EW M A ( t ) = Λ∆¯ θ ( t ) + (1 − Λ) M EW M A ( t − , (4)with the moving covariance matrix: V ( t ) = Λ Cov (∆ θ ( t )) + (1 − Λ) V ( t − , where Λ is a smoothing coefficient that controls how much a new observation can influ-ence the overall mean. Lower values of Λ are more conservative in that new observationsdo not have much influence on the mean as higher values allow new observations tohave a greater influence on the mean. This is used to control false signals as any pro-cess being monitored will have some natural variation that may cause the observationto be signalled. Typical values of Λ include 0 . , . , . , .
25 and 0 . , , ,
0) indicating no shift in the process. In the multivariatecase, a test statistic based on Hotelling’s T can be formed as T ( t ) = M EW M A ( t ) T V ( t ) − M EW M A ( t ) . (5)When n p is large, T ( t ) in Equation (5) is approximately χ (4) which has a 0.95-quantile equal to 9.48. Hence any T ( t ) > .
48 would be deemed as a significant changein the parameter differences and thus a significant change in the SEIRD process. Notethat, in our case n p = 1 ,
000 and hence deemed large enough for the T approximation.
6. Data Analysis: State of Qatar
The World Health Organization, Johns Hopkins University and other agencies main-tain data sets on the daily number of confirmed infected cases, deaths, and recoveriesfor every country. In this work, we study the evolution of the pandemic in the Stateof Qatar. All data for Qatar were obtained from Johns Hopkins University and arefreely accessible via the Johns Hopkins COVID-19 GitHub repository [19]. The GitHubsite includes daily cumulative number of confirmed infections, cumulative number ofrecovered and cumulative number of deaths starting 22 January 2020.The goal of the data analysis is to demonstrate and assess the proposed modelingapproach and its use in monitoring the pandemic as the varying coefficient approachwill allow for the model parameters to adjust quickly to changes in the data generationprocess. In model (1), the Recovered and Death states are cumulative with no outgoingtransitions whereas the Infected state has transitions from Exposed and to Recoveredand Death states. Hence the data for confirmed infections are cumulative and includethe numbers from both the Recovered and Death states. As such, if CI ( t ) denotesthe confirmed infections at time t , then the number of Infected subjects at time t isdefined as I ( t ) = CI ( t ) − R ( t ) − D ( t ) . Figure 1.
Plot of Active Infections (a), Recovered (b) and Deaths (c) for the State of Qatar data.
From here on, the term “Active Infections” will be used to denote this derived variableversus the cumulative Infected provided in the data.Figure 1 shows the plots of daily Active Infections, Recovered and Deaths data forthe State of Qatar since 29 February 2020. The Active Infections start very low butencounter a large jump around day 12 due to increased testing. The Active Infectionsthen seem to plateau until day 30, after which there is an extreme growth in ActiveInfections. The Active Infections start to go down after day 90. The patterns in theRecovered and Deaths are similar and reflect the time of infection before recovery ordeath. Both graphs show a very slight increase in the number of recovered or deadsubjects until about day 90, after which a steady increase is noticed.
The State of Qatar is one of the countries heavily affected by the COVID-19 pandemic.Since its first infected case way back on February 29, 2020, Qatar has become one ofthe highest infected countries in the Middle East with the total number of confirmedcases standing at 148,258 as of January 26, 2021. The total number of deaths in Qatarso far stands at 248 cases, which is low relative to the total number of infected cases,which is an indication of the country’s highly effective healthcare system.Qatar prepared an excellent flexible plan for risk management, grounded on nationalrisk assessment, taking into account of the global risk assessment done by WHO andfocusing on reinforcing capacities to reduce or eliminate health risks from COVID-19. Along with the well-organized healthcare system, the country was very quick torespond to the global pandemic. The country implemented many preventive measuresvery early on in the pandemic, including border control for early detection of cases.This included, were but not limited to, installing thermal screening for passengerswho entered the country at Hamad International Airport and at seaports as early asJanuary 2020, with the first quarantine facilities opening on February 1 [1].On March 9, 2020 (day 10), Qatar closed all universities and schools and placed atravel ban on 15 countries: Bangladesh, China, Egypt, India, Iran, Iraq, Italy, Lebanon,Nepal, Pakistan, the Philippines, South Korea, Sri Lanka, Syria and Thailand. OnMarch 14, 2020 (day 15), Qatar expanded its travel ban to include three new coun-tries: Germany, Spain and France [11, 20]. The Ministry of Municipal and Environmenton March 21, 2020, closed all parks and public beaches to curb the spread of coron-avirus. On March 23, 2020 (day 24), the Ministry of Commerce and Industry decidedto temporarily close all restaurants, cafes, food outlets, and food trucks for the pub-lic. Also, the Ministry of Commerce and Industry decided to close all unnecessary8usinesses on March 27, 2020 (day 28) [11, 20].As the number of infected cases continued to rise, on April 8, 2020 (day 40), theMinistry of Public Health (MoPH) announced that Primary Health Care Coopera-tion will be designating two health centers, one in Umm-Salal and one in GharrafatAl-Rayyan, for screening, testing, and quarantining COVID-19 patients. MoPH alsoannounced a hotline for psychological aid on April 9, 2020 (day 41).These interventions made by the government change the dynamics of the pandemicand hence, need to be considered while setting up a a real-time monitoring systemof the infection, recovery and death rates. In the next subsection, we illustrate theproposed model as a data-driven forecasting model for use by stakeholders in theState of Qatar to monitor the COVID-19 pandemic.
For Qatar, the prior distribution specification is: α ( t ) ∼ Exp (2 / β ( t ) ∼ Exp (1 / γ ( t ) ∼ Exp (1 /
14) and η ( t ) ∼ Exp (1 / γ ( t ) the mean is 1 /
14 which corresponds to a 14-day infectionduration.The SEIRD model was run with the following initial values for Qatar: S (0) =2 , , E (0) = 3, I (0) = 1, R (0) = 0 and D (0) = 0. The SEIRD model and thesampling process given in Algorithm 1 were coded in MATLAB R2020a and were runon a PC with Intel Core i7-7700 CPU at 3.60GHz with 8GB of RAM. At each timestep the sampler was run with n c = 10 , n p = 1 ,
000 and n b = 10. Hence, each ofthe n p samples had a batch of n b samples generated in its neighborhood resulting in n p candidate particles at the beginning of the next time step. The model computationtime is about 60 minutes for T = 135 time steps. Note that, the number of individualsin each compartment in the model is much smaller due to the population size. Thismeans that many of the computations are faster especially when dealing with largefactorials associated with the Poisson distribution.Figure 2 shows the coefficient estimates across time, α ( t ) (panel a), β ( t ) (panelb), γ ( t ) (panel c) and η ( t ) (panel d) for the Qatar data set. Of particular interest isthe time frame from day 90 to day 95 in which the active infections I ( t ) exhibiteda large drop (see Figure 1). Notice that the distribution for α ( t ) becomes incrediblyconcentrated in this time frame as evidenced by the narrow credible intervals. This isfurther exhibited in β ( t ) and η ( t ). However, by examining γ ( t ) during this time frameone sees that, a large spike in the recovery rate with very narrow credible intervals aswell. And for γ ( t ) after about 80 days the recovery rate begins to increase dramaticallywith another spike around day 115.Figure 3 shows the model fitted to the data with 95% posterior predictive boundsfor Active Infections (panel a), Recovered (panel b) and Deaths (panel c). First note isthat all three models appear to fit the data extremely well based on visual inspection.Particularly notice that around day 90 when the dramatic drop in Active infectionsoccurs this can be contrasted with Figure 2.To assess the fit of the model a P seudo − R was calculated as: P seudo − R = 1 − P nt =1 h I ( t ) − b I ( t ) i + P nt =1 h R ( t ) − b R ( t ) i + P nt =1 h D ( t ) − b D ( t ) i P nt =1 h I ( t ) − I ( t ) i + P nt =1 h R ( t ) − R ( t ) i + P nt =1 h D ( t ) − D ( t ) i Figure 2.
Plots of α ( t ) (a), β ( t ) (b), γ ( t ) and η ( t ) across time with associated 95% credible intervals for theState of Qatar data. (a) (b) (c) Figure 3.
Plots of the data, fitted model with 95% posterior prediction bounds for Active Infections (a),Recovered (b) and Deaths (c) for the State of Qatar data. I ( t ), R ( t ) D ( t ) are the sample means of I , R and D , respectively, across time(and hence are not a function of time), b I ( t ), b R ( t ) and b D ( t ) are the medians of theposterior predictive distributions for I , R and D , respectively, at each time ( and henceare functions of time),Since uncertainty quantification is important, the proportion of observations thatfall into the predictive bands was calculated as follows: b P fit = P nt =1 I { I ( t ) ∈ b C I ( t ) } + P nt =1 I { R ( t ) ∈ b C R ( t ) } + P nt =1 I { D ( t ) ∈ b C D ( t ) } n where b C I ( t ) , b C R ( t ) and b C E ( t ) denote the 95% predictive intervals for I ( t ), R ( t ) and D ( t ), respectively and I {A} is an indicator function taking the value one if event A istrue.The Pseudo-R = 0 . b P fit = 0 . T ( t ) for the Qatar data with a control limit set at 9.48(red dashed line) and a smoothing parameter Λ = 0 .
2. Notice that until time 40 theprocess seems to be pretty stable as indicated by the T ( t ) being below the controllimit. After day 40 there are several time points signalling a change in the process ondays: 40, 44, 47, 64, 65, 67, 68, 69, 71, 77, 95, and 123. The early days (40,44,47) caneasily be seen to agree with the change in both infection rate α ( t ) and death rate η ( t )in Figure 2 (a) and (d). Notice in these plots the high volatility with spikes in α ( t )and a clear shift upwards in η ( t ) at the same times. During the days 64 to 77 thereappears to be a very large amount of volatility in the infection rate α ( t ) which is alsoevidenced in Figure 3 (a) for active infections. Notice how the active infections havelarge increases one day and smaller increases the next day. The method also picks upthe spike in infection rates as well as the dramatic increase in recovery rate γ ( t ) attime 95. At day 123, Figure 2 shows a spike in the rate of exposed to actively infected, β ( t ) (panel b) and another shift in recovery rate γ ( t ) (panel c).To study the sensitivity to the smoothing parameter in the MEWMA chart acrosssignals several other analyses were performed with Λ = 0 . , Λ = 0 . , Λ = 0 .
25, andΛ = 0 .
3, which are commonly chosen values for this parameter. When Λ = 0 . .
15 themonitoring process signaled at days: 40, 44, 47, 64, 68, 69, 71, 77 and 123. Notice thatfor both Λ = 0 . .
15 day 95 is not signaled, which upon inspection of almostall the charts there is a shift in the process. Hence, these parameter choices wouldbe considered too conservative in this case. When Λ = 0 .
25 the following days weresignalled: 40, 44, 47, 64, 65, 67, 68, 69, 71, 77, 95 and 123. This is the same resultas when Λ = 0 .
2. When Λ = 0 . α ( t ) across thistime frame. Overall, when less conservative values for Λ are chosen, the days signalledare quite reasonable. In could be argued that in a pandemic situation a more sensitivemonitoring process would be beneficial to public policy makers as it can signal when11 igure 4. Plot of the Hotelling’s T statistic through time. Horizontal line corresponds to the 95% controllimits an effective intervention has been introduced.
7. Discussion
This work provides a novel tool for monitoring and capturing changes in a pandemicevolution process via monitoring changes in parameters of mathematical epidemiolog-ical models, such as the Susceptible, Exposed, Infected, Recovered, Death (SEIRD)model using the Multivariate Exponentially Weighted Moving Average (MEWMA)process monitoring technique. A Bayesian approach is taken for the parameter esti-mation with a sampling algorithm that allows for both quick updating of the SEIRDmodel but also provides samples that can be monitored by the MEWMA regime.This sampling algorithm uses the notion of Sampling Importance Resampling, butaugments the particles at each step to avoid particle depletion. This quick updatingallows for the process monitoring scheme to “signal” quickly if there is a change in themodel parameters. The method is then used to monitor the evolution of the COVID-19pandemic in the State of Qatar.Despite the proliferation of forecasting models for the evolution of the COVID-19pandemic, their accuracy achieved can be compromised and comparisons can be com-plicated due to numerous factors, e.g., their construction methods, distinct healthcaresystems adopted by different countries/regions, different political decisions or policiesmade, distinct testing and reporting mechanisms [21]. Hence, using the forecasts givenby a particular forecasting model for critical decision making is challenging. The pro-posed approach takes a different perspective and enables decision-makers to work witha tailored SEIRD model, assess the effectiveness of the policies/decisions made, andadopt interventions and/or prevention strategies consistently over time.The State of Qatar example illustrates the proposed method’s ability to performdaily monitoring of a pandemic. The proposed model fits the data very well with a pseudo − R = 0 . R ,12hey would have a negligible effect on the fit. Furthermore, the model does not containcompartments for subjects who recovered without being confirmed infections. Sincethis is not observed, one can only speculate on the impact that additional data wouldhave on the model fit; however, it would be very small. As seen in Figures 2 and 3, theproposed method successfully picks up the day to day fluctuations in the pandemicevolution process in Qatar via the estimated time-varying model parameters. Notethat the pandemic’s overall state can also be monitored by tracking the T statisticover time (see Figure 4). For Qatar, the method signals the first change in the processaround day 40. This change can be attributed to several government interventions suchas closing parks and public beaches on day 24, closing all unnecessary businesses onday 28 and announcing two major health centers catered towards COVID-19 patientson day 40. The method also signals multiple days beyond day 40, all of which seemreasonable upon further inspection. Thus, the proposed method gives decision-makersthe ability to evaluate planned interventions as well as discover new changes to theprocess and respond accordingly. This method can also be extended for monitoring aprocess at the state/county level by incorporating a spatial covariance and using themixed model approach.Since the augmented sampling regime allows posterior samples to be saved fromthe previous day, updating is performed on a daily basis and only requires the newdata and the previous day’s samples. Thus the entire SEIRD model need not be fitfrom the beginning of the series. Furthermore, the MEWMA is quickly calculatedfrom the posterior samples and can quickly signal those managing the pandemic.Note that the method is not tied to the SEIRD model given in Equation (1), as theaugmented sampler and MEWMA monitoring protocol are generic. Our motivationhere has been using a system where the reproduction number fails to include all therelevant parameters. In systems where the reproduction number is dependent on allparameters, the reproduction number could be added as a dimension to the monitoringprotocol as well. In situations where the reproduction number is meaningful, thiscould be another dimension that could “signal” serious changes in long-term processoutcomes. References [1] Al Khal A., Al-Kaabi S., Checketts RJ.,
Qatar’s response to COVID-19 pandemic , HeartViews, 32 (2020), pp. 21-129.[2] Bagal, D. K., Rath, A., Barua, A., & Patnaik, D.,
Estimating the parameters of susceptible-infected-recovered model of COVID-19 cases in India during lockdown periods , Chaos, Soli-tons & Fractals, 140 (2020), pp. 110-154.[3] Basu, D., Salvatore, M., Ray, D., Kleinsasser, M., Purkayastha, S., Bhattacharyya,R., & Mukherjee, B.,
A comprehensive public health evaluation of lockdown as a non-pharmaceutical intervention on COVID-19 spread in India: National trends masking statelevel variations , (2020). medRxiv.[4] Chinazzi M. Davis, J.T., Ajelli, M., Gioannini, C., Litvinova, M., Merler, S., y Piontti,A.P., Mu, K., Rossi, L., Sun, K. and Viboud, C.,
The effect of travel restrictions on thespread of the 2019 novel coronavirus (COVID-19) outbreak , Science (2020).[5] Chowell G., Brauer F.,
The basic reproduction number of infectious diseases: Computa-tion and estimation using compartmental epidemic models , Mathematical and statisticalestimation approaches in epidemiology (2009), Springer.[6] Clancy, D., & O’Neill, P. D.,
Bayesian estimation of the basic reproduction number instochastic epidemic models , Bayesian Analysis, 3(4), 737-757 (2008).
7] Deb, S., & Majumdar, M.,
A time series method to analyze incidence pattern and estimatereproduction number of COVID-19 , arXiv preprint arXiv:2003.10655 (2020).[8] Fisher D., Heymann D.,
Q&A: The novel coronavirus outbreak causing COVID-19 , BMCMedicine. doi: 10.1186/s12916-020-01533-w (2020).[9] Ghanam, Ryad, Edward Boone, and Abdel-Salam Abdel-Salam,
COVID-19: SEIRDModel for Qatar COVID-19 Outbreak , Letters in Biomathematics, June (2020). https://lettersinbiomath.journals.publicknowledgeproject.org/index.php/lib/article/view/323 [10] Giuliani D. Dickson, M.M., Espa, G. and Santi, F.,
Modelling and predicting thespread of Coronavirus (COVID-19) infection in NUTS-3 Italian regions , arXiv preprintarXiv:2003.06664 (2020).[11] Hamad Medical Corporation. (2020). Major Risks to Business Continuity. Doha, Qatar.[12] Han, Z., Li, T., & You, J.,
These Unprecedented Times: The Dynamic Pattern Of COVID-19 Deaths Around The World , (2020). arXiv preprint arXiv:2011.02824.[13] Jewell, C. P., Kypraios, T., Neal, P., & Roberts, G. O.,
Bayesian analysis for emerginginfectious diseases , Bayesian analysis, 4(3), 465-496 (2009).[14] Kyrychko, Y.N., Blyuss, K.B. & Brovchenko, I.,
Mathematical modelling of the dy-namics and containment of COVID-19 in Ukraine , Scientific Reports 10 (2020), 19662.https://doi.org/10.1038/s41598-020-76710-1[15] Loli Piccolomini E, Zama F.,
Monitoring Italian COVID-19spread by a forced SEIRD model , PLoS ONE 15(2020), e0237417. https://doi.org/10.1371/journal.pone.0237417 [16] Lowry, C.A., Woodall, W.H., Champ, C.W. and Rigdon, S.E.,
A multivariate exponen-tially weighted moving average control chart , Technometrics, 34 (1992), pp. 46-53.[17] May, Robert M.; Anderson, Roy M.,
Infectious diseases of humans: dynamics and control ,Oxford: Oxford University Press, 1991. ISBN 0-19-854040-X.[18] McIntosh K.,
Coronavirus disease 2019 (COVID-19), Up-To-Date , Bulletin-Association of Canadian Map Libraries and Archives (ACMLA), 164 (2020), pp. 47-51.[20] Ministry of Public Health. (2019). Qatar National Preparedness and Response Plan forCommunicable Diseases. Doha, Qatar.[21] Nikolopoulos, K., Punia, S., Sch¨afers, A., Tsinopoulos, C. and Vasilakis, C.,
Forecastingand planning during a pandemic: COVID-19 growth rates, supply chain disruptions, andgovernmental decisions , European Journal of Operational Research 290 (2021), pp. 99-115.[22] Rezabakhsh A. Ala, A. and Khodaei, S.H.,
Novel Coronavirus (COVID-19): A NewEmerging Pandemic Threat , Journal of Research in Clinical Medicine 8(1), pp. 5-6 (2020)[23] Roy, A., & Karmakar, S.,
Bayesian semiparametric time varying model for count data tostudy the spread of the COVID-19 cases , (2020). arXiv preprint arXiv:2004.02281.[24] Vynnycky, E.; White, R. G.,
An Introduction to Infectious Disease Modelling , Oxford:Oxford University Press, eds. (2010). ISBN 978-0-19-856576-5.[25] Wang, L., Zhou, Y., He, J., Zhu, B., Wang, F., Tang, L., ... & Song, P. X. (2020),
Anepidemiological forecast model and software assessing interventions on the COVID-19 epi-demic in China , Journal of Data Science, 18(3), 409-432 (2020).[26] World Health Organization,
Novel Coronavirus (2019-nCoV) situation reports .[27] Wu F. Zhao, S., Yu B. Chen, Y.M. Wang, W. Song, Z.G. Hu, Y. Tao, Z.W. Tian, J.H.Pei, Y.Y. and Yuan, M.L.,
A new coronavirus associated with human respiratory diseasein China , Nature (2020), pp. 265-269., Nature (2020), pp. 265-269.