烏魯木齊市新冠肺炎疫情數(shù)據(jù)的貝葉斯分析

古麗斯坦·庫爾班尼牙孜田茂再

(1.新疆財經(jīng)大學(xué) 統(tǒng)計與數(shù)據(jù)科學(xué)學(xué)院，新疆烏魯木齊 830012;2.中國人民大學(xué) 應(yīng)用統(tǒng)計科學(xué)研究中心，中國人民大學(xué) 統(tǒng)計學(xué)院，北京 100972)

§1 引言

新型冠狀病毒肺炎(Corona Virus Disease 2019，COVID-19)，簡稱“新冠肺炎”，是一種由2019新型冠狀病毒引起的急性呼吸道傳染病.新冠肺炎出現(xiàn)以后持續(xù)蔓延至許多國家，嚴重威脅著人類的生命健康，因此引起了世界各國的廣泛關(guān)注與高度重視.在全球抗擊新冠肺炎的過程中，眾多學(xué)者對疫情傳播和發(fā)展趨勢進行了大量研究.Guan和Hu[1]等收集截止2020年1月29日中國30個省，自治區(qū)，直轄市552家醫(yī)院的1099例新冠肺炎確診病例，對新冠肺炎的潛伏期，死亡率進行估計.劉勛，孟秋雨[2]等收集湖北省截至2020年2月10日新冠肺炎疫情數(shù)據(jù)，利用OpenGeoda軟件對湖北省新冠肺炎發(fā)病的空間聚集性特征進行分析.周濤[3]等對我國COVID-19基本再生數(shù)值進行預(yù)測.白堯等[4]利用SEIAR模型對陜西省疫情風(fēng)險進行了評估.

國家衛(wèi)健委組織每天發(fā)布新冠肺炎每日新增確診病例與每日新增無癥狀感染者病例，這為新冠肺炎的建模和分析提供了良好的基礎(chǔ).它們在描述疫情的傳播特征，探討疫情傳播的影響因素，在數(shù)據(jù)驅(qū)動下預(yù)測疫情的可能發(fā)展趨勢，評估防疫措施的有效性等方面起著重要作用.因此從理論層面探究這兩組數(shù)據(jù)的分布特征，通過基于數(shù)學(xué)，統(tǒng)計學(xué)以及計算機科學(xué)的完整有效的分析方法來對數(shù)據(jù)進行建模，從樣本入手挖掘出隱藏的，有價值的信息在新冠疫情分析中具有很重要的意義，并且對新冠疫情分析方法的重要補充.

目前國內(nèi)外新冠疫情的有關(guān)研究主要集中在疫情發(fā)展趨勢的預(yù)測[5-9]，疫情對經(jīng)濟[10]，旅游[11]等方面的影響與沖擊等層面.利用貝葉斯方法探索疫情數(shù)據(jù)的分布特征方面的研究甚少.基于此，本文根據(jù)新疆烏魯木齊市2020年7月15日至8月15日新冠肺炎疫情每日新增確診病例與2020年7月15日至8月16日每日新增無癥狀感染者病例的實際數(shù)據(jù)，分別基于泊松分布，變點泊松分布與過度發(fā)散泊松分布建立統(tǒng)計模型，并利用貝葉斯方法對模型進行統(tǒng)計推斷，旨在從理論層面揭示這兩組數(shù)據(jù)的內(nèi)在統(tǒng)計特征和統(tǒng)計規(guī)律，能為抗擊疫情提供一定的理論依據(jù).

統(tǒng)計推斷中最常用的方法是古典統(tǒng)計與貝葉斯統(tǒng)計.其中古典統(tǒng)計是基于大樣本推斷總體，把未知參數(shù)看作一個未知的固定量，缺點是過于依賴樣本信息，如果樣本有偏或樣本稀疏，則會使參數(shù)估計出現(xiàn)偏差.而貝葉斯統(tǒng)計是把未知參數(shù)也看作是隨機變量，并充分利用先驗信息，更加符合人類認知自然的過程.特別是對于小樣本數(shù)據(jù)而言，基于大樣本推斷總體的傳統(tǒng)統(tǒng)計模型受到更多挑戰(zhàn).而貝葉斯層次模型通過引入先驗信息，充分考慮各種不確定性，可在一定程度上克服小樣本缺陷[12-15].從數(shù)據(jù)自身的特征來看，新冠肺炎每日新增確診病例與每日新增無癥狀感染者病例每天只有一個觀測值，因此不滿足大樣本性質(zhì);再者每日新增確診病例與每日新增無癥狀感染者病例也有可能不滿足獨立同分布假定.綜上所述，直接利用古典統(tǒng)計方法，這兩組數(shù)據(jù)可能會突破大樣本與獨立同分布的兩個前提條件.而利用貝葉斯方法在理論上可能會更加嚴密可靠.因此為了突出貝葉斯方法在疫情數(shù)據(jù)分析中的重要作用，本文通過貝葉斯方法和分層貝葉斯方法對烏魯木齊市每日新增確診病例與每日新增無癥狀感染者病例進行分析.

本文的安排如下:§2分別基于泊松分布，變點泊松分布和過度發(fā)散泊松分布對烏魯木齊市每日新增確診病例與每日新增無癥狀感染者病例進行統(tǒng)計建模，并利用貝葉斯方法對模型進行統(tǒng)計推斷;§3利用DIC值對三種模型的擬合優(yōu)度進行評價;§4是本文的結(jié)論.

§2 烏魯木齊市新冠肺炎數(shù)據(jù)的貝葉斯分析

本文研究的數(shù)據(jù)來源為2020年7月15日至2020年8月16日全國健康委員會公布的烏魯木齊市每日新增確診病例與每日新增無癥狀感染者病例數(shù)據(jù).每日新增確診病例的樣本區(qū)間為7月5至8月15日，樣本個數(shù)為n32;每日新增無癥狀感染者病例的樣本區(qū)間為7月5至8月16日，樣本個數(shù)為n33.樣本期間烏魯木齊市新冠肺炎每日新增確診病例與無癥狀感染者病例變化趨勢如下圖1所示.下文對這兩組數(shù)據(jù)的分布特征進行統(tǒng)計建模.

圖1 每日新增確診病例與無癥狀感染者病例的變化趨勢圖

2.1 基于泊松分布模型的貝葉斯分析

泊松分布是概率論中常用的一種離散型概率分布，適合于描述單位時間內(nèi)隨機事件發(fā)生的次數(shù)，如機器出現(xiàn)故障的次數(shù);汽車站內(nèi)候客人數(shù);自然災(zāi)難發(fā)生的次數(shù)等.由于新冠肺炎每日新增確診病例與每日新增無癥狀感染者病例是計數(shù)數(shù)據(jù)，且為離散隨機變量，所以從數(shù)據(jù)的這個特點出發(fā)，先利用泊松分布進行建模.

2.1.1 基于泊松分布模型的每日新增確診病例貝葉斯分析

假設(shè)烏魯木齊市總?cè)丝跒镹(N355.2萬人)，每日新增確診病例發(fā)病率為p1，每日新增確診病例Y1，···，Yn是獨立同分布的離散隨機變量，且Yi服從二項分布b(N，p1).因為N很大，p1很小，根據(jù)泊松定理[16]，隨機變量Yi的分布可以用參數(shù)為λNp1的泊松分布來近似，即

其中λ代表每日平均新增確診病例.

為了獲取參數(shù)λ的后驗分布，需要先給定相應(yīng)的先驗分布.根據(jù)國內(nèi)外研究經(jīng)驗，λ的先驗分布假定為伽馬分布，即π(λ)Ga(α，β).樣本的聯(lián)合分布函數(shù)和π(λ)相乘可以得到聯(lián)合后驗分布

因此每日平均新增確診病例λ的條件后驗分布為

其中Ga(·)代表伽馬分布.基于參數(shù)λ的條件后驗分布，本文采用MCMC方法來對后驗分布進行抽樣.MCMC算法的基本思想是:通過從后驗分布中獲取隨機樣本，建立一條平穩(wěn)分布為后驗分布的馬爾科夫鏈，然后根據(jù)這些樣本對未知參數(shù)進行統(tǒng)計推斷.MCMC算法中常用的抽樣方法有Metropolis-Hasting算法，Gibbs算法，獨立抽樣算法等.其中最簡單，應(yīng)用最廣泛的MCMC方法是Gibbs算法.當(dāng)參數(shù)的后驗分布為標準分布時，Gibbs算法直接從參數(shù)的后驗分布中抽取樣本[17].根據(jù)(4)式可知，參數(shù)λ的后驗分布為常見的伽馬分布，所以可以用Gibbs抽樣方法直接對驗分布進行抽樣.本文利用Winbugs軟件進行操作[18].假定參數(shù)λ的的初始值為10，超參數(shù)α4，β2，馬爾科夫迭代鏈長度為10000.表1給出了參數(shù)估計的最終結(jié)果，圖2是MCMC方法數(shù)據(jù)模擬得到的參數(shù)λ的軌跡圖，核密度圖和自相關(guān)函數(shù)圖.

表1 每日平均新增確診病例的貝葉斯估計

從表1可以看出，每日平均新增確診病例的貝葉斯估計25.21，跟實際值25.93比較接近，參數(shù)估計的MC誤差趨向于零，說明用模型的擬合效果較好.由圖2可以看出，抽樣值基本都均值附近波動，核密度估計圖近似正態(tài)，自相關(guān)函數(shù)值趨近于零，說明產(chǎn)生的馬爾科夫鏈呈現(xiàn)收斂狀態(tài)，根據(jù)λ與p1的關(guān)系計算出來的樣本期間新增確診病例平均發(fā)病率p1的估計值為2.2436×10-4，這跟實際發(fā)病率p12.2375×10-4的比較接近.總之，模型來擬合每日新增確診病例比較合理.

圖2 參數(shù)λ的迭代軌跡，核密度與自相關(guān)函數(shù)圖

2.1.2 基于泊松分布模型的每日新增無癥狀感染者病例貝葉斯分析

同理，假設(shè)樣本期間烏魯木齊市新增無癥狀感染者發(fā)病率為p1′，平均新增無癥狀感染者病例為ηNp1′.每日新增無癥狀感染者病例Z1，···，Zm服從參數(shù)為η的泊松分布，即

其中z(z1，···，zm)是每日新增確診病例的樣本觀測值.設(shè)η的先驗分布為伽馬分布Ga(α1，β1).得到的每日平均新增無癥狀感染者病例η條件后驗分布為

利用Gibbs算法對η的后驗分布進行抽樣，初始值取為η5，超參數(shù)設(shè)定為α2，β2.表2給出了η的最終參數(shù)估計結(jié)果，圖3是MCMC方法數(shù)據(jù)模擬下得到參數(shù)η的迭代軌跡圖，核密度圖與自相關(guān)函數(shù)圖.

表2 η的貝葉斯估計結(jié)果

圖3 參數(shù)η 的迭代軌跡，核密度與自相關(guān)函數(shù)圖

從表2的結(jié)果可以看出，每日平均新增無癥狀感染者病例η的貝葉斯估計為11.24，估計結(jié)果跟樣本均值11.51非常接近，參數(shù)估計的MC誤差趨向于零，說明模型的擬合效果很好.由圖3可以看出，參數(shù)η的迭代軌跡基本穩(wěn)定，核密度估計圖近似正態(tài)，自相關(guān)函數(shù)值趨近于零，因此可以認為迭代已經(jīng)收斂.根據(jù)η與p1′的關(guān)系得到的每日新增無癥狀感染者平均發(fā)病率的估計為1.0250×10-4，這跟無癥狀感染者實際發(fā)病率p1′1.0839×10-4很接近，因此利用貝葉斯方法估計模型參數(shù)是有效的.

2.2 基于變點泊松分布模型的貝葉斯分析

在上文的建模過程中假定每日新增確診病例和每日新增無癥狀感染者病例在整個樣本期間是同分布的，但實際數(shù)據(jù)有可能不滿足這個假定.為了更深入地探究數(shù)據(jù)的分布特征，本節(jié)利用變點泊松分布對數(shù)據(jù)進行建模.在變點模型中先把樣本按照出現(xiàn)時間的先后順序排列，如果在某個未知的時刻開始樣本的數(shù)字特征和分布特征突變，那這個未知的時刻稱為變點.變點檢測分為對分布參數(shù)的變點檢測和分布本身的變點檢測[19].本文基于分布參數(shù)的變點檢測方法，只有存在一個變點的情形下分別對兩組數(shù)據(jù)進行建模.

2.2.1 基于變點泊松分布模型的每日新增確診病例貝葉斯分析

假設(shè)τ是每日新增確診病例的一個變點，新增確診病例的發(fā)病率為p2，每日新增確診病例Yi服從參數(shù)為λNp2的泊松分布.模型的第二層對λ進行建模.假設(shè)，當(dāng)i1，···，τ時，λexp(β1);當(dāng)iτ+1，···，n時λexp(β1+β2).把這兩個模型結(jié)合可以寫成

其中δ()是示性函數(shù)，若()里的表達式為非負其值為1，否則為零.參數(shù)β1與β2的先驗分布設(shè)定為獨立先驗，并假定為π(βk)～ N(0，0.001)(k1，2).假定變點發(fā)生在樣本出現(xiàn)的時間間隔內(nèi)，變點參數(shù)τ的先驗分布為U(1，n)，n32.利用MCMC方法對未知參數(shù)β1，β2和τ進行貝葉斯推斷，進行10000次迭代，參數(shù)初始值的設(shè)定為τ5，(β1，β2)(0，0)，表3列出了得到的參數(shù)估計結(jié)果.圖4，圖5給出了參數(shù)β1，β2和τ的迭代軌跡圖核密度估計圖.

表3 每日新增確診病例變點泊松模型參數(shù)的估計結(jié)果

圖4 參數(shù)β1，β2與τ的迭代軌跡圖

圖5 參數(shù)β1，β2與τ的核密度圖

2.2.2 基于變點泊松分布模型的每日新增無癥狀感染者病例貝葉斯分析

同理，假設(shè)τ1是每日新增無癥狀感染者病例的一個變點，每日新增無癥狀感染者病例zj～poisson(η)(j1，···，m).假定當(dāng)j1，···，τ1時，ηexp(γ1).當(dāng)jτ1+1，···，m時ηexp(γ1+γ2).因此有

其中δ()是示性函數(shù)，若()里的表達式為非負其值為1，否則為零.假定γ1與γ2的先驗分布互相獨立，并假定為π(γk)～ N(0，0.001)(k1，2)，變點參數(shù)τ1的先驗分布為U(1，m)，m33.利用MCMC方法對未知參數(shù)進行貝葉斯推斷，進行10000次迭代.參數(shù)初始值的設(shè)定為(γ1，γ2)(0，0)，τ15.最終得到的參數(shù)估計結(jié)果如下表4所示.

表4 每日新增無癥狀病例變點泊松模型的估計結(jié)果

圖6 參數(shù)γ1，γ2 與τ1的迭代軌跡圖

從表4的結(jié)果可以發(fā)現(xiàn)，的置信區(qū)間不包括零，說明是顯著不為零，意味著15.5是每日新增無癥狀感染者病例的一個變點.的符號是負的，可以認為開始每日新增無癥狀感染者病例開始減少.根據(jù)以上結(jié)果計算可得，7月15日至29日，每日平均新增無癥狀感染者的貝葉斯估計為e2.85117.305，發(fā)病率的估計為7.4083×10-5，這跟樣本均值η17.333 和實際發(fā)病率為p′7.4162×10-5非常接近.7月30日到8月16日每日平均新增無癥狀感染者的貝葉斯估計為e(2.851-0.958)6.6392，實際值為η6.6667，發(fā)病率的估計為3.4122×10-5，實際發(fā)病率為p′3.4229×10-5，可見模型的估計效果很好.參數(shù)估計的迭代軌跡圖和核密度圖分別在如下圖7，圖8所示.由圖7，圖8可以看出，當(dāng)?shù)螖?shù)到達一定數(shù)量后參數(shù)γ1，γ2與τ1的抽樣值基本都在均值附近波動，可以認為產(chǎn)生的鏈具有較好的收斂性.從核密度估計來看γ1，γ2的核密度基本正態(tài)，變點參數(shù)τ1的核密度大致可以看成均勻分布，說明參數(shù)的先驗分布均為共軛先驗分布.

圖7 參數(shù)γ1，γ2與τ1的核密度圖

2.3 基于過度發(fā)散泊松分布模型的貝葉斯分析

由實際樣本數(shù)據(jù)計算可得，樣本期間烏魯木齊市新冠肺炎每日新增確診病例的均值為E(y)25.937，方差為Var(y)717.546;每日新增無癥狀感染者病例的均值為E(Z)11.515，方差Var(Z)97.695，顯然這兩組數(shù)據(jù)的均值與方差不相等，不滿足泊松分布的均值與方差相等的特征，因此這兩組數(shù)據(jù)可以認為過度發(fā)散數(shù)據(jù)[20].本節(jié)基于過度發(fā)散泊松分布對這兩組數(shù)據(jù)進行建模.

2.3.1 基于過度發(fā)散泊松分布模型的每日新增確診病例貝葉斯分析

假設(shè)每日新增確診病例的發(fā)病率為pi，每日新增確診病例Yi(i1，···，n)服從參數(shù)為λiNpi的泊松分布，具體為

分層模型的第二層對λi進行建模，假定λi服從形狀參數(shù)為α，尺度參數(shù)為β的伽馬分布.即

其中α，β為超參數(shù)，分層模型的第三層假定α的先驗分布為指數(shù)分布，β的先驗分布為伽馬分布.具體而言

其中A，B，C為確定的常數(shù).令y(y1，···，yn)，λ(λ1，···，λn)，把參數(shù)先驗分布乘以似然函數(shù)，可以得到聯(lián)合后驗分布為:

從以上的聯(lián)合后驗密度函數(shù)可以推出，λi，β和α的條件后驗分布分別為

可見參數(shù)λ與β的后驗分布均有標準的分布形式，直接可以用Gibbs算法進行抽樣，但參數(shù)α的后驗分布沒有標準的分布形式，需要采用MH 算法進行抽樣，因此本文采用Gibbs抽樣和MH算法相結(jié)合的MCMC方法[21-22]估計模型的未知參數(shù)θ(λ，α，β).MH算法抽樣過程中所需要的提議分布選取為0到1上的均勻分布U(0，1).通過抽樣產(chǎn)生的序列(θ(0)，θ(1))，···，θ(k)是一個馬爾科夫鏈，它的平穩(wěn)分布收斂于聯(lián)合后驗密度p(θ|y)[20].抽樣次數(shù)設(shè)定為N10000，未知參數(shù)的初始值設(shè)定為λi(0)5(i1，···，32)，(α(0)，β(0))(2，1)，超參數(shù)設(shè)定為A2，B2，C1.表5給出了最終得到的參數(shù)估計結(jié)果.因為篇幅有限，下面的圖8，圖9只給出了λ的最大分量λ16與超參數(shù)α，β對應(yīng)的迭代軌跡圖和核密度估計圖.

表5的結(jié)果表明，每日平均新增確診病例λi(i1，···，n)的估計值與實際值非常接近，標準差也比較小，參數(shù)估計的MC誤差趨向于零，利用貝葉斯方法對模型估計是有效的.由圖8，圖9可以發(fā)現(xiàn)，當(dāng)?shù)螖?shù)到達一定數(shù)量后參數(shù)α，β和λ16的迭代軌跡基本穩(wěn)定，可以認為所產(chǎn)生的馬爾科夫鏈收斂.對于核密度而言，λi(i1，···，32)的核密度近似為正態(tài)分布，參數(shù)α和β的核密度近似可以看成伽馬分布，這說明參數(shù)β的先驗分布是共軛先驗分布.

圖8 參數(shù)λ16，α與β的迭代軌跡圖

圖9 參數(shù)λ16，α與β的核密度圖

表5 每日平均新增確診病例的貝葉斯估計結(jié)果

2.3.2 基于過度發(fā)散泊松分布模型的每日新增無癥狀感染者病例貝葉斯分析

假設(shè)烏魯木齊市7月15日至8月16日每日新增無癥狀感染者病例發(fā)病率為pj(j1，···，m)，每日新增無癥狀感染者病例zj(j1，...，m)服從參數(shù)為ηjNpj的泊松分布，取ηi的先驗分布為Ga(α2，β2)，超參數(shù)α2的先驗分布為exp(A1)，β2的先驗分布為Ga(B1，C1)，A1，B1，C1為給定的常數(shù).記η(η1，...，ηm)，通過跟上文一樣的思路與步驟，利用MCMC方法對未知參數(shù)(η，α2，β2)進行估計，總共進行N10000 次迭代，未知參數(shù)的初始值設(shè)定為ηj(0)8(j1，···，m)，(α2(0)，β2(0))(2，1)，常數(shù)A13，B12，C12.最終得到的參數(shù)估計結(jié)果如下表6所示.因為篇幅有限，圖10，圖11中只給出了η的最大分量η11與超參數(shù)α2，β2對應(yīng)的迭代軌跡圖，核密度圖.表6中的參數(shù)估計結(jié)果表明，每日新增無癥狀感染者ηj(j1，···，m)與發(fā)病率pj的估計值跟實際值很接近;估計量的MC 誤差趨近于零，因此可以認為模型對實際數(shù)據(jù)的擬合效果很好.由圖10，圖11可以看出，當(dāng)?shù)螖?shù)到達一定數(shù)量后參數(shù)η11與α2與β2的迭代軌跡基本穩(wěn)定，所產(chǎn)生的的馬爾科夫鏈收斂.參數(shù)η11的核密度近似正態(tài)，α2和β2核密度近似為伽馬分布，β2的先驗分布是共軛先驗分布.

圖10 參數(shù)η11，α2與β2的迭代軌跡圖

圖11 參數(shù)η11，α2與β2的核密度圖

§3 模型的評價

由上文的結(jié)果可以看出，對于每日新增確診病例和每日新增無癥狀感染者病例而言，三個模型均有較好的擬合效果，為了選擇其中的最優(yōu)模型，利用Spiegelhalter等[23]提出的DIC準則來對三種模型的擬合優(yōu)度進行評價.DIC(Spiegelhalter[23])是Akaike信息標準(AIC)在貝葉斯統(tǒng)計方面的推廣.它的具體形式為

表7 各個模型的DIC值計算結(jié)果

由表7的結(jié)果可知，對于烏魯木齊新冠肺炎每日新增確診病例與每日新增無癥狀感染者病例而言，基于過度發(fā)散泊松分布模型的DIC值最小，其次為變點泊松模型，泊松模型的DIC值最大，這說明過度發(fā)散泊松分布模型能夠為最準確地刻畫出這兩組數(shù)據(jù)的分布特征，因此對這兩組數(shù)據(jù)利用過度發(fā)散泊松分布來建模是最合理的.

§4 結(jié)論

每日新增確診病例和每日新增無癥狀感染者病例對新冠肺炎傳播過程的描述，分析和探討疫情傳播的影響因素方面重要作用，并且為新冠肺炎的建模和分析提供了良好的基礎(chǔ).因此本文從烏魯木齊市新冠肺炎疫情每日新增確診病例與無癥狀感染者病例的實際數(shù)據(jù)出發(fā)，以探索數(shù)據(jù)的分布特征為目的，分別基于泊松分布，變點泊松分布，過度發(fā)散泊松分布對數(shù)據(jù)進行統(tǒng)計建模，利用貝葉斯方法估計模型中的未知參數(shù)，并利用DIC值對模型的擬合優(yōu)度進行評價.結(jié)果表明，三種建模方法對兩組數(shù)據(jù)都具有較好的適用性，相比之下，基于過度發(fā)散泊松分布的模型在擬合優(yōu)度方面均優(yōu)于其它兩個模型，因此對這兩組數(shù)據(jù)利用過度發(fā)散泊松分布來進行建模是最為合理的，其次為基于變點泊松分布的模型，泊松分布模型的擬合優(yōu)度相對較差.