以貝氏馬可夫蒙地卡羅方法估計隨機微分方程呼吸道傳染病模式

Stochastic Differential Equation Model for Respiratory Infectious Disease Using Bayesian Markov Chain Monte Carlo (MCMC) Method

楊雅閔(Ya-Min Yang)

指導教授 ： 陳秀熙

國立臺灣大學/公共衛生學院/流行病學與預防醫學研究所/博士(2017年)

https://doi.org/10.6342/NTU201700245

摘要

研究背景
過去呼吸道傳染病的數學模式大多以決定性微分方程為基礎，由於此類傳染病牽涉到為數眾多的易感族群，因此與其相關的參數如基礎再生數相對穩定。然而，近來新興傳染病(如中東呼吸症候群冠狀病毒感染症)及現行流行性感冒的監測則卻因病毒基因變異頻繁而增加其複雜性，也因此，需要利用隨機性微分方程模式處理此類傳染病在族群中的異質性，如個案共病症有無及家戶間的異質性等。
就方法學觀點而論，利用隨機微分方程模型過去往往遭遇求解困難。本論文利用貝氏馬可夫蒙地卡羅方法以電腦模擬方式得到相關參數的事後分佈，除了可以推導基礎再生數，估計傳染病絕跡的機率之外，尚可進行局部和全球的漸進穩定性和平衡點的評估分析，並據此發展傳染病動態過程的電腦輔助應用工具。

研究目的
本論文研究目的在應用隨機微分方程模式於傳染病SIR或其他相關模式，發展貝氏馬可夫蒙地卡羅電腦模擬程式用以估計呼吸道傳染病之基礎再生數及相關統計量如傳染病絕跡機率以作為傳染病監視之用，本論文以新興中東呼吸道症候群及現存的流行性感冒為例。特定研究目的如下：
(1) 發展貝氏馬可夫蒙地卡羅程式多階層多區間模式(包括SIR模式及出生-死亡模式型)以評估共病症等變項對感染者復原率的影響；
(2) 發展貝氏馬可夫蒙地卡羅程式用於微分方程多區間隨機過程，以中東呼吸道症候群為例估計基礎再生數及傳染病絕跡機率，並與傳統最大概似方法(非微分方程為基礎)做比較；
(3) 利用(2)的結果評估局部和全球的漸進穩定性和平衡點；
(4) 利用(2)的結果發展傳染病動態過程的電腦輔助，並考量部份免疫力及疫苗的影響。

材料與方法
1. 實證資料
(1) 流行性感冒
本研究流行性感冒資料來自健保申報資料，個案層級資料包括性別、年齡、是否施打疫苗及感冒就診日期，此外，亦收集家戶資料用以建構多階層隨機模型。
(2) 中東呼吸症候群冠狀病毒感染症(MERS)
本研究MERS資料取自於Rambaut 依據世界衛生組織公告的各國感染症資訊及文獻資料整理而得的線上資料庫。
2. 模式建構及統計分析
　　本研究發展一系列貝氏馬可夫蒙地卡羅電腦語法，在SAS環境下完成以下模型建構：
(1) 多階層二階段馬可夫模型
(A) 建立以貝氏多階層之離散狀態連續時間馬可夫模型為基礎之二階段多階層模型，評估影響多國家流行的中東呼吸道症候群死亡風險之因素。
(B) 建立以貝氏架構為基礎之貝塔-二項分佈迴歸模型運用於流行性感冒家戶資料，評估個人之罹病危險與基礎再生數(R0)。
(2) 本研究發展出一系列微分方程隨機過程傳染病數學模式，如SIR模式及生死過程等，並應用於中東呼吸道症候群之流行資料，用以推估其基礎再生數，並透過負二項分佈的性質推算傳染病絕跡的機率。
(3) 運用(2)估計所得的參數，進行中東呼吸道症候群局部和全球的漸進穩定性和平衡點的評估。
(4) 在微軟Excel的環境之下發展傳染病動態過程的電腦輔助，並考量部份免疫力及疫苗的影響。

結果
1. 發展多階層二階段馬可夫模型運用於傳染病實證資料
(1) 中東呼吸症候群的流行
MERS 致死率為32.1% (95%CI: 29.9-34.3%)，死亡風險在第一周平均增加率為每日13%，其後2周則為3%。在考慮共病以及年齡之影響並考量國家間的異質性後，具有共病者其死亡風險為未有共病者之4倍(aHR: 3.6，95% CI: 2.46-5.33)。
(2) 流行性感冒
罹病危險隨傳染世代增加而增加，與第一世代相較，第二到第四世代之罹病危險分別為35%、52%以及約兩倍。施打疫苗者之罹病危險較低(OR: 0.96，95%CI: 0.9-1.04)。R¬0於第一世代之1.76增加至至第四世代之3.64。在同一世代中，施打疫苗者相較於未施打者其R¬0 較低。以第一世代為例，施打疫苗者與未施打者其R0分別為1.70 (95% CI: 1.35-2.04) 與1.76 (95% CI: 1.40-2.12)。
2. 利用隨機過程傳染病數學模式估計中東呼吸症候群的流行
(1) 以非微分方程為基礎的貝氏蒙地卡羅馬可夫鏈模型估計中東呼吸症候群的流行，結果發現在不同國家，不同城市的不同時間點，若疾病流行已經結束，則所得的基礎再生數約略相同，皆在1左右。然而若是疾病仍在流行，基礎再生數則會大於1。吉達(2014/3-6)、利雅德(1, 2014/3-6)、利雅德(2, 2015/6-9)和南韓(2014/5-8)的基礎再生數分別為1.003、2.14、1.04和1.04。第60天中東呼吸症候群可滅絕的預測機率分別為72%、29%、71%及75%。
(2) 使用以微分方程為基礎的貝氏蒙地卡羅馬可夫鏈模型解三間隔的微分方程時，基礎再生數在在吉達(2014/3-6)、利雅德(1, 2014/3-6)、利雅德(2, 2015/6-9)和南韓(2014/5-8)分別為2.3、10.8、2.4及2.6。第60天中東呼吸症候群可滅絕的預測機率分別為41%、9%、36%及38%。
3. 局部和全球的漸進穩定性和平衡點的評估分析
不論利用非微分方程為基礎或以隨機微分方程為基礎的貝氏蒙地卡羅估計在四個區域均顯示平衡點僅針對無病個案(且估計結果不穩定)，因此中東呼吸症候群的流行不會被滅絕，且可能變成地方性流行或者流行性傳染病。
4. 傳染病動態過程的電腦輔助應用工具設計
電腦輔助應用工具之發展及呈現以中東呼吸道症候群及流行性感冒為例，藉由傳染病動力曲線呈現兩種傳染病在不同傳染動力模型下，模擬族群在可感染、感染及復原間隨時間的動態變化，基礎有效再生數亦可產生以評估疾病傳播快慢。疾病在傳染期間之可感染、感染及復原人數亦可透過此應用工具預測各階段人數。

結論
本論文發展以貝氏馬可夫蒙地卡羅方法為基礎的傳染病模式，包括多階層模式及隨機微分方程模型，並將之應用於流行性感冒及中東呼吸道症候群實證資料，推導其基礎再生數及傳染病絕跡機率，最後則根據上述研究架構開發傳染病動態過程電腦輔助工具以作為疾病監測之依據。

關鍵字

隨機過程 ； 傳染病 ； 中東呼吸症候群冠狀病毒感染症 ； 基礎再生數 ； 局部和全球的漸進穩定 ； 平衡分析

並列摘要

Background
Mathematical modelling for elucidating respiratory infectious disease used to rely on deterministic differential equation model as preventive strategies (such as vaccination, isolation, and quarantine) for containing such a kind of infectious disease often involves a large susceptible population that renders the parameter of interest i.e. basic reproduction number (R0) rather stable. Nonetheless, the surveillance of the emerging infectious disease such as MERS-CoV and the existing disease like influenza that is highly variable in mutation suggests that the deterministic one may need to be supplanted by the stochastic one because heterogeneity observed in these types of infectious diseases. For example, the sequale of susceptible with co-morbid diseases would be different from that of those without. The cause of heterogeneity may be also attributed to small proportion of susceptible clustering in the same setting (such as household), which makes the variation change from cluster to cluster.
From the viewpoint of methodology, the use of stochastic equation is faced with intractable solution or estimation of a string of differential equations. The recently proposed Bayesian Markov Chain Monte Carlo (MCMC) provides an opportunity to solve this issue but this has been never addressed before. The development of feasible and efficient computer algorithms for these ODEs play an important role in the estimation of basic reproduction number, the evaluation of the local and the global stability with respect to the equilibrium process, and the development of computer-aided tool of dynamic epidemic process of infection.

Objectives
The objectives of the thesis are therefore to apply stochastic ODE to the SIR or related models and to develop computer algorithms with Bayesian MCMC underpinning so as to to estimate the basic reproduction number and other interests of statistics such as the probability of extinction for the surveillance of respiratory infectious diseases with two illustrations, the emerging MERS and the commonly seen influenza. The specific aims are to
(1) develop Bayesian MCMC computer algorithms for multi-level compartment model to assess the effect of covariate of interest (such as co-morbidity) on the infected state and the recovery state using the SIR model or birth-and-death process;
(2) develop Bayesian MCMC computer algorithms for ODE-based compartment stochastic process in comparison with conventional maximum likelihood (non-ODE based) approach to data of MERS for estimating infection rate, the discharge rate, and death rate so as to derive and the probability of extinction through negative binomial transformation;
(3) evaluate local and global asymptotic stability and equilibrium based on (2);
(4) develop a computer-aided tool for dynamic epidemic process of infection based on (2) under consideration partial immunity and vaccination.

Materials and Methods
1. Empirical Data
(1) Influenza epidemic
Data on whether subjects developed influenza-like symptoms were derived from the National Health Insurance database. Information on the characteristics of the subjects, such as age, sex, vaccination status, and date of diagnosis were collected. Household information was also collected for the construction of multi-level stochastic model. .
(2) MERS-CoV
Data on reported cases of MERS cases since April, 2012 were derived from a web-based line list maintained by Rambaut with reference to the published literatures and the announcement of the disease information system of each country and WHO.
2. The Development of Bayesian MCMC Computer Algorithms
This thesis developed a series of Bayesian MCMC computer algorithms with SAS environment for the following proposed model.
(1) Multilevel Two-state Process
(A) We proposed a Bayesian multilevel model with discrete-state and continuous-time Markov model to assess the risk of death of MERS infective cases across continent countries.
(B) We developed a Beta-binomial model with Bayesian underpinning to assess the effect of individual level factors on the risk of contracting influenza and to derive the individualized basic reproduction number (R0) considering randomness of household.
(2) Stochastic processes with ordinary differential equation like the SIR model or birth-and-death process were propose and applied to epidemic of MERS CoV so as to derive and the probability of extinction through negative binomial transformation.
3. Using the estimated parameters obtained above, local and global asymptotic stability and equilibrium analysis of MERS-CoV were assessed.
4. Finally, a computer-aided tool for dynamic epidemic process of infection under consideration of partial immunity and vaccination by Microsoft Excel 2016 was developed.

Results
1. Bayesian MCMC Multilevel Compartment Model
(1) MERS-CoV
The fatality rate of MERS cases was 32.1% (95% confidence interval: 29.9%, 34.3%). Notably, the incremental change of daily death rate was most prominent during the first week since disease onset with an average increase of 13%, but then stabilized in the remaining two weeks when it only increased 3% on average. After adjusting for age and making allowance for the heterogeneity across countries with random-effect, MERS patients with comorbidity had around 4 times the risk for fatal infection than those without (adjusted hazard ratio of 3.60 (95% confidence interval: 2.46, 5.33)).
(2) Influenza
The risk of developing into influenza cases was increased with generation, compared with the first generation, the increased proportions for the second to fourth generation were 35%, 52% and around two-fold, respectively. Those who were vaccinated had a lower risk of developing into influenza cases (Odds ratio (OR): 0.96, 95% CI: 0.90-1.04)).
　　The R0 increased from around 1.76 in the first generation to around 3.64 for the fourth generation. The vaccinated group had a lower R0 than those without for the same generation. For example, the R0 estimates during the first generation were estimated as 1.70 (95% CI: 1.35-2.04) for the vaccinated group and 1.76 (95% CI: 1.40-2.12) for the unvaccinated group, respectively.
2. Bayesian MCMC Stochastic Process Applied to MERS-CoV Epidemic
(1) Using non-ODE-based MCMC method, the estimated basic reproduction number estimates were around unity and similar in different countries/cities in different time periods.. The reproduction number in Jeddah (2014/3-6), Riyadh (1, 2014/3-6), Riyadh (2, 2015/6-9), and South Korea (2014/5-8) were 1.003, 2.14, 1.04, and 1.04, respectively. The corresponding probabilities of extinction at day 60 were 72%, 29%, 71%, and 75%.
(2) The estimated reproduction numbers with ODE-based MCMC method in Jeddah (2014/3-6), Riyadh (1, 2014/3-6), Riyadh (2, 2015/6-9), and South Korea (2014/5-8) were 2.3, 10.8, 2.4, and 2.6, respectively. The corresponding probabilities of extinction at day 60 were 41%, 9%, 36%, and 38%.
3. Evaluation of local and global asymptotic stability and equilibrium analysis
As all the estimated reproduction number were greater than unity regardless of non-ODE- or ODE-based methods and the equilibrium point was only for disease free but unstable, the MERS-CoV infection would not be eradicated and may become endemic or epidemic.
4. Computer-aided tool of dynamic epidemic process of infection
The developed application was demonstrated with illustrations of MERS and influenza A. Evolution of susceptible, infective and recovery by different dynamic epidemic models was characterized by kinetic epidemic curves. Basic and effective reproduction number (R0) was calculated to provide an quantitative measure for assessing the spread of disease. The predicted number of different compartments in dynamic process of infection during the epidemic period could be obtained.

Conclusions
The novel computer algorithms for stochastic ordinary differential equation (ODE) model and multi-level model with Bayesian MCMC underpinning has been developed and applied to influenza and MERS-CoV to demonstrate how to derive basic reproduction number and the extinction probability in a stochastic manner so as to develop computer-aided tools for the surveillance of dynamic epidemics of respiratory infectious disease.

並列關鍵字

stochastic process ； infectious disease ； MERS-CoV ； basic reproduction number ； local and global stability ； equilibrium analysis

參考文獻

AL-TAWFIQ, J. A., HINEDI, K., GHANDOUR, J., KHAIRALLA, H., MUSLEH, S., UJAYLI, A. & MEMISH, Z. A. (2014). Middle East respiratory syndrome coronavirus: a case-control study of hospitalized patients. Clin Infect Dis 59, 160-5.

ALSAHAFI, A. J. & CHENG, A. C. (2016). The epidemiology of Middle East respiratory syndrome coronavirus in the Kingdom of Saudi Arabia, 2012–2015. International Journal of Infectious Diseases 45, 1-4.

ANDERSON, R. M. (1991). Discussion: the Kermack-McKendrick epidemic threshold theorem. Bull Math Biol 53, 3-32.

ANDERSON, R. M. & MAY, R. M. (1985). Age-related changes in the rate of disease transmission: implications for the design of vaccination programmes. J Hyg (Lond) 94, 365-436.

ANDERSSON, M., EKDAHL, K., MOLSTAD, S., PERSSON, K., HANSSON, H. B. & GIESECKE, J. (2005). Modelling the spread of penicillin-resistant Streptococcus pneumoniae in day-care and evaluation of intervention. Stat Med 24, 3593-607.