王 婷 曾 平，2△

霍亂傳染數(shù)據(jù)貝葉斯ZIP分析的后驗(yàn)預(yù)測(cè)評(píng)價(jià)*

王 婷1曾 平1，2△

1926年英國(guó)流行病學(xué)家麥克德里克(McKendrick)應(yīng)用數(shù)學(xué)模型研究了1906年印度孟買一個(gè)村莊流行性霍亂傳染的數(shù)據(jù)[1]。為了描述霍亂傳播過(guò)程，麥克德里克首先考慮了簡(jiǎn)單的Poisson分布，在發(fā)現(xiàn)Poisson分布擬合效果很差后，麥克德里克又考慮了零截尾Poisson分布(zero-truncated Poisson)，Irwin采用EM算法改進(jìn)了麥克德里克的方法[2]。但是無(wú)論哪種分布，麥克德里克都沒能很好處理霍亂數(shù)據(jù)中存在的大量沒有患者的家庭，即霍亂感染者個(gè)數(shù)為0。麥克德里克也注意到這個(gè)問題，并給出了一個(gè)十分合理的解釋：霍亂病菌經(jīng)水源傳播，村莊中有許多水井，也許只有某些水井被病菌污染，而一些水井沒有被污染，是干凈清潔的。因此，飲用被病菌污染井水的家庭中有人患病，而飲用干凈井水的村民自然不會(huì)患病，這就揭示了為什么數(shù)據(jù)中存在大量沒有霍亂患者的家庭。根據(jù)麥克德里克的理由，Meng認(rèn)為可采用了零膨脹Poisson模型(zero inflated Poisson，ZIP)來(lái)分析霍亂傳染數(shù)據(jù)[3]。一個(gè)很自然的問題是，ZIP模型對(duì)該數(shù)據(jù)的擬合優(yōu)度如何？如何做出評(píng)價(jià)？為此，本文主要研究基于后驗(yàn)預(yù)測(cè)分布的貝葉斯ZIP模型評(píng)價(jià)，并通過(guò)該數(shù)據(jù)展示相應(yīng)的方法。

貝葉斯ZIP模型

霍亂傳染數(shù)據(jù)描述了孟買一個(gè)村莊223戶家庭流行性霍亂傳染情況[2]，見表1。表中x表示被感染家庭中霍亂患者人數(shù)，nx表示患者人數(shù)為x的家庭個(gè)數(shù)。數(shù)據(jù)顯示大約75%的家庭沒有霍亂患者，大約14%的家庭有1個(gè)患者，7%的家庭有2個(gè)患者，患者人數(shù)超過(guò)3個(gè)的家庭數(shù)不足4%?；魜y傳染數(shù)據(jù)的ZIP模型可表示如下[4-6]：

(1)

p表示家庭中成員不可能感染霍亂的概率，I為指示函數(shù)，當(dāng)x=0取值為1，否則為0。ZIP分布實(shí)際是退化到0的分布(概率為p)和Poisson分布(概率為1-p)的混合分布。

(2)

∝表示兩邊相差一個(gè)常數(shù)因子。采用基于MCMC(Markov Chain Monte Carlo)的隨機(jī)游走M(jìn)etropolis模擬算法[7]，通過(guò)構(gòu)造平穩(wěn)分布為后驗(yàn)分布的Markov鏈來(lái)得到參數(shù)的后驗(yàn)樣本(pt,μt)，Monte Carlo樣本量為105。圖1-2給出了參數(shù)p和μ的直方圖和等高線圖[8]，表2為后驗(yàn)統(tǒng)計(jì)量。隨機(jī)游走M(jìn)etropolis模擬在SAS9.2的MCMC程序下完成[7]，后續(xù)統(tǒng)計(jì)分析在R2.11.1下完成[8-9]。

*:數(shù)據(jù)來(lái)源于《現(xiàn)代數(shù)學(xué)手冊(cè)-隨機(jī)數(shù)學(xué)卷》。貝葉斯ZIP模型擬合時(shí)選擇參數(shù)的中位數(shù)，分別為p=0.5930和μ=0.9470。

后驗(yàn)預(yù)測(cè)模型評(píng)價(jià)

模型后驗(yàn)評(píng)價(jià)是建立在后驗(yàn)預(yù)測(cè)分布基礎(chǔ)上的。后驗(yàn)預(yù)測(cè)分布(posterior predictive distribution)是指在觀察到數(shù)據(jù)X之后未來(lái)可能觀察到數(shù)據(jù)X*的分布[10-11]。稱后驗(yàn)是因?yàn)閿?shù)據(jù)X*的分布是建立在已觀察到數(shù)據(jù)X之上的條件分布，稱預(yù)測(cè)是因?yàn)樵摲植际轻槍?duì)尚沒有觀察到的數(shù)據(jù)而言。X*的后驗(yàn)預(yù)測(cè)分布可表示為：

f(X*|X)=∫f(X*|X,p,μ)f(p,μ|X)dpdμ, =∫f(X*|p,u)f(p,μ|X)dpdμ,(3)

上式中第二個(gè)等式成立是基于f(X*|X,p,u)=f(X*|p,μ)的事實(shí)，也即是X*的后驗(yàn)預(yù)測(cè)分布是參數(shù)p,μ的條件分布，X只通過(guò)參數(shù)p,μ來(lái)影響X*的分布，和已經(jīng)觀察到的數(shù)據(jù)X無(wú)關(guān)，在已知p,μ時(shí)X*和X條件獨(dú)立。從后驗(yàn)預(yù)測(cè)分布中產(chǎn)生的數(shù)據(jù)又稱為重復(fù)數(shù)據(jù)(replicated data)，本文中X指223戶家庭霍亂發(fā)病人數(shù)，X*指預(yù)測(cè)霍亂患者人數(shù)。后驗(yàn)預(yù)測(cè)分布量化了在已經(jīng)觀察到數(shù)據(jù)X后在未來(lái)能夠再次觀察到X的可能性，因此不但可以用來(lái)預(yù)測(cè)未知數(shù)據(jù)，還能進(jìn)行模型檢查評(píng)價(jià)。

從后驗(yàn)預(yù)測(cè)分布隨機(jī)抽取樣本X*，然后將這些樣本和觀察到的數(shù)據(jù)X比較，如果模型對(duì)數(shù)據(jù)擬合得好，那么重復(fù)數(shù)據(jù)應(yīng)該看上去和數(shù)據(jù)X一致，兩者之間的明顯差別表明選擇模型的偏離?？梢酝ㄟ^(guò)定義一個(gè)檢驗(yàn)量T來(lái)量化這種偏離程度[11-12]，T(X*)表示應(yīng)用預(yù)測(cè)數(shù)據(jù)得到的檢驗(yàn)量，為一個(gè)隨機(jī)變量，T(X)表示應(yīng)用觀察數(shù)據(jù)得到的檢驗(yàn)量，為一個(gè)確定的數(shù)值。在霍亂傳染數(shù)據(jù)的ZIP模型中，合適的檢驗(yàn)量是數(shù)據(jù)中0的家庭數(shù)T(X*)=n0。后驗(yàn)檢驗(yàn)的策略是，如果ZIP模型擬合良好，那么從預(yù)測(cè)分布f(X*|X)產(chǎn)生的0的家庭數(shù)T(X*)=n0和T(X)=n0應(yīng)該接近，T(X*)和T(X)之間系統(tǒng)的差別提示ZIP模型擬合不足。本文中T(X)=168，圖3給出了霍亂傳染數(shù)據(jù)和9個(gè)預(yù)測(cè)數(shù)據(jù)的頻數(shù)分布圖，其中預(yù)測(cè)數(shù)據(jù)由公式(3)隨機(jī)產(chǎn)生。圖3顯示預(yù)測(cè)數(shù)據(jù)分布和實(shí)際觀察數(shù)據(jù)圖形相近。

霍亂傳染數(shù)據(jù)ZIP模型的后驗(yàn)?zāi)Ｐ蜋z驗(yàn)具體過(guò)程如下[11-12]：

(1)從后驗(yàn)分布f(p,μ|X)抽取(pt,μt),t=1，…,104,這一過(guò)程已經(jīng)由隨機(jī)游走M(jìn)etropolis模擬完成；

(3)計(jì)算后驗(yàn)預(yù)測(cè)樣本X*中患者為0的家庭數(shù)T(X*)；

(4)對(duì)所有104組參數(shù)(pt,μt)重復(fù)(1)～(3)步，則產(chǎn)生T(X*)t,t=1,…,104；

如果ZIP能擬合霍亂數(shù)據(jù)，則p不應(yīng)該出現(xiàn)極端值，例如很大(大于0.95)或很小(小于0.05)。T(X*)的直方圖見圖4，參考線為168的位置。T(X*)≥T(X)和T(X*)≤T(X)的比例分別是0.507和0.537，兩者之和不為1是因?yàn)椴坏仁街卸即嬖诘忍?hào)的原因。這表明如果霍亂患者傳染數(shù)據(jù)服從ZIP分布，那么觀察到高達(dá)75%的家庭沒有患者的可能很大，這種情況不大可能是偶然發(fā)生的。最后的結(jié)果顯示T(X*)在135～196之間，P25和P75分別為162和174，均數(shù)和中位數(shù)分別為167.5和168，標(biāo)準(zhǔn)差為9.04。特殊地，如果用p=0.593和μ=0.947預(yù)測(cè)霍亂發(fā)病人數(shù)，0的預(yù)測(cè)值為167，概率為75.01%，和實(shí)際168十分一致，除此之外，其他預(yù)測(cè)值和實(shí)際觀察到的發(fā)病人數(shù)也很接近。這些數(shù)值表示ZIP模型很好地?cái)M合了實(shí)際數(shù)據(jù)，能夠有效地處理數(shù)據(jù)中存在的大量0記錄的情況。

如果選擇Poisson模型，先驗(yàn)選擇共軛gamma分布，指定α=β=10-3，貝葉斯Poisson模型T(X*)的分布如圖5。T(X*)介于112～184之間，P25和P75分別為145和158，均數(shù)和中位數(shù)分別為151.7和152，標(biāo)準(zhǔn)差為9.38。T(X*)≥T(X)和T(X*)≤T(X)的比例分別是0.044和0.966，這表示從Poisson模型中不大可能觀察到高達(dá)75%的零記錄，Poisson模型對(duì)數(shù)據(jù)擬合不夠，至少不能處理數(shù)據(jù)中零過(guò)多的問題。

討 論

Poisson分布和其他模型對(duì)霍亂傳染數(shù)據(jù)擬合效果差的一個(gè)主要原因在于不能很好地處理那些沒有患者的家庭。由于霍亂病毒只污染了村里部分水井，那么飲用健康井水的村民無(wú)疑是不會(huì)患病的，即便是飲用霍亂污染的井水，村民也不一定表現(xiàn)出明顯的患病癥狀，正是這兩種原因使得很多家庭沒有霍亂患者，其中前者是Poisson分布不能解釋的0的來(lái)源。我們采用ZIP模型對(duì)這個(gè)混合總體的不同人群分別建立模型，給前者指定為退化到0點(diǎn)的分布，后者指定為Poisson分布，ZIP模型就是這兩個(gè)分布的混合，準(zhǔn)確地捕捉到了數(shù)據(jù)結(jié)構(gòu)所呈現(xiàn)的信息。貝葉斯ZIP模型結(jié)果顯示，預(yù)測(cè)的家庭患病人數(shù)和實(shí)際觀察數(shù)據(jù)十分接近，尤其是準(zhǔn)確地預(yù)測(cè)了沒有患者的家庭數(shù)。

為了評(píng)價(jià)麥克德里克問題ZIP模型的擬合優(yōu)度，

我們采用后驗(yàn)預(yù)測(cè)分布對(duì)ZIP模型做出評(píng)價(jià)。后驗(yàn)預(yù)測(cè)檢查需要定義一個(gè)合適的檢驗(yàn)量T，T不但可以是數(shù)據(jù)的函數(shù)T(X)，如本文中所選擇的檢驗(yàn)量，這和經(jīng)典統(tǒng)計(jì)假設(shè)檢驗(yàn)統(tǒng)計(jì)量類似，還可以是數(shù)據(jù)和參數(shù)的函數(shù)T(X,θ)，后驗(yàn)預(yù)測(cè)分布重復(fù)數(shù)據(jù)的T(X*)和實(shí)際數(shù)據(jù)的T(X)之間的差異反映了模型和數(shù)據(jù)之間的吻合程度。后驗(yàn)預(yù)測(cè)檢查的原理在于，如果實(shí)際數(shù)據(jù)違背了模型的重要假設(shè)，那么T(X*)相對(duì)T(X)就會(huì)表現(xiàn)出具有統(tǒng)計(jì)學(xué)意義的差別來(lái)。如果得到的后驗(yàn)預(yù)測(cè)p值接近0或1，暗示實(shí)際數(shù)據(jù)具有極端的檢驗(yàn)量T，而選擇的模型很可能是不合適的。在后驗(yàn)評(píng)價(jià)中采用重復(fù)數(shù)據(jù)的0的比例作為檢驗(yàn)量，這也是選擇ZIP模型的主要理由。因此，如果對(duì)ZIP模型是恰當(dāng)?shù)?，那么由后?yàn)預(yù)測(cè)分布所產(chǎn)生的0就應(yīng)該和實(shí)際觀察到的0的比例一致。后驗(yàn)檢查顯示，ZIP模型能夠很好地?cái)M合實(shí)際數(shù)據(jù)，而Poisson模型不能很好地處理數(shù)據(jù)中存在的0過(guò)多問題。除此之外，差別檢驗(yàn)量還可選擇χ2值或變異系數(shù)。

[1]Http://en.wikipedia.org/wiki/Anderson_Gray_McKendrick.

[2]徐利治.現(xiàn)代數(shù)學(xué)手冊(cè)-隨機(jī)數(shù)學(xué)卷.武漢:華中科技大學(xué)出版社,1999：357-390.

[3]Meng XL.The EM Algorithm and Medical Studies:A Historical Link.Statistical Methods s in Medical Research,1997,6(1):3-23.

[4]曾平,劉桂芬,曹紅艷.零膨脹模型在心肌缺血節(jié)段數(shù)影響因素研究中的應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2008,25(5):464-466.

[5]曾平.零過(guò)多計(jì)數(shù)資料回歸模型及其醫(yī)學(xué)應(yīng)用.碩士學(xué)位論文,太原:山西醫(yī)科大學(xué),2009.

[6]Lambert D.Zero-inflated Poisson Regression with an Application to Defects in Manufacturing.Technimetrics,1992,34(1):1-14.

[7]SAS Institute Inc.Preliminary Capabilities for Bayesian Analysis in SAS/STAT Software.Cary,NC,USA,2006.

[8]Albert J.Bayesian Computation with R.Second Edition.New York:Springer,2009.

[9]R Development Core Team.R:A language and environment for statistical computing.R Foundation for Statistical Computing,Vienna,Austria,2014.URL http://www.R-project.org.

[10]Ntzoufras I.Bayesian Modeling Using WinBUGS.New York:John Wiley & Sons,2009.

[11]Gelman A,Carlin JB,Stern HS,et al.Bayesian Data Analysis.Second Edition.London:Chapman & Hall,2004.

[12]Gelman A,Meng XL,Stern HS.Posterior Predictive Assessment of Model Fitness via Realized Discrepancies.Statistica Sinica,1996,6(4):733-807.

[13]Meng XL.Posterior predictive p-values.Annals of Statistics,1994,22(3):1142-1160.

(責(zé)任編輯：鄧 妍)

*國(guó)家自然科學(xué)基金項(xiàng)目(81402765)；國(guó)家統(tǒng)計(jì)局全國(guó)統(tǒng)計(jì)科學(xué)研究項(xiàng)目(2014464)；江蘇省教育廳高校哲學(xué)社會(huì)科學(xué)研究基金項(xiàng)目(2013SJB790059,2013SJD790032)

1.徐州醫(yī)學(xué)院公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)教研室(221004)

2.南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系

△ 通信作者：曾平，E-mail：zpstat@xzmc.edu.cn