廣義Gamma分布簇廣義線性混合模型的參數估計

謝遠濤，楊 娟,徐梅笛

（1.對外經濟貿易大學a.保險學院;b.國際經濟貿易學院，北京 100029；2.中國人民大學統計學院，北京 100872）

0 引言

廣義Gamma分布簇廣義線性混合模型與廣義線性混合模型乃至廣義線性模型的關鍵區別是，變量不再滿足指數分布簇，因此，統計推斷問題變得異常復雜，本文主要討論廣義Gamma分布簇廣義線性混合模型的參數估計問題。

自1972年Nelder和Wedderburn引進廣義線性模型以來，學者一直嘗試在更廣義的分布中建立模型，如指數分布簇，單參數Gamma分布（Clayton and Cuzick,1985），正則穩態分布（Hougaard,1986），逆高斯分布(Hougaard,1986)，對數正態分布（McGilchrist and Aisbett,1991）。

同時，不斷引入異質性來增強模型的適應性。Zeger,Liang and Albert(1988)集中在縱向數據系統討論了隨機效應模型的建模問題。Gurka et al.(2006)把Box-Cox變換推廣到線性混合模型，He and Lawless(2005)建立了一類位置-尺度模型。隨機效應模型與重復觀測效應結合起來分析，最終產生了廣義線性混合模型，模型參見Jiang（2007）。

國內對廣義線性混合模型的研究文獻極少，謝遠濤和楊娟（2010）構造了廣義Gamma分布簇廣義線性混合模型，充分利用響應變量之間的相關性結構，在變量滿足廣義Gamma分布簇時能有效降低模型誤設的風險，但是作者沒有展開分析參數的估計問題。

參數估計等推斷的前提是構造似然函數，那么對隨機效應的處理，無論是條件似然法還是極大似然法都會遇到很多麻煩。

Zeger,Liang and Albert(1988)集中在縱向數據系統討論了隨機效應模型的建模問題。關于隨機效應部分建模，這里有兩種不同的思路。第一種思路是把ui當作擾動（nuisance）變量，僅僅使用不包含ui信息的那部分數據來對特定系數進行推斷。這被稱為總體平均法（PA，population-averaged）。第二種思路被稱為個體設定方法（SS，subject-specific），該法重點關注單個個體隨機效應u的估計，以及它與總體參數也即固定效應β之間的關系。這里把ui看作來自某一分布的若干獨立樣本，同時估計固定效應βi和隨機效應ui。

在隨機常數項模型中，在這兩種推斷方法之間的選擇導致使用橫截信息與使用縱向信息的不同。第一種思路僅僅使用了縱向信息，也即，在個體內進行比較來估計ui。而第二種思路中假定ui滿足一個分布，我們可以同時考慮縱向數據信息和橫截信息，相應的每個數據源的權重由ui的變動性來確定。

這兩種思路的選取將直接影響似然函數的構造。第一種思路產生了偏似然估計，包括條件似然估計和邊際似然估計；第二種思路產生了完全似然估計，后者要復雜很多。

一些學者選擇偏似然方法（Cox,1975）進行邊際模型估計，偏似然是邊際似然和條件似然的推廣。而Fraser（1968）、Kalbfleisch and Sprott(1974)早已對邊際似然和條件似然進行了研究。

完全似然方法也有不少學者研究，其中最關鍵的一步是對似然函數的推廣。Wedderburn(1974)提出擬似然（quasi-likelihood）函 數 ，在 McCullagh(1981)、McCullagh and Nelder(1989)的文獻中有系統論述。

擬似然的提出為一大類模型提供了推斷的基礎。對于大多數非正態模型，需要使用數值積分（如Crouch and Spiegelman,1990）。在階數取更高的時侯，可以使用Mon-te Carlo積分方法。Li and Raghunathan（1991）在先驗分布中使用重要重復抽樣（或者Gibbs抽樣來回避數值積分（Zeger and Karim，1991）。對于求解極大似然估計，最常用的是EM算法（Dempster等，1977）。

到了后來，一些基于近似的方法逐漸提出并占主流地位。這類近似方法包括Laplace近似方法（Tierney等1989）和Taylor序列展開式方法。后文會有詳細的論述。

這些參數估計理論非常豐富，但前提是，這些估計思想只適用于正態分布，或者指數分布簇，對于含隨機效應項的廣義Gamma分布簇，尚缺乏研究結果。

本文將重點討論廣義Gamma分布簇廣義線性混合模型的參數估計問題。這類模型的參數估計理論沒有公開發表，這是本文的創新點。另外，本文在推斷上廣泛應用廣義線性混合模型的研究成果，這是本文的特點。

1 模型與對數擬似然函數的構建

廣義Gamma分布簇廣義線性混合模型可以用矩陣形式來表述如下：

其中，η為線性預測，等于連接函數g(μ)，μ為條件均值。X是固定效應解釋變量矩陣，β是固定效應系數向量，Z是隨機效應設計矩陣，u是隨機效應系數向量，均值為0，協方差矩陣為G，后文中常常用到的假定是隨機效應項滿足MVN（0，G）分布，e是擾動項，協方差矩陣為R。記為：

設該協方差結構中包含的未知參數為φ。此框架下，可以同時設定隨機效應的方差結構G和重復觀測效應方差結構R。

如果yi是來自廣義Gamma分布的隨機樣本，構造對數擬似然函數寫成標準參數形式為：

2 參數估計

本文討論的參數估計方法是基于近似的一類估計方法。主要根據擬似然、偽似然、邊際似然構建極大似然估計和受限極大似然估計。6種參數估計的方法，也即罰擬似然、邊際擬似然、偽似然的ML估計和REML估計，盡管理論假設不同，分析思路不同，但是最終的估計形式驚人的一致。

罰擬似然（PQL）方法基本思想是把完全似然函數用Laplace方法來近似，通過對對數似然函數關于參數δ偏導并令其等于零來構造得分方程，然后我們可以通過求解一些得分方程來找到極大似然估計。

邊際擬似然（MQL）方法的最大優點是穩健性，即使方差成分誤設，廣義估計方程仍然能提供固定效應的一致性估計，當然前提條件是一階矩設定正確它的一個特點是可重復性和容易實施，通過重復調用常規正態分布觀測值方差分析軟件就可以用來估計非正態分布模型，如重復調用廣義最小二乘來估計廣義線性模型。因此這種方法很流行。

Breslow and Clayton(1993)同時對兩種方法進行系統介紹和比較，為把得分方程組的求解轉化為迭代方式求解提供了理論支持，具體求解使用的方法是迭代（再）加權最小二乘算法（Nelder and Wedderburn,1972）。Vonesh等（2002）同時擴展了罰擬似然方法和廣義估計方程，在估計方程中引入了一二階條件矩。Lin等（1997）進一步把協變量效應引入到方差成分中。

Wolfinger and O’Connell（1993）使用的偽似然方法在實際編程應用中有重要地位。偽似然（Carroll and Ruppert，1988,P71）是指“偽稱其他回歸參數β已知并且等于當前的估計值進一步基于正態性假定得到參數θ的極大似然估計。”偽似然實施過程中，是輪流估計參數β和θ，迭代直到收斂。如果在估計隨機效應協方差矩陣G和殘差協方差矩陣R的參數時使用了ML方法，那就得到PL估計量;如果在估計隨機效應協方差矩陣G和殘差協方差矩陣R的參數時使用了REML方法，那就得到REPL估計量。偽似然方法提供了SS推斷和PA推斷的統一的框架，并且把PQL和MQL作為特例包含進來。而且，該方法可以允許自定義G和R的協方差結構以適應更一般的模型。更重要的一點是，PL或REPL方法可以很方便地利用SAS的固有線性混合模型模塊來實現。

Davidian and Giltinan(1993)討論了兩類非線性混合模型的估計方法，一種他們稱為pooled two stage（PTS）法，另外一種他們稱為線性化混合效應（LME）方法。PTS法對從每個個體得到的估計進行合成；LME方法是Vonesh and Carter（1992）線性化方法的推廣。Wolfinger and O’Connell（1993）偽似然方法與LME方法非常相似，不同之處在于，LME對于隨機效應的線性化是0，但Wolfinger and O’Connell（1993）的方法以及Lindstrom and Bates（1990）的方法都對參數的當前估計進行線性化。

總結一下，Breslow and Clayton(1993)利用擬似然總結了PQL和MQL兩種方法，但最終以特例形式在Wolfinger and O’Connell(1993)的偽似然方法中體現出來。同時，偽似然方法與LME方法思路非常相似，只是在線性近似處理上有所不同。

綜合這些模型，得出本模型參數估計的基本思路如下：

首先固定τ，偽稱其估計已知τ=τ?，然后討論其余參數的估計，最后基于搜索算法求的使所有參數估計最優的τ?。

具體說來，通過使用Laplace近似的方法（Tierney等1989）或者Taylor序列近似，可以得到目標似然函數的線性近似或者分析近似，構造出得分方程組，求解都可以運用Fisher得分法得到迭代（再）加權最小二乘估計。也即可以通過下面的方程組以迭代方式求解β和u：

然后可以得到隨機效應的收縮估計：

迭代（再）加權最小二乘法迭代求解的時候要求G?非奇異。但是在迭代運算的過程中，常常會出現迭代估計失敗造成方差成分取邊界值0。如果奇異，需要使用Cholesky分解進行調整(Henderson 1984)。

設是?的Cholesky下三角根那么迭代（再）加權最小二乘法方程組變為：

而且，這個結論還可以推廣，如果把V用一致性估計替換，以上漸進性依然成立。

把這些參數估計帶入到目標似然函數，得到調整后參數?的近似概括擬似然函數，利用這個公式關于?偏導并令其為0，得到得分方程，基于此可以求解

在建模過程中假定尺度參數φ=1可以得到Breslow and Clayton的罰擬似然方法（penalized quasi-likelihood,PQL）；而 Wolfinger and O’Connell方法為偽似然方法（pseudo-likelihood,PL）或者受限偽似然方法（restricted pseudo-likelihood,REPL），在建模過程中假定尺度參數φ是未知的，因此利用PL得到尺度參數的極大似然估計或者利用REPL方法獲得尺度參數的受限極大似然估計。從這種意義上說，罰擬似然方法只是偽似然方法中尺度參數φ=1的特例。

利用前面的方法可以得出固定τ下的參數估計，假定為，把它們帶回聯合似然函數，得到聯合似然函數的估計式

因為y*=yτ，在τ＞0時為單調函數，那么τ的估計可以直接對聯合似然函數極大化來求解：

本文設計非對稱有記憶二分搜索算法：

圖1 直接搜索算法流程圖

縱觀算法，運算不外乎兩類：一是τ的變化，二是步長d的變化。我們把τ的變化（加減d個步長）定義為延展，把τ不變化而步長改變定義為收縮。一旦運算進入收縮階段，就會收斂到最小值，當然，這個最小值可能是局部最小值。

初始化時可以選擇τ=1，初始化步長d，d的取值不易大于1，防止τ取負，同時，d的取值不易太小，以防解收斂于局部最小值；該算法對往正方向移動還是往負方向移動是非對稱的。k為步率，反映了調整步長的程度。當τ的解往正方向或負方向移動時，步長不改變，恒定為；但是，一旦發現無論是正方向移動還是負方向移動都失敗時（極值點出現在當前估計值的(τ-d,τ+d)區間上），程序會自動調整步長為原來步長的k倍，以防止出現來回震蕩的情況。

3 結論

本文提出了廣義Gamma分布簇廣義線性混合模型的參數估計方法，該方法充分利用了指數分布簇廣義線性混合模型參數估計的結果，并且參數估計表達式上與6種廣義線性混合模型參數估計結果形式相同，但是模型假設和表示的含義有異。而且，該方法可以通過反復調用廣義線性混合模型中的參數估計模塊來實現，具有漸進正態性，因此有較好的易用性。

[1]謝遠濤，楊娟.廣義Gamma分布簇廣義線性混合模型的構建[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).