謝遠(yuǎn)濤,楊 娟,徐梅笛
(1.對外經(jīng)濟(jì)貿(mào)易大學(xué)a.保險(xiǎn)學(xué)院;b.國際經(jīng)濟(jì)貿(mào)易學(xué)院,北京 100029;2.中國人民大學(xué)統(tǒng)計(jì)學(xué)院,北京 100872)
廣義Gamma分布簇廣義線性混合模型與廣義線性混合模型乃至廣義線性模型的關(guān)鍵區(qū)別是,變量不再滿足指數(shù)分布簇,因此,統(tǒng)計(jì)推斷問題變得異常復(fù)雜,本文主要討論廣義Gamma分布簇廣義線性混合模型的參數(shù)估計(jì)問題。
自1972年Nelder和Wedderburn引進(jìn)廣義線性模型以來,學(xué)者一直嘗試在更廣義的分布中建立模型,如指數(shù)分布簇,單參數(shù)Gamma分布(Clayton and Cuzick,1985),正則穩(wěn)態(tài)分布(Hougaard,1986),逆高斯分布(Hougaard,1986),對數(shù)正態(tài)分布(McGilchrist and Aisbett,1991)。
同時(shí),不斷引入異質(zhì)性來增強(qiáng)模型的適應(yīng)性。Zeger,Liang and Albert(1988)集中在縱向數(shù)據(jù)系統(tǒng)討論了隨機(jī)效應(yīng)模型的建模問題。Gurka et al.(2006)把Box-Cox變換推廣到線性混合模型,He and Lawless(2005)建立了一類位置-尺度模型。隨機(jī)效應(yīng)模型與重復(fù)觀測效應(yīng)結(jié)合起來分析,最終產(chǎn)生了廣義線性混合模型,模型參見Jiang(2007)。
國內(nèi)對廣義線性混合模型的研究文獻(xiàn)極少,謝遠(yuǎn)濤和楊娟(2010)構(gòu)造了廣義Gamma分布簇廣義線性混合模型,充分利用響應(yīng)變量之間的相關(guān)性結(jié)構(gòu),在變量滿足廣義Gamma分布簇時(shí)能有效降低模型誤設(shè)的風(fēng)險(xiǎn),但是作者沒有展開分析參數(shù)的估計(jì)問題。
參數(shù)估計(jì)等推斷的前提是構(gòu)造似然函數(shù),那么對隨機(jī)效應(yīng)的處理,無論是條件似然法還是極大似然法都會(huì)遇到很多麻煩。
Zeger,Liang and Albert(1988)集中在縱向數(shù)據(jù)系統(tǒng)討論了隨機(jī)效應(yīng)模型的建模問題。關(guān)于隨機(jī)效應(yīng)部分建模,這里有兩種不同的思路。第一種思路是把ui當(dāng)作擾動(dòng)(nuisance)變量,僅僅使用不包含ui信息的那部分?jǐn)?shù)據(jù)來對特定系數(shù)進(jìn)行推斷。這被稱為總體平均法(PA,population-averaged)。第二種思路被稱為個(gè)體設(shè)定方法(SS,subject-specific),該法重點(diǎn)關(guān)注單個(gè)個(gè)體隨機(jī)效應(yīng)u的估計(jì),以及它與總體參數(shù)也即固定效應(yīng)β之間的關(guān)系。這里把ui看作來自某一分布的若干獨(dú)立樣本,同時(shí)估計(jì)固定效應(yīng)βi和隨機(jī)效應(yīng)ui。
在隨機(jī)常數(shù)項(xiàng)模型中,在這兩種推斷方法之間的選擇導(dǎo)致使用橫截信息與使用縱向信息的不同。第一種思路僅僅使用了縱向信息,也即,在個(gè)體內(nèi)進(jìn)行比較來估計(jì)ui。而第二種思路中假定ui滿足一個(gè)分布,我們可以同時(shí)考慮縱向數(shù)據(jù)信息和橫截信息,相應(yīng)的每個(gè)數(shù)據(jù)源的權(quán)重由ui的變動(dòng)性來確定。
這兩種思路的選取將直接影響似然函數(shù)的構(gòu)造。第一種思路產(chǎn)生了偏似然估計(jì),包括條件似然估計(jì)和邊際似然估計(jì);第二種思路產(chǎn)生了完全似然估計(jì),后者要復(fù)雜很多。
一些學(xué)者選擇偏似然方法(Cox,1975)進(jìn)行邊際模型估計(jì),偏似然是邊際似然和條件似然的推廣。而Fraser(1968)、Kalbfleisch and Sprott(1974)早已對邊際似然和條件似然進(jìn)行了研究。
完全似然方法也有不少學(xué)者研究,其中最關(guān)鍵的一步是對似然函數(shù)的推廣。Wedderburn(1974)提出擬似然(quasi-likelihood)函 數(shù) ,在 McCullagh(1981)、McCullagh and Nelder(1989)的文獻(xiàn)中有系統(tǒng)論述。
擬似然的提出為一大類模型提供了推斷的基礎(chǔ)。對于大多數(shù)非正態(tài)模型,需要使用數(shù)值積分(如Crouch and Spiegelman,1990)。在階數(shù)取更高的時(shí)侯,可以使用Mon-te Carlo積分方法。Li and Raghunathan(1991)在先驗(yàn)分布中使用重要重復(fù)抽樣(或者Gibbs抽樣來回避數(shù)值積分(Zeger and Karim,1991)。對于求解極大似然估計(jì),最常用的是EM算法(Dempster等,1977)。
到了后來,一些基于近似的方法逐漸提出并占主流地位。這類近似方法包括Laplace近似方法(Tierney等1989)和Taylor序列展開式方法。后文會(huì)有詳細(xì)的論述。
這些參數(shù)估計(jì)理論非常豐富,但前提是,這些估計(jì)思想只適用于正態(tài)分布,或者指數(shù)分布簇,對于含隨機(jī)效應(yīng)項(xiàng)的廣義Gamma分布簇,尚缺乏研究結(jié)果。
本文將重點(diǎn)討論廣義Gamma分布簇廣義線性混合模型的參數(shù)估計(jì)問題。這類模型的參數(shù)估計(jì)理論沒有公開發(fā)表,這是本文的創(chuàng)新點(diǎn)。另外,本文在推斷上廣泛應(yīng)用廣義線性混合模型的研究成果,這是本文的特點(diǎn)。
廣義Gamma分布簇廣義線性混合模型可以用矩陣形式來表述如下:
其中,η為線性預(yù)測,等于連接函數(shù)g(μ),μ為條件均值。X是固定效應(yīng)解釋變量矩陣,β是固定效應(yīng)系數(shù)向量,Z是隨機(jī)效應(yīng)設(shè)計(jì)矩陣,u是隨機(jī)效應(yīng)系數(shù)向量,均值為0,協(xié)方差矩陣為G,后文中常常用到的假定是隨機(jī)效應(yīng)項(xiàng)滿足MVN(0,G)分布,e是擾動(dòng)項(xiàng),協(xié)方差矩陣為R。記為:
設(shè)該協(xié)方差結(jié)構(gòu)中包含的未知參數(shù)為φ。此框架下,可以同時(shí)設(shè)定隨機(jī)效應(yīng)的方差結(jié)構(gòu)G和重復(fù)觀測效應(yīng)方差結(jié)構(gòu)R。
如果yi是來自廣義Gamma分布的隨機(jī)樣本,構(gòu)造對數(shù)擬似然函數(shù)寫成標(biāo)準(zhǔn)參數(shù)形式為:
本文討論的參數(shù)估計(jì)方法是基于近似的一類估計(jì)方法。主要根據(jù)擬似然、偽似然、邊際似然構(gòu)建極大似然估計(jì)和受限極大似然估計(jì)。6種參數(shù)估計(jì)的方法,也即罰擬似然、邊際擬似然、偽似然的ML估計(jì)和REML估計(jì),盡管理論假設(shè)不同,分析思路不同,但是最終的估計(jì)形式驚人的一致。
罰擬似然(PQL)方法基本思想是把完全似然函數(shù)用Laplace方法來近似,通過對對數(shù)似然函數(shù)關(guān)于參數(shù)δ偏導(dǎo)并令其等于零來構(gòu)造得分方程,然后我們可以通過求解一些得分方程來找到極大似然估計(jì)。
邊際擬似然(MQL)方法的最大優(yōu)點(diǎn)是穩(wěn)健性,即使方差成分誤設(shè),廣義估計(jì)方程仍然能提供固定效應(yīng)的一致性估計(jì),當(dāng)然前提條件是一階矩設(shè)定正確它的一個(gè)特點(diǎn)是可重復(fù)性和容易實(shí)施,通過重復(fù)調(diào)用常規(guī)正態(tài)分布觀測值方差分析軟件就可以用來估計(jì)非正態(tài)分布模型,如重復(fù)調(diào)用廣義最小二乘來估計(jì)廣義線性模型。因此這種方法很流行。
Breslow and Clayton(1993)同時(shí)對兩種方法進(jìn)行系統(tǒng)介紹和比較,為把得分方程組的求解轉(zhuǎn)化為迭代方式求解提供了理論支持,具體求解使用的方法是迭代(再)加權(quán)最小二乘算法(Nelder and Wedderburn,1972)。Vonesh等(2002)同時(shí)擴(kuò)展了罰擬似然方法和廣義估計(jì)方程,在估計(jì)方程中引入了一二階條件矩。Lin等(1997)進(jìn)一步把協(xié)變量效應(yīng)引入到方差成分中。
Wolfinger and O’Connell(1993)使用的偽似然方法在實(shí)際編程應(yīng)用中有重要地位。偽似然(Carroll and Ruppert,1988,P71)是指“偽稱其他回歸參數(shù)β已知并且等于當(dāng)前的估計(jì)值進(jìn)一步基于正態(tài)性假定得到參數(shù)θ的極大似然估計(jì)。”偽似然實(shí)施過程中,是輪流估計(jì)參數(shù)β和θ,迭代直到收斂。如果在估計(jì)隨機(jī)效應(yīng)協(xié)方差矩陣G和殘差協(xié)方差矩陣R的參數(shù)時(shí)使用了ML方法,那就得到PL估計(jì)量;如果在估計(jì)隨機(jī)效應(yīng)協(xié)方差矩陣G和殘差協(xié)方差矩陣R的參數(shù)時(shí)使用了REML方法,那就得到REPL估計(jì)量。偽似然方法提供了SS推斷和PA推斷的統(tǒng)一的框架,并且把PQL和MQL作為特例包含進(jìn)來。而且,該方法可以允許自定義G和R的協(xié)方差結(jié)構(gòu)以適應(yīng)更一般的模型。更重要的一點(diǎn)是,PL或REPL方法可以很方便地利用SAS的固有線性混合模型模塊來實(shí)現(xiàn)。
Davidian and Giltinan(1993)討論了兩類非線性混合模型的估計(jì)方法,一種他們稱為pooled two stage(PTS)法,另外一種他們稱為線性化混合效應(yīng)(LME)方法。PTS法對從每個(gè)個(gè)體得到的估計(jì)進(jìn)行合成;LME方法是Vonesh and Carter(1992)線性化方法的推廣。Wolfinger and O’Connell(1993)偽似然方法與LME方法非常相似,不同之處在于,LME對于隨機(jī)效應(yīng)的線性化是0,但Wolfinger and O’Connell(1993)的方法以及Lindstrom and Bates(1990)的方法都對參數(shù)的當(dāng)前估計(jì)進(jìn)行線性化。
總結(jié)一下,Breslow and Clayton(1993)利用擬似然總結(jié)了PQL和MQL兩種方法,但最終以特例形式在Wolfinger and O’Connell(1993)的偽似然方法中體現(xiàn)出來。同時(shí),偽似然方法與LME方法思路非常相似,只是在線性近似處理上有所不同。
綜合這些模型,得出本模型參數(shù)估計(jì)的基本思路如下:
首先固定τ,偽稱其估計(jì)已知τ=τ?,然后討論其余參數(shù)的估計(jì),最后基于搜索算法求的使所有參數(shù)估計(jì)最優(yōu)的τ?。
具體說來,通過使用Laplace近似的方法(Tierney等1989)或者Taylor序列近似,可以得到目標(biāo)似然函數(shù)的線性近似或者分析近似,構(gòu)造出得分方程組,求解都可以運(yùn)用Fisher得分法得到迭代(再)加權(quán)最小二乘估計(jì)。也即可以通過下面的方程組以迭代方式求解β和u:
然后可以得到隨機(jī)效應(yīng)的收縮估計(jì):
迭代(再)加權(quán)最小二乘法迭代求解的時(shí)候要求G?非奇異。但是在迭代運(yùn)算的過程中,常常會(huì)出現(xiàn)迭代估計(jì)失敗造成方差成分取邊界值0。如果奇異,需要使用Cholesky分解進(jìn)行調(diào)整(Henderson 1984)。
設(shè)是?的Cholesky下三角根那么迭代(再)加權(quán)最小二乘法方程組變?yōu)椋?/p>
而且,這個(gè)結(jié)論還可以推廣,如果把V用一致性估計(jì)替換,以上漸進(jìn)性依然成立。
把這些參數(shù)估計(jì)帶入到目標(biāo)似然函數(shù),得到調(diào)整后參數(shù)?的近似概括擬似然函數(shù),利用這個(gè)公式關(guān)于?偏導(dǎo)并令其為0,得到得分方程,基于此可以求解
在建模過程中假定尺度參數(shù)φ=1可以得到Breslow and Clayton的罰擬似然方法(penalized quasi-likelihood,PQL);而 Wolfinger and O’Connell方法為偽似然方法(pseudo-likelihood,PL)或者受限偽似然方法(restricted pseudo-likelihood,REPL),在建模過程中假定尺度參數(shù)φ是未知的,因此利用PL得到尺度參數(shù)的極大似然估計(jì)或者利用REPL方法獲得尺度參數(shù)的受限極大似然估計(jì)。從這種意義上說,罰擬似然方法只是偽似然方法中尺度參數(shù)φ=1的特例。
利用前面的方法可以得出固定τ下的參數(shù)估計(jì),假定為,把它們帶回聯(lián)合似然函數(shù),得到聯(lián)合似然函數(shù)的估計(jì)式
因?yàn)閥*=yτ,在τ>0時(shí)為單調(diào)函數(shù),那么τ的估計(jì)可以直接對聯(lián)合似然函數(shù)極大化來求解:
本文設(shè)計(jì)非對稱有記憶二分搜索算法:
圖1 直接搜索算法流程圖
縱觀算法,運(yùn)算不外乎兩類:一是τ的變化,二是步長d的變化。我們把τ的變化(加減d個(gè)步長)定義為延展,把τ不變化而步長改變定義為收縮。一旦運(yùn)算進(jìn)入收縮階段,就會(huì)收斂到最小值,當(dāng)然,這個(gè)最小值可能是局部最小值。
初始化時(shí)可以選擇τ=1,初始化步長d,d的取值不易大于1,防止τ取負(fù),同時(shí),d的取值不易太小,以防解收斂于局部最小值;該算法對往正方向移動(dòng)還是往負(fù)方向移動(dòng)是非對稱的。k為步率,反映了調(diào)整步長的程度。當(dāng)τ的解往正方向或負(fù)方向移動(dòng)時(shí),步長不改變,恒定為;但是,一旦發(fā)現(xiàn)無論是正方向移動(dòng)還是負(fù)方向移動(dòng)都失敗時(shí)(極值點(diǎn)出現(xiàn)在當(dāng)前估計(jì)值的(τ-d,τ+d)區(qū)間上),程序會(huì)自動(dòng)調(diào)整步長為原來步長的k倍,以防止出現(xiàn)來回震蕩的情況。
本文提出了廣義Gamma分布簇廣義線性混合模型的參數(shù)估計(jì)方法,該方法充分利用了指數(shù)分布簇廣義線性混合模型參數(shù)估計(jì)的結(jié)果,并且參數(shù)估計(jì)表達(dá)式上與6種廣義線性混合模型參數(shù)估計(jì)結(jié)果形式相同,但是模型假設(shè)和表示的含義有異。而且,該方法可以通過反復(fù)調(diào)用廣義線性混合模型中的參數(shù)估計(jì)模塊來實(shí)現(xiàn),具有漸進(jìn)正態(tài)性,因此有較好的易用性。
[1]謝遠(yuǎn)濤,楊娟.廣義Gamma分布簇廣義線性混合模型的構(gòu)建[J].統(tǒng)計(jì)研究,2010,(10).
[2]Nelder J A,Wedderburn R W M.Generalized Linear Models[J].Jour?nal of the Royal Statistical Society,Ser.A,1972,(135).
[3]Clayton D G,Cuzick J.Multivariate Generalizations of the Proportion?al Hazards Model(with Discussion)[J].Journal of the Royal Statistical Society,Ser A,1985,(148).
[4]Hougaard P.A Class of Multivariate Failure Time Distributions[J].Biometrika,1986,(73).
[5]Hougaard P.Survival Models for Heterogeneous Populations Derived from Stable Distributions[J].Biometrika,1986,(73).
[6]McCullagh P,Nelder J A.Generalized Linear Models(2ndEdtion)[M].London:Chapman&Hall/CRC,1989.
[7]Zeger S L,Liang K Y,Albert P S.Models for Longitudinal Data:a Generalized Estimating Equation Approach[J].Biometrics,1988,(44).
[8]Geert Verbeke,Geert Molenberghs.Linear Mixed Models for Longitu?dinal[M].New York:Springer,2000.
[9]Gurka M J,Edwards L J,Muller K E,Kupper L L.Extending the Box-Cox Transformation to the Linear Mixed Model[J].Journal of the Royal Statistical Society,Ser,A,2006,(169).
[10]He W Q ,Lawless J F.Bivariate Location-Scale Models for Regres?sion Analysis,with Applications to Lifetime Data[J].Journal of the Royal Statistical Society,Ser,B,2005,(67).
[11]Jiang J.Linear and Generalized Linear Mixed Models and their Ap?plications[M].New York:Springer,2007.
[12]McCullagh P,Nelder J A.Generalized Linear Models(2ndEditon)[M].London:Chapman and Hall,1989.
[13]Breslow N E,Clayton D G.Approximate Inference in Generalized Linear Mixed Models[J].JASA,1993,(88).
[14]Wolfinger R ,O’Connell M.Generalized Linear Mixed Models:a Pseudo-likelihood Approach[J].J.Statistics,Computation simula?tion,1993,(48).
[15]Carroll R J ,Ruppert D.Transformation and Weighting in Regres?sion[M].London:Chapman and Hall,1988.
[16]Harville D A.Maximum Likelihood Approaches to Variance Compo?nent Estimation and to Related Problems[J].Journal of the American Statistical Association,1977,(72).