基于廣義線性混合模型的參數估計

吳琳

摘要：本文主要討論了廣義線性混合模型（即GLMM）模型的兩種不同的參數估計方法，基于懲罰擬似然參數估計方法和MCMC貝葉斯參數估計方法。本文考慮借助開源軟件R對癲癇病數據集epil進行統計推斷，并與其他模型和參數估計方法進行比較。

關鍵詞：廣義線性混合模型;方差分量;參數估計;glmmPQL;MCMCglmm

Abstract： This paper mainly discusses two different parameter estimation methods of GLMM model， one is the penalty quasi likelihood parameter estimation method and the other is the MCMC parameter estimation method based on Bayesian， and this paper considers the use of open source software R to estimate the parameters of a specific data set （epil） Compared with other models and parameter estimation methods， it is concluded that the application of GLMM is very extensive and necessary.

Key Words：GLMM;variance component;parameter estimation;glmmPQL;MCMCglmm

在對實際問題建模時，當我們感興趣的變量是一個正態分布隨機變量時，通常考慮使用線性模型（LM），要求響應變量滿足正態性和獨立性，而實際數據會違背相互獨立的假設，這些數據之間往往存在某種相關結構，由此發展了線性混合模型（LMM），在LMM中，放松了獨立性要求，考慮響應變量之間的相關性。在后續的研究中發現除了滿足正態分布的隨機變量之外，還存在很多非正態分布，如泊松分布，二項分布等。對于這類分布，考慮對其均值進行建模，發展成廣義線性模型（GLM），相比LM，GLM放松了要求響應變量的正態性，但是仍然需要滿足獨立性，但實際中數據之間的相關性不可忽略。忽略相關性將會導致偏差或者得出錯誤的結論。為了體現出數據之間的相關性，在廣義線性模型的連接函數上加上了一個隨機變量，也就是本文所研究的廣義線性混合模型（GLMM）[3]。理論上，廣義線性混合模型是廣義線性模型和線性模型的結合，同時它也是線性混合模型的拓展。廣義線性混合模型適合用來模擬存在相關性的區組效應的數據或縱向數據，此時不再要求響應變量滿足正態性和獨立性，模型設定更加靈活，更能反映出響應變量之間的相關性，這是廣義線性混合模型的一個顯著優點。

對于GLMM的參數估計，以往很多學者提出各種方法。針對GLM的參數估計首先構造似然函數，對隨機效應進行處理，無論是對條件似然函數的積分處理還是對邊際似然函數的極大似然估計[2]的過程中都存在很大的困難，尤其是在對隨機效應的積分處理往往伴隨各種高維積分。后來的研究中，有學者提出擬似然（quasi—likelihood）函數[1]，詳細參考McCullagh[3]，擬似然的提出這類模型提供了推斷的基礎，后來的學者提出了基于懲罰你似然函數的參數估計方法。對于多數非正態模型，使用Monte-Carlo積分方法[2]。除此之外，還可以在先驗分布[4]中使用重要性抽樣或者Gibbs抽樣避免數值積分，而對于極大似然估計方法最常用的是EM算法和MCEM算法。

參數估計方法比較豐富，但目前針對廣義線性模型的參數估計方法主要是基于擬似然或者懲罰擬似然角度出發的極大似然估計，而從貝葉斯后驗角度出發進行的參數估計目前還并不多，一般進行貝葉斯參數估計借助于WinBUGS或OpenBUGS等軟件進行，遇到較大數據集的情況時運行速度較慢，本文借助于R軟件中的MCMCglmm包中的MCMCglmm函數[5]可以很好的避免運行速度過慢的缺陷，同時可以達到一個很好的參數估計的效果。本文主要針對R中自帶的數據集epil進行懲罰擬似然參數估計和馬爾可夫蒙特卡洛參數估計兩種估計方法的結果比較。

1 模型介紹和分析

廣義線性混合模型（GLMM）是廣義線性模型（GLM）的拓展，GLM需要包含三個部分： （1）響應變量的分布 ;（2）線性預測器 ;（3）連接函數 。模型可表示為：

其中， 為Y的分布函數， 為Y的期望， 是單調可導的連接函數， 為線性預測器。而對廣義線性模型中參數 進行參數估計時比較常用的一種方法是迭代加權最小二乘（IWLS）算法，具體算法的細節可以參考McCullagh 和 Nelder（1978）的工作。

而廣義線性混合模型（GLMM）與廣義線性模型相似，需要以下四個組成部分，除了GLM的三部分還需要添加一個隨機效應： ，其中 為 的方差;更具體的說： ，其中 ，且同時有 。令 ，則可以得到樣本的似然函數：

其中 表示給定隨機項 時Y的條件概率密度函數， 為b的概率密度函數，感興趣的是模型中參數 的估計，即估計出固定效應 和方差分量 。

2 數據集介紹和描述性分析

本文所考慮使用的數據集為R中MASS包中的癲癇病數據集epil，由Thall和Vail（1990）年給出的59個癲癇病患者的2周內的發病次數，患者被隨機分配到兩個不同的實驗組，一個組進行新藥治療（Trt=1），一個組進行對照的安慰劑治療（Trt=0），具體介紹可在R中進行查看。

對數據集epil進行描述性分析，新藥治療組的樣本相關性比安慰劑組的樣本相關性更大，同時從藥物治療組和安慰劑對照組不同時期發病次數箱線圖（圖1）無法斷定癲癇犯病次數隨時期的變化產生而產生明顯的變化，但是可以明顯地看出新藥治療組的發病次數比安慰劑對照組的發病次數水平更低。

3 參數估計

接下來考慮兩種不同的參數估計方法來對廣義線性混合模型（GLMM）的參數進行估計，基于懲罰擬似然方法（PQL）的極大似然估計方法，具體的細節可參考干曉蓉（2007）的論文工作[6];基于貝葉斯后驗分布出發的參數估計方法，具體的細節可以參考Jarrod Hadfield的工作。

3.1 PQL參數估計法

擬似然方法并不要求響應變量是某一個具體已知的分布，只需知道響應變量的均值和方差，在大樣本情形下可近似為正態分布。然而為了減少方差估計的無偏性，采用添加懲罰項的方法來提高估計的精度。假設已知響應變量 的條件分布的均值和方差為：

，

連接函數 的逆函數為 ，即： ，其中 ，有 ， ，有 ，其中I為對角陣，對角線元素為單位陣， 。得到擬似然函數：

借助于拉普拉斯近似，進一步可以得到懲罰擬似然函數為：

針對上式進行求導可以得到方程并進一步得到參數的極大似然估計量。接下來考慮在癲癇病數據集上進行實例應用。我們考慮在epil數據集上應用GLMM模型通過PQL方法來估計參數，首先考慮建立模型：

利用R軟件中的glmmPQL進行參數估計，考慮兩種情況，將period的四個階段其為數值型變量和因子型變量。

3.2 貝葉斯參數估計

貝葉斯統計推斷是從后驗分布出發，通常對未知參數有一定的先驗信息或無信息先驗，當給定先驗先驗信息，可求后驗分布，進而完成貝葉斯參數計。

參考Jarrod Hadfield的工作，考慮第i個數據的概率為： ，假設 服從泊松分布，則有：

其中 為泊松分布規范參數。解釋變量之間的線性有下列線性關系：

其中X為與固定效應有關的設計矩陣，Z是與隨機效應有關的設計矩陣，e為殘差，假設 來自多元正態分布，滿足：

利用MCMC抽樣方法我們可以得到：

按照 更新參數 ， ，其中C為固定模型的系數矩陣，滿足：

其中 ，B為固定效應的先驗方差， 和 為從多元正態分布中隨機抽取的樣本：

且 ， 即為從下列分布中抽樣的樣本： 。

可借助Jarrod Hadfield[5]編寫的MCMCglmm包中的MCMCglmm函數實現廣義線性混合模型的參數估計，由于對癲癇病數據集缺乏先驗信息，故不設置先驗信息，或者設置一個無信息先驗信息。利用MCMCglmm函數對數據集進行擬合，考慮每10個間隔進行一次抽樣，共抽取13000個樣本，除去預燒期的3000個抽樣樣本，進行貝葉斯參數估計，同樣考慮period的兩種情形，模型的DIC值分別為1153.947和1158.283，兩種情形下隨機效應的參數估計結果為：0.252和0.249.

根據表1和表2中的glmmPQL參數估計結果得出結論：新藥的回歸系數是顯著的，年齡并不顯著，基礎發病次數是顯著的，但基礎發病次數與新藥的交叉作用不顯著，同時從整體上看階段是一個顯著的回歸項，具體到四個階段的每一階段，只有第四階段的作用是顯著的，其它的回歸項的顯著性和估計結果基本相似。同時兩種不同情形下的隨機效應估計為0.443和0.444，兩者相差不大。這與前面的數據描述性分析是相互對照的。

表1和表2中同樣給出了貝葉斯框架下的MCMCglmm參數估計結果，當將period視為數值型變量時，此時隨機效應的估計值為 ：0.252，可信區間為（0.129，0.402）。當將period視為因子型變量，此時隨機效應的估計值為 ：0.250，可信區間為（0.131，0.409）。

看出無論是將period看成數值型變量還是因子型變量時，各個時期都不顯著，這與PQL方法下的結果是不同的。同時在MCMCglmm包中還可以給出隨機效應的迭代軌跡圖像以及隨機效應的概率密度函數圖像，可以通過迭代軌跡圖判斷是否收斂以及基于后驗分布的統計推斷是否合理。

4 結論

文章主要考慮了兩種廣義線性混合模型的參數估計方法，PQL方法下的廣義線性模型方差分量比MCMC方法下的要略微偏大一些，但兩種方法下大部分的回歸結果是相似的，這也再次證明了這兩種方法的參數估計結果是比較準確的，同時第二種方法下的貝葉斯參數估計方法的效率較高，運行時間比利用WinBUGS或OpenBUGS等軟件的效率要高很多，這也是這個方法的一個顯著的優勢。同時本文也將一般的參數估計方法上拓寬了，將貝葉斯的參數估計方法應用過來，且參數估計的結果表明二者估計的結果是相似且準確的。

參考文獻

[1] 曹紅艷，劉桂芬，曾平，等.預測性偽似然法和貝葉斯法廣義線性混合模型估計[J].中國藥物與臨床，2008，8（011）：861-863.

[2] 陳永成，STEPHEN W.RAUDENBUSH.用蒙特卡羅方法求取廣義線性混合模型之最大似然估計：應用于求取乳癌死亡率之小區域估計[J].應用概率統計，2006，22（1）：69-80.

[3] McCulloch CE，Searle SR.Generalized，linear and mixed models.New York：Wiley，2001.

[4] Gelman A.Prior distribution for variance parameters in hierar-

chical models. Bayesian Analysis， 2006， 1（3）：515-534.

[5] Hadfield J D . MCMC Methods for Multi-Response Generalized Linear Mixed Models： The MCMCglmm R Package[J]. Journal of Statistical Software， 2010， 33（02）：1-22.

[6] 干曉蓉.廣義線性混合模型[J].昆明理工大學學報（理工版），2007，32（4）：107-113.