趙為華，張日權(quán)

(1. 南通大學(xué) 理學(xué)院，江蘇 南通226019； 2. 華東師范大學(xué) 金融與統(tǒng)計學(xué)院，上海200241)

?

Beta-Binomial回歸模型及其應(yīng)用

趙為華1，張日權(quán)2

(1. 南通大學(xué) 理學(xué)院，江蘇 南通226019； 2. 華東師范大學(xué) 金融與統(tǒng)計學(xué)院，上海200241)

摘要：在成敗型試驗中或滿意度支持率調(diào)查中，Beta-Binomial分布常被用來刻畫具有偏大離差的計數(shù)型比例數(shù)據(jù)，由此提出Beta-Binomial 回歸模型，研究參數(shù)的最大似然估計方法并基于Newton-Raphson算法給出參數(shù)估計的迭代方法；重點討論模型中回歸參數(shù)和相關(guān)性參數(shù)存在的檢驗問題，提出Score檢驗方法并通過數(shù)值模擬研究Score檢驗統(tǒng)計量的檢驗功效問題；實例分析證明Beta-Binomial 回歸模型的有用性。

關(guān)鍵詞：比例數(shù)據(jù)；Beta-Binomial 分布；偏大離差；相關(guān)性

一、引 言

在許多生物學(xué)和醫(yī)學(xué)試驗或某項政策的社會滿意度調(diào)查中，常常得到的數(shù)據(jù)是計數(shù)型比例數(shù)據(jù)的形式，即Y/m，其中m是試驗次數(shù)(調(diào)查的人數(shù))，Y是試驗中成功的次數(shù)(支持政策的人數(shù))。當m固定時，隨機變量Y可以看成是m個相互獨立的Bernoulli隨機變量的和，即Y=W1+W2…+Wm服從二項分布(binomialdistribution)，記為 b(m,p)。然而，在一些實際問題中 W1,W2,…,Wm之間的獨立性并不能得到滿足，例如在抽樣調(diào)查中，當兩個樣本來自于同一個家庭或行業(yè)時，回答的結(jié)果就不再具有獨立性，而往往存在相關(guān)性， Y=W1+W2+…+Wm就不再服從二項分布。此時，如果仍然基于二項分布預(yù)測試驗成功率或政策支持率的話，會導(dǎo)致估計偏差較大、預(yù)測不準確等弱點。

E(Y)=mp

Var(Y)=mp(1-p)[1+(m-1)δ]

(1)

很明顯當δ=0，式(1)退化成二項分布的方差，而當δ>0時可以用來刻畫偏大離差計數(shù)數(shù)據(jù)，而當δ<0時可以用來刻畫偏小離差計數(shù)數(shù)據(jù)(實際情況中很少發(fā)生)。

為刻畫具有相關(guān)性Bernoulli隨機變量和的分布，許多研究者提出使用混合分布(mixture)或復(fù)合分布來描述響應(yīng)計數(shù)數(shù)據(jù)，其中最常用的方法是Beta-Binomial分布，即Y在給定參數(shù)p時服從二項分布 Y|p～b(m,p)，同時p服從beta分布，p～beta(a,b)，其中a>0和b>0為兩個未知參數(shù)[1]。眾所周知，beta分布是二項分布p的共軛先驗分布，且Y的分布列為：

(y=1,2，…,m)

(2)

以及Beta-Binomial分布的期望和方差為：

E(Y)=ma(a+b)-1

Var(Y)=mab(a+b)-2[1+(m-1)(a+b+1)-1]

(3)

令θ=a/(a+b)-1，δ=(a+b+1)-1，由式(1)可知，Beta-Binomial分布可以用來刻畫具有偏大離差的計數(shù)數(shù)據(jù)且Wj之間的相關(guān)性為δ=(a+b+1)-1。已有文獻已對Beta-Binomial模型展開過一些研究，分別將之應(yīng)用到生物學(xué)、生物化學(xué)和社會調(diào)查學(xué)等領(lǐng)域中[2-4]；申躍等基于經(jīng)驗Bayes估計方法將Beta-Binomial應(yīng)用到商品售后抱怨率問題中[5]；Paulo主要介紹了如何使用Stata軟件估計Beta-Binomial分布中的未知參數(shù)[6]；趙為華等則應(yīng)用EM算法研究Beta-Binomial模型中未知參數(shù)的估計方法，并將之應(yīng)用到學(xué)生考試不及格率和優(yōu)秀率的估計問題[7]。

然而，已有文獻主要研究Beta-Binomial模型及其應(yīng)用問題，很少研究基于Beta-Binomial分布的回歸建模及其統(tǒng)計推斷問題。本文將基于Newton-Raphson算法研究Beta-Binomial回歸模型的參數(shù)估計問題，并重點討論回歸參數(shù)和相關(guān)性參數(shù)的檢驗問題。

二、Beta-Binomial 回歸模型及其參數(shù)估計

(4)

其中β=(β1,β2,…,βp)T是p維回歸系數(shù)。 如果假定自變量中的第一個分量xi1≡1，則模型(4)中包含了常數(shù)項。

根據(jù)式(2)，在忽略與參數(shù)無關(guān)的常數(shù)后，Beta-Binomial回歸模型的對數(shù)似然函數(shù)可表示為：

l(β,φ)=

(5)

其中φ=δ(1-δ)-1，當δ=0亦即φ=0時，Beta-

Binomial的對數(shù)似然退化成Binomial回歸的對數(shù)似然。

由式(5)，得到參數(shù)β和φ的得分函數(shù)為：

(6)

進一步，可得觀測Fisher信息陣為：

(7)

(i=1,2,…,n)

記η=(βT,φ)T。由式(6)和(7)，根據(jù)Newton-Raphson算法得到迭代公式為：

η(k+1)=η(k)+I-1S|η(k)

三、參數(shù)的Score檢驗

在實際回歸建模時，需要評價模型的正確性和模型中自變量的重要性。為此下面應(yīng)用Score檢驗統(tǒng)計量對相關(guān)性參數(shù)φ和回歸系數(shù)β的重要性進行檢驗。Score檢驗統(tǒng)計量最大的優(yōu)點是只需要在原假設(shè)成立情況下計算參數(shù)估計。

對于相關(guān)性參數(shù)φ，本文提出檢驗問題：

H01:φ=0?H11:φ>0

(8)

(9)

其中Iφφ為觀測Fisher信息陣的逆矩陣對應(yīng)參數(shù)φ的分塊。通過一些簡單推導(dǎo)，檢驗統(tǒng)計量SC1漸近服從自由度為1的卡方分布χ2(1)。

H02:βb=0?H12:βb≠0

(10)

(11)

四、數(shù)值模擬

下面通過隨機模擬來說明估計方法的有效性以及檢驗統(tǒng)計量的檢驗功效。響應(yīng)變量數(shù)據(jù)由如下模型隨機產(chǎn)生：

yi|pi～b(mi,pi)，pi～beta(θiω,(1-θi)ω)

表1　不同樣本量下參數(shù)β和φ的偏差、絕對值偏差及其標準差估計表

表1的結(jié)果表明，估計量的偏差都較小，且隨著樣本量n的增加估計的絕對值偏差越來越接近于0。此外，不同樣本量下估計量的標準差估計和樣本標準差估計非常接近，且都隨著n的增大接近于0。圖 1 顯示了回歸系數(shù)β和相關(guān)性參數(shù)φ=1/ω估計的均方誤差(MSE)箱線圖，易見隨著樣本量的增大，本文的估計效果是令人滿意的。

下面研究Score檢驗統(tǒng)計量的檢驗功效問題。在數(shù)據(jù)生成時，其他參數(shù)保持不變的情況下，分別取β5=0，0.1，0.2，0.3，0.4，0.5和1時，考察檢驗統(tǒng)計量(11)的檢驗功效，即在顯著性水平0.05下，計算1 000次模擬中拒絕原假設(shè)H0:β5=0的比例，結(jié)果見表2和圖2。

圖1　參數(shù)估計MSE的箱線圖

圖2　參數(shù)β5的檢驗功效圖

從表2和圖2易見，在相同的樣本量下，隨著參數(shù)β5的取值遠離0，檢驗的功效顯著增加；另外一方面，隨著樣本量的增大，檢驗的功效迅速接近于1，且在原假設(shè)正確時(β5=0)，檢驗的功效非常接近于名義水平0.05。

相關(guān)性參數(shù)φ的檢驗功效見表3和圖3，其中模擬中ω分別取1 000，500，100，50，20，10，5，3，對應(yīng)的相關(guān)性φ=1/1 000，1/500，1/250,1/100，1/50，1/20，1/10，1/3由弱到強(φ=1/1 000可近似地認為不相關(guān))。模擬結(jié)果再次表明：參數(shù)接近于原假設(shè)時，檢驗的功效與名義水平非常接近，當遠離原假設(shè)時，檢驗的功效隨著樣本量的增大迅速提高。

表2　不同樣本量下β5的檢驗功效表

表3　不同樣本量下φ的檢驗功效表

圖3　參數(shù)φ的檢驗功效圖

五、實例分析

Wheatley等研究和報道了歐洲胡蘿卜蠅對植物胡蘿卜生長的破壞數(shù)據(jù)[8]。該數(shù)據(jù)共有36 組數(shù)據(jù)，主要是用來研究噴灑農(nóng)藥的劑量( 共5 種不同劑量)和噴灑殺蟲劑農(nóng)藥的類型(即二嗪農(nóng)殺蟲劑(diazinon)、乙拌磷殺蟲劑(disulfoton)以及不噴灑農(nóng)藥)對控制胡蘿卜蠅破壞胡蘿卜植物生長是否有顯著影響。以mi表示進入研究的植物數(shù)目，yi表示受到破壞的植物數(shù)目，xi1表示噴灑農(nóng)藥的劑量，xi2和xi3表示0-1型分類變量，即：

(i=1,2，…,36)

由于該數(shù)據(jù)存在明顯的偏大離差，Hardin等提出了具有偏大離差的計數(shù)回歸模型擬合該數(shù)據(jù)[9]160-161。這里，使本文提出的Beta-Binomial回歸模型(4)再次分析該數(shù)據(jù)，其結(jié)果見表4(為比較起見，表中亦包含Binomial回歸模型擬合結(jié)果)。

表4　歐洲胡蘿卜蠅數(shù)據(jù)的擬合結(jié)果表

從表4 不難看出，噴灑殺蟲劑的種類以及是否噴灑農(nóng)藥有顯著差異，且系數(shù)β1和β2的估計是負值，說明噴灑對有效防治胡蘿卜蠅的破壞性有著重要作用；另外一方面，系數(shù)β3的檢驗p值在兩個模型下均大于0.05，說明增加噴灑劑量對防治胡蘿卜蠅破壞作用不敏感。同時，模型的擬合優(yōu)度BIC值清晰地表明Beta-Binomial模型好于Binomial回歸模型。

為進一步檢驗Beta-Binomial模型的有用性，我們計算Score檢驗統(tǒng)計量(9)，其結(jié)果為：

SC1=224.794 0>>χ2(1)=3.841 5

再次清楚地說明數(shù)據(jù)中存在著偏大離差，此時直接使用Binomial模型擬合數(shù)據(jù)會導(dǎo)致擬合不充分、推斷效果較差等問題。

六、總結(jié)

Beta-Binomial分布是一種混合分布模型，能夠處理具有偏大離差的計數(shù)型比例數(shù)據(jù)。本文基于Beta-Binomial分布提出了Beta-Binomial回歸模型，給出了模型中參數(shù)的最大似然估計方法，并基于Score經(jīng)驗統(tǒng)計量給出了重要變量和相關(guān)性參數(shù)的檢驗問題。數(shù)值模擬和實例分析充分說明了所提方法的有效性。進一步，為刻畫數(shù)據(jù)的相關(guān)性和異質(zhì)性，提高估計和預(yù)測的準確性，在Beta-Binomial模型中可以對均值參數(shù)和相關(guān)性參數(shù)進行同時回歸建模，并研究回歸系數(shù)的重要性檢驗問題以及對模型進行擬合優(yōu)度檢驗。

參考文獻：

[1]Crowder M. Beta-Binomial ANOVA for Proportions[J]. Applied Statistics, 1978, 27.

[2]Otake M, Prenticce R. The Analysis of Chromosomally Aberrant Cells Based on a Beta-Binomial Distribution[J]. Radiation Research, 1984, 98(3).

[3]Prentice R. Binary Regression Using an Extended Beta-Binomial Distribution, with Discussion of Correlation Induced by Covariate Measurement Errors[J]. Journal of the American Statistical Association, 1986, 81(394).

[4]Silverman B, Wilson J. A Beta-Binomial Model for Library Survey Data[J]. Journal of Documentation, 1987, 43(2).

[5]申躍, 趙平. Beta-Binomial 模型和經(jīng)驗貝葉斯方法在抱怨行為傾向性研究中的應(yīng)用[J]. 統(tǒng)計研究, 2004, 21( 11).

[6]Paulo, Guimaraes. A Simple Approach to Fit the Beta-Binomial Model[J]. Stata Journal, 2005, 5(3).

[7]趙為華, 李澤安, 陸志峰. Beta-Binomial 模型基于EM 算法的應(yīng)用研究[J]. 大學(xué)數(shù)學(xué), 2010, 26(3).

[8]Wheatley G, Freeman, G. A Method of Using the Proportion of Undamaged Carrots or Parsins to Estimate the Relative Population Densities of Carrot Fly( Psila Rosae) Larvae, and Its Practical Application[J]. Annals of Applied Biology, 1982,100 (2) .

[9]Hardin J, Hilbe J. Generalized Linear Models and Extensions[M]. 2nd. College Station，USA: Stata Press, 2007.

(責(zé)任編輯：崔國平)

Beta-Binomial Regression and Its Application

ZHAO Wei-hua1, ZHANG Ri-quan2

(1. School of Science, Nantong University, Nantong 226019, China; 2. School of Finance and Statistics, East China Normal University, Shanghai 200241, China)

Abstract:Beta-Binomial distribution is frequently used to describe the over-dispersion count data in the success-failure type test or satisfaction ratings survey problem. In this paper, we address the Beta-Binomial regression model, and the maximum likelihood estimate method is investigated, which can be implemented by the Newton-Raphson algorithm. Furthermore, the Score test methods are proposed to test the important covariates and the existing of the correlation parameter, and their test powers are investigated by simulations. Finally, the real data analysis is used to illustrate the usefulness of Beta-Binomial regression model.

Key words:proportional data; Beta-Binomial distribution; over-dispersion; correlation

中圖分類號：C81∶O212

文獻標志碼：A

文章編號：1007-3116(2016)03-0009-05

作者簡介：趙為華，男，江蘇海門人，理學(xué)博士，副教授，碩士生導(dǎo)師，研究方向：分位數(shù)回歸建模及其應(yīng)用；

基金項目：教育部人文社科青年基金項目《比例數(shù)據(jù)的分位數(shù)回歸建?！?14YJC910007)

收稿日期：2015-10-23

張日權(quán)，男，山西大同人，理學(xué)博士，教授，博士生導(dǎo)師，研究方向：分位數(shù)回歸和半?yún)?shù)模型。

【統(tǒng)計理論與方法】