關(guān)于傳染性負二項分布發(fā)病率的貝葉斯推斷

戴 成，陳博鈺，田茂再

(中國人民大學 統(tǒng)計學院，北京 100872)

0 引言

在流行病學研究領(lǐng)域中，當所研究疾病的發(fā)病率很小時，傳統(tǒng)方法為采用負二項分布即逆抽樣方法對發(fā)病率進行研究，并假定發(fā)病率固定不變。但當該種疾病具有傳染性時，可能由于多種因素而導致發(fā)病率產(chǎn)生變動，如新近發(fā)生的甲型H1N1流感的在墨西哥爆發(fā)，鄰近的美國居民獲知該信息便可能采取相應(yīng)措施來預防流感病毒，會使得相應(yīng)的染病率低于在墨西哥的發(fā)病率。在這種情況下，若采用傳統(tǒng)方法假設(shè)發(fā)病率不變，則不能有效處理 (Bradlow,Hardie,and Fader,(2002);Philip J and Eric T(2002))。

為了有效地對發(fā)病率進行統(tǒng)計推斷，本文考慮將發(fā)病率考慮為一個變動的變量，采用分層貝葉斯模型進行處理。在貝葉斯理論中，選擇參數(shù)的先驗分布起到舉足輕重的作用，一般經(jīng)驗為令先驗分布與后驗分布達到共軛 (Howard Raiffa and Robert Schlaifer(1961)和 George Alfred Barnard(2005))，較常見的是當正態(tài)分布的均值和方差均未知時，選擇正態(tài)逆Gamma分布作為其先驗概率分布(Gelman,Carlin,Stern和Rubin(2003))。

此外，beta分布由于能夠為伯努利分布、二項分布、負二項分布和幾何分布等提供有效的先驗概率分布，因此在貝葉斯推斷過程中被大量的使用。盡管如此，beta分布仍存在一定的弊端，如當尾部較薄時，分布會迅速上升的性質(zhì)將導致一些擬合問題 (詳見 James J和Melvin R(1984))。

為了克服標準beta分布的弊端，本文考慮將廣義第二類beta分布作為發(fā)病率的先驗概率分布引入分層貝葉斯模型，來對發(fā)病率進行統(tǒng)計推斷和計算貝葉斯可信區(qū)間。

1 建立模型

此處將X定義為在第r次試驗成功前失敗的試驗次數(shù)，p為實驗成功的概率，0＜1＜p。顯然變量X服從負二項分布NB(r,p)，其概率質(zhì)量函數(shù)為

通常假定p為常量，但在傳染性疾病的研究領(lǐng)域中，這種假設(shè)與真實情況相悖，發(fā)病率p會隨著疾病的蔓延而發(fā)生變動，此處將其納入貝葉斯思想的研究，即假定p的先驗概率密度分布來推斷其后驗概率分布情況。這里，我們考慮引入廣義第二類beta分布作為p的先驗分布。

1.1 廣義第二類beta分布

我們將p考慮成一個在0～1之間變動的變量，對于負二項分布的概率 p，一般將其先驗分布選擇為beta分布。beta分布存在一定的缺點，因此曾有學者將指數(shù)函數(shù)等引入beta分布而將其廣義化以優(yōu)化其分布性質(zhì)。

Saralees Nadarajah和 Samuel Kotz在超幾何分布函數(shù)的基礎(chǔ)上提出了一種新的廣義第二類beta分布，包含三個參數(shù)(a,b,γ)，其概率密度函數(shù)為

其中2F1(1-γ,a;a+b;x)為高斯超幾何函數(shù)，其函數(shù)表達式為

2F1γ＞0；(z)k 為遞增階乘，

Saralees Nadarajah和Samuel Kotz(2003)給出了廣義第二類beta分布的矩估計函數(shù):

1.2 負二項分布X的邊際密度函數(shù)

在p的先驗概率分布為廣義第二類beta分布(2)式的情況下，由負二項分布(1)式產(chǎn)生的X的邊際密度函數(shù)見(4)式(詳細推導略)：

其中3F2(1-γ,a,a+b+r;a+b,a+b+r+x+1;1)為廣義超幾何函數(shù)， 對于 Re(a+b+γ)＞0絕對收斂 (參見 Earl D.Rainville(1971))。

此處我們可以通過高斯超幾何函數(shù)的特殊性質(zhì)來獲得關(guān)于 (3)式的一些相對特殊的情形 (參見Prudnikov et al.(1990)和 Gradshteyn&Ryzhik(2000))，具體如下：

（1）如果 a+b+γ=1，則(3)式簡化為

利用高斯加和定理，我們可以獲得

其中(c-a-b)＞0,c≠0,-1,-2,-3,…；由(6)推知

(2)如果 γ=1 或 a=0，則有

由此 (3)式簡化為

(3)還有其他的一些特殊情形如 (a)a+b-1;;(b)γ=0;(c)b=1,…。

當(a,b,γ)在條件(a)下時，(2) 式簡化為

若a+b-1為整數(shù)，則有

由(11)式可見，用其來簡化(3)式過于復雜(同理可推得(b),(c),(d)條件下的情形)。

1.3 結(jié)構(gòu)性質(zhì)

上面我們推導出負二項分布第r次成功前試驗失敗次數(shù)X的邊際密度函數(shù)，以下我們將探討它的一些性質(zhì)，包括均值、方差以及矩生成函數(shù)。

1.3.1 均值與方差

（1）均值：由雙期望定理可知，

均值=E（X）=E（E（X|p））

其中E(X|p)為給定p下隨機變量X的條件期望，已知X|p服從負二項分布(1)式，其均值和方差分別為

(2)方差:根據(jù)(16),我們可由X的條件方差公式得到X的方差

由(14)和(16)式可得

1.3.2 矩生成函數(shù)

由于超幾何函數(shù)3F2(1-γ,a,a+b+r;a+b,a+b+r+x+1;1)為無窮序列，其總和難于計算導致(3)式的的導數(shù)無法直接求得。(4)式是由條件負二項分布的積分而得，已知X|p的矩生函數(shù)為

由于p服從廣義第二類beta分布，由此推得X的矩生成函數(shù)為(推導過程略)

2 分層貝葉斯模型

由貝葉斯定理我們可以得到所關(guān)心的變量——發(fā)病率p的后驗概率分布，然后可根據(jù)其后驗概率分布對其進行統(tǒng)計推斷，構(gòu)造bayes可信區(qū)間。在上面我們已經(jīng)得到X的邊際密度函數(shù)，在給定（x,a,b,γ）的情況下，p的后驗概率密度分布為(推導過程略)

此處關(guān)于超參數(shù)(a,b,γ)，我們假設(shè)其相互獨立且均服從廣義均勻分布，即

則有 f(a,b,γ)=cacbcγ，由此可得

由此可得

由于(23)式中包含了超幾何分布的無窮級數(shù)求和，使得計算過程過于繁瑣，因此考慮采用計算機密集計算MCMC方法來對其進行計算機實現(xiàn)。通過R軟件編程，我們可以極大地簡化運算過程，得到（a,b,γ）的極大似然估計，由此來對發(fā)病概率p進行統(tǒng)計推斷。

我們所建立的分層貝葉斯模型具體為：

第一層：超參數(shù)（a,b,γ）服從廣義均勻分布(見(21)式)；

第二層：發(fā)病率p服從廣義第二類beta分布，參數(shù)為（a,b,γ）(見(2)式)；

第三層：逆抽樣得到的X服從負二項分布(見(1)式)。

由此我們可以得到關(guān)于超參數(shù)（a,b,γ）的聯(lián)合后驗概率分布函數(shù)(見(22)和(23)式)，在給定樣本的情況下，通過簡單的數(shù)值計算我們可以得到（a,b,γ）的估計，從而得到對p做統(tǒng)計推斷的結(jié)果。

在構(gòu)造p的置信區(qū)間方面，假設(shè)1-α為置信水平，求出滿足(24)的 pu和 pl,

從而可求得 p 的 1-α 置信區(qū)間[pl,pu]，其中，對(a,b,γ)的極大似然估計及p的置信區(qū)間的算法都是采用Gibbs抽樣下MCMC完成。

3 結(jié)論

本文考慮了當有傳染性的疾病蔓延時，發(fā)病率會相應(yīng)產(chǎn)生變動，使得傳統(tǒng)方法即假定發(fā)病率不變的負二項分布模型不能夠有效地處理這種情況。

為了有效地估計具有傳染性的疾病發(fā)病率，本文中將發(fā)病率考慮為變量，采用分層貝葉斯的方法來對發(fā)病率p進行統(tǒng)計推斷。為此，將能夠克服傳統(tǒng)beta分布缺點的廣義第二類Beta分布引入作為負二項分布發(fā)病率的先驗分布，構(gòu)建了分層貝葉斯模型，并提供了p的置信區(qū)間的構(gòu)造方法，通過計算機密集計算及MCMC方法來對其進行計算機實現(xiàn)。

[1]Bradlow,Eric T.,Bruce G.S.Hardie,Peter S.Fader.Bayesian Inference for the Negative Binomial Distribution Via Polynomial Expansions[J].Journal of Computational and Graphical Statistics,2002,11(1).

[2]Barlow,R.E,Proschan,F.Statistical Theory of Reliability and Life Testing:Probability Models[M].New York:Holt,Rinehart and Winston，1975.

[3]David L,Libby,Melvin R,Novick.Multivariate Generalized Beta Distributions with Applications to Utility Assessment[J].Journal of Educational Statistics,1982，7（4）.

[4]D.H.Young.Some Results for the Order Statistics of the Negative Binomial Distribution Under Slippage Configuration[J].Biometrika,1973，60(1).

[5]Gelman,A.,Carlin,J.,Stern,H.,Rubin,D.Bayesian Data Analysis(2ndEdition)[M].New York:Chapman and Hall/CRC Press，2003.

[6]Gradshteyn,I.S.,Ryzhik,I.M.Table of Integrals,Series,and Products(6thEdition)[M].San Diego:Academic Press，2000.

[7]Gupta,A.K,Nadarajah,S.Handbook of Beta Distribution and Its Applications[M].New York:Marcel Dekker，2004.

[8]Howard Raiffa,Robert Schlaifer.Applied Statistical Decision Theory.Division of Research,Graduate School of Business Administration[M].Boston:Harvard University，1961.

[9]James J,Chen,Melvin R.Novick.Bayesian Analysis for Binomial Models with Generalized Beta Prior Distributions[J].Journal of E-ducational Statistics，1984，9（2）.

[10]Jeff Miller，et al.Earliest Known Uses of Some of the Words of Mathematics,"Conjugate Prior Distributions"[J].Electronic Document,Revision of November,2005,（13）.

[11]Johnson,N.L.,Kotz,S.,Balakrishnan,N.Continuous Univariate Distributions(2ndEdition)[M].New York:John Wiley and Sons，1995.

[12]Newbold,E.Practical Applications of the Statistics of Repeated events,Particularly to IndustrialAccidents[J].J.Roy.Statist,1927，(90).

[13]Saralees Nadarajah,Samuel Kotz.A Generalized Beta Distribution II[J].Joural of Econometrice,2003，3.