邢伊琦,吳劉倉,聶興鋒

(昆明理工大學(xué)理學(xué)院,云南 昆明650093)

1.引言

經(jīng)典的線性模型可以描述許多研究領(lǐng)域的現(xiàn)象,但是,在生物、醫(yī)學(xué)、保險等領(lǐng)域的研究人員發(fā)現(xiàn),經(jīng)典的線性模型會遇到很多方法上的困難,在此基礎(chǔ)上,廣義線性模型也就應(yīng)運而生了.廣義線性模型是經(jīng)典線性模型的推廣.1972年Nelder和Wedderburn首先提出廣義線性模型[1],McCullagh和Nelder在1983年出版了系統(tǒng)論述此專題的專著《Generalized Linear Models》.在一些經(jīng)濟領(lǐng)域和工業(yè)產(chǎn)品質(zhì)量改進中,對均值和方差同時建模是非常有必要的,所以雙重廣義線性模型近年來引起了許多學(xué)者的關(guān)注.1984年,Pregibon在文章中首先提出了對散度參數(shù)建模的廣義線性模型[2];WANG和ZHANG[3]研究了雙重廣義線性模型中僅均值模型的變量選擇;WU和LI[4]研究了逆高斯分布下聯(lián)合均值和散度模型的變量選擇;陳海露[5]研究了雙重廣義線性模型的參數(shù)估計與變量選擇;吳劉倉等[6]研究了基于Box-Cox變換下聯(lián)合均值與散度廣義線性模型的極大似然估計;胡江等[7]研究了基于Pena距離的廣義線性回歸模型的影響分析等.

在響應(yīng)變量分布未知但已知其前兩階矩存在的情況下,Wedderburn提出了擬似然方法來進行參數(shù)估計.擬似然的方法是假定總體前兩階矩陣存在,然后通過對其對數(shù)似然方程求極值得到參數(shù)的估計值.陳希孺[8]在廣義線性模型中對擬似然方法有詳細(xì)的闡述;吳劉倉等[9]研究了缺失數(shù)據(jù)下雙重廣義線性模型的參數(shù)估計;袁巧莉等[10]研究了混合雙重廣義線性模型的參數(shù)估計等.

我們知道,統(tǒng)計診斷在數(shù)據(jù)分析中占有舉足輕重的地位,主要目的就是找出數(shù)據(jù)中的異常點或強影響點,常用的統(tǒng)計量有似然距離、Cook距離等,Pena距離[11]這一診斷統(tǒng)計量是美國統(tǒng)計學(xué)教授Daniel Pena在2005年首次提出的,并對其在線性模型上的影響分析做了詳細(xì)的研究,這種方法與之前的方法有較大差別,之前的方法是研究刪除一(組)點,對回歸分析的影響以及對預(yù)測值的影響,或者是某個(組)樣本點的微小擾動對參數(shù)估計的影響或是對模型預(yù)測的影響.而Pena距離這一統(tǒng)計量則是研究的是樣本中的某一個(組)點受其余各個(組)點的影響,簡單來說,就是樣本中各點刪除后,對某一特定的點的回歸值或預(yù)測值的影響.本文基于Pena距離,采用偽似然和擴展擬似然的方法估計參數(shù),并通過數(shù)據(jù)刪除模型的參數(shù)估計和統(tǒng)計診斷,比較了刪除模型和未刪除模型對應(yīng)的統(tǒng)計量之間的差異.通過Monte Carlo模擬驗證,本文提出方法的有效性.最后,通過實例研究,表明本文所提出的模型和方法是實用可行的.

2.雙重廣義線性模型下的極大似然估計

我們先給出雙重廣義線性模型(Double Generalized Linear Model,DGLM)為:

其中yi為被解釋變量,xi=(xi1,xi2,...,xip)T,zi=(zi1,zi2,...,zip)T為解釋變量,β=(β1,β2,...,βp)T為均值模型中的未知參數(shù),γ=(γ1,γ2,...,γp)T為散度模型中的未知參數(shù).xi,zi兩個解釋變量可能完全相同,部分相同或者完全不同,但相同的解釋變量在均值模型和方差模型中有不同的影響方式.g(μi)=xTi β是均值模型,h(?i)=zTi γ是散度模型,g(·)和h(·)是聯(lián)系函數(shù),V(·)是關(guān)于均值的方差函數(shù).

定理2.1對于模型(2.1),其Pena距離為:

證根據(jù)文[11],我們定義Pena距離如下:

對上式兩邊求方差，有：

上述兩式相減,得：

故

故Pena距離為：

3.雙重廣義線性模型的統(tǒng)計診斷

數(shù)據(jù)刪除是統(tǒng)計診斷中最常用的也是最基本的方法之一,比較刪除第i個點前后模型參數(shù)估計量之間的差異,能得出一些的結(jié)論.用這些結(jié)論,我們能評價我們估計方法的好壞,詳細(xì)內(nèi)容參考韋博成[12]等的文獻或書刊.模型(2.1)的數(shù)據(jù)刪除模型可表示為：

對于未刪除模型(2.1)和刪除模型(3.1),為檢驗第i個數(shù)據(jù)點在整個數(shù)據(jù)集中是否為異常點或強影響點,可通過比較刪除第i個點前后統(tǒng)計推斷結(jié)果的變化,看出這個點是否為異常點或者強影響點,而統(tǒng)計推斷結(jié)果的變化可由統(tǒng)計推斷量來得到.

對于一組隨機變量y1,y2,···,yn的分布是未知的,但知道其期望和方差存在,期望我們用E(y)表示,那么方差為:

從式中,我們知道E(y)是一階原點矩,Var(y)是二階中心矩,我們也可以認(rèn)為方差是二階原點矩減去期望的平方.對于本文選用的模型,我們選用擴展擬似然算法和偽似然算法來進行我們的參數(shù)估計.本文采用的擴展擬似然函數(shù)(EQL)Q+為：

本文中采用的偽似然函數(shù)(PL)Qp為:

用Gauss-Newton迭代法可得到參數(shù)極大似然估計的估計值.設(shè)未刪除模型的參數(shù)估計值用,表示,則刪除模型的參數(shù)估計值用,表示,則有=(,)T.

由Gauss-Newton迭代法可得到參數(shù)極大似然估計的估計值和,但如果解釋變量的維數(shù)為二維或者高于二維，這時參數(shù)和均為向量,難以比較大小.這時我們就可以用Cook距離來刻畫參數(shù)的變化,Cook距離定義如下:

其中H=(xT,zT)T為解釋變量,p為對應(yīng)解釋變量的維數(shù),為未刪除模型方差的估計值.

在分析具體數(shù)據(jù)時,先計算出各點的Cook距離,通過畫散點圖,找出其中特別大的Di,對應(yīng)數(shù)據(jù)點可能就是異常點或強影響點.

Pena距離與Cook距離相比較，前者研究的是刪除一個(組)點后對估計值或預(yù)測值的影響,而后者則研究的是樣本中的某一點受其余各點的影響,簡單的說,就是研究樣本中各點刪除后,對某一特定的點的估計值或預(yù)測值的影響,Pena距離定義如下：

其中H=X(XTX)?1XT稱為帽子矩陣,p為相應(yīng)解釋變量的維數(shù),為刪除第i個點后模型方差的估計值.是刪除第j個點后第i個點的參數(shù)估計值.具體分析時,同樣是先算出刪除各點后某一點的Si,畫出散點圖,其中Si較大的就可能是異常點或強影響點.

4.Monte Carlo模擬

為了比較Pena距離和Cook距離的診斷效果,本文我們采用Extra-Poisson模型進行模擬.Extra-Poisson 模型如下

根據(jù)模型(4.1)產(chǎn)生模擬數(shù)據(jù),其中yi是根據(jù)雙重廣義回歸模型產(chǎn)生的相互獨立的響應(yīng)變量,解釋變量xi和zi相互獨立產(chǎn)生于U(0,1).給定β0和γ0的真值分別為β0=(0,1,1)T,γ0=(0,1,1)T.將第170號,190 號樣本點的被解釋變量的值做改變,即從樣本點中人為的制造兩個異常點,然后應(yīng)用本文研究的Pena距離以及Cook 距離進行診斷,根據(jù)異常點的診斷情況來判斷本文提出的方法是否行之有效.模擬結(jié)果如圖1-4所示：

從圖中我們可以看出,無論是用PL或者EQL方法,第170號點以及190號點均被診斷出來了,這說明本文提出的方法是可行并且有效的,下面用具體的實例進一步說明.

圖1 PL的Cook距離CD散點圖

圖2 PL的Pena距離PD散點圖

圖3 EQL的Cook距離CD散點圖

圖4 EQL的Pena距離PD散點圖

5.實例分析

這里我們用一組檢驗?zāi)撤N工業(yè)用發(fā)動機性能試驗的數(shù)據(jù),該試驗使用的原料是柴油和從有機原料中通過蒸餾產(chǎn)生的氣體的混合物,在各種不同的速度x(計量單位：百轉(zhuǎn)/分鐘)下，測量發(fā)動機的馬力y.

表1 發(fā)動機性能數(shù)據(jù)

我們建立模型

利用偽似然和擴展擬似然的估計方法得到的參數(shù)做統(tǒng)計診斷，得到圖5-8的結(jié)果.

圖5 發(fā)動機性能數(shù)據(jù)PL的Cook距離CD散點圖

圖6 發(fā)動機性能數(shù)據(jù)PL的Pena距離PD散點圖

圖7 發(fā)動機性能數(shù)據(jù)EQL的Cook距離CD散點圖

圖8 發(fā)動機性能數(shù)據(jù)EQL的Pena距離PD散點圖

從圖5-8的結(jié)果可以看出,用偽似然和擴展擬似然的方法估計參數(shù)做統(tǒng)計診斷的效果大致相同.我們以圖5和圖6用偽似然的方法估計的參數(shù)做的統(tǒng)計診斷為例,由圖5可知第2號點、16號點和24號點可能為異常點或強影響點,由圖6可知第2號點和24號點可能為異常點或強影響點.由文[12]中的例5.4 可知,第2、24號點為異常點或強影響點,比起Cook距離,Pena距離很好地診斷出了這個點.

6.結(jié)論

本文針對分布未知但其一階矩和二階矩存在的隨機變量,建立了雙重廣義線性模型,運用擴展擬似然和偽似然方法進行參數(shù)估計再用Pena距離和Cook距離進行統(tǒng)計診斷,得到在一定的條件下Pena距離優(yōu)于Cook距離的結(jié)論.最后,通過Monte Carlo模擬和實例研究的結(jié)果驗證,說明了本文所提出的模型與方法的有效性和實用性.