吳琳

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

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

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

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

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

1 模型介紹和分析

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

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

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

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

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

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

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

3 參數(shù)估計(jì)

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

3.1 PQL參數(shù)估計(jì)法

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

，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

參考文獻(xiàn)

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

[2] 陳永成，STEPHEN W.RAUDENBUSH.用蒙特卡羅方法求取廣義線性混合模型之最大似然估計(jì)：應(yīng)用于求取乳癌死亡率之小區(qū)域估計(jì)[J].應(yīng)用概率統(tǒng)計(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].昆明理工大學(xué)學(xué)報(bào)（理工版），2007，32（4）：107-113.