基于廣義線性混合模型的參數(shù)估計

吳琳

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

關鍵詞：廣義線性混合模型;方差分量;參數(shù)估計;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

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

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

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

1 模型介紹和分析

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

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

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

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

2 數(shù)據(jù)集介紹和描述性分析

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

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

3 參數(shù)估計

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

3.1 PQL參數(shù)估計法

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

，

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

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

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

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

3.2 貝葉斯參數(shù)估計

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

參考文獻

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

[2] 陳永成，STEPHEN W.RAUDENBUSH.用蒙特卡羅方法求取廣義線性混合模型之最大似然估計：應用于求取乳癌死亡率之小區(qū)域估計[J].應用概率統(tǒng)計，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.