關于傳染性負二項分布發病率的貝葉斯推斷

戴 成，陳博鈺，田茂再

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

0 引言

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

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

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

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

1 建立模型

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

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

1.1 廣義第二類beta分布

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

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

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

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

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

1.2 負二項分布X的邊際密度函數

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

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

此處我們可以通過高斯超幾何函數的特殊性質來獲得關于 (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為整數，則有

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

1.3 結構性質

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

1.3.1 均值與方差

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

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

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

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

由(14)和(16)式可得

1.3.2 矩生成函數

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

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

2 分層貝葉斯模型

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

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

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

由此可得

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

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

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

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

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

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

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

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

3 結論

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

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

[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.