貝葉斯偏正態(tài)線型混合模型的應(yīng)用

王明高

（1.山東工商學(xué)院 統(tǒng)計(jì)學(xué)院，山東 煙臺(tái) 264005；2.中國(guó)人民大學(xué) 統(tǒng)計(jì)學(xué)院，北京 100872）

0 引言

線性混合模型經(jīng)常用于分析重復(fù)觀測(cè)的縱向數(shù)據(jù)，在很多領(lǐng)域如農(nóng)業(yè)、生物、衛(wèi)生、經(jīng)濟(jì)、保險(xiǎn)精算等都有廣泛應(yīng)用。線性混合模型在線性解釋變量中包含隨機(jī)效應(yīng)，隨機(jī)效應(yīng)不但可以考查同類觀測(cè)值之間的相關(guān)性結(jié)構(gòu)，還可以考慮不同類別的異質(zhì)性。該模型具有廣泛的應(yīng)用，主要是他們建模的靈活性，可以處理對(duì)稱數(shù)據(jù)和非對(duì)稱數(shù)據(jù)，而且存在比較有效的軟件來(lái)擬合線性混合模型。

線性混合模型是線性模型的擴(kuò)展，假設(shè)一個(gè)具有m個(gè)類別的數(shù)據(jù)，每類數(shù)據(jù)具有ni個(gè)個(gè)體i=1，…，m，假設(shè)Yi是一個(gè)對(duì)樣本i連續(xù)性測(cè)量的ni維向量。則Yi所具有的線性混合模型一般表達(dá)式為：

雖然上面討論的模型對(duì)縱向數(shù)據(jù)建模具有很大的靈活性，但是對(duì)于偏離假設(shè)分布情況，他們?nèi)狈ο鄳?yīng)的穩(wěn)健性。假設(shè)隨機(jī)效應(yīng)和響應(yīng)變量(隨機(jī)誤差)都服從正態(tài)分布并不適合所有數(shù)據(jù)，比如響應(yīng)變量為計(jì)數(shù)數(shù)據(jù)、二元數(shù)據(jù)和偏態(tài)數(shù)據(jù)等，并且隨機(jī)效應(yīng)不能保證服從對(duì)稱的正態(tài)分布。對(duì)標(biāo)準(zhǔn)線性混合模型的擴(kuò)展在一些文獻(xiàn)中涉及，主要是用不同的分布來(lái)代替隨機(jī)效應(yīng)的正態(tài)分布假設(shè)，如Zhang(2001)[1]Ghidey(2004)[2]等。這些文章中有一個(gè)共同點(diǎn)就是對(duì)誤差都假設(shè)服從正態(tài)分布。在涉及連續(xù)性測(cè)量時(shí)，對(duì)個(gè)體之間的分布假設(shè)為正態(tài)分布在大多數(shù)的情況下比較合理，在其他一些情況下需要進(jìn)行轉(zhuǎn)換來(lái)滿足正態(tài)分布假設(shè)，比如對(duì)數(shù)據(jù)進(jìn)行對(duì)數(shù)變換。雖然這些方法可以給出比較合適的實(shí)證結(jié)果，對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)換可能使原來(lái)數(shù)據(jù)潛在生成機(jī)制提供的信息減少，另外數(shù)據(jù)轉(zhuǎn)換適用于單個(gè)個(gè)體，不能保證聯(lián)合正態(tài)性，還有轉(zhuǎn)換后的變量比較難以解釋。所以對(duì)數(shù)據(jù)進(jìn)行建模時(shí)要根據(jù)數(shù)據(jù)的分布形式，選擇合適的模型進(jìn)行評(píng)估，如果存在更加合適的理論模型我們盡量避免對(duì)數(shù)據(jù)進(jìn)行變換。

基于以上原因，本文探討對(duì)線性混合模型隨機(jī)效應(yīng)和誤差效應(yīng)的非正態(tài)假設(shè)。介紹更加靈活的分布處理非正態(tài)分布的這種情況，而不去進(jìn)行特別的數(shù)據(jù)變換。

1 偏正態(tài)分布

本文主要研究偏正態(tài)分布，偏正態(tài)分布在很多文獻(xiàn)中都有涉及，在這里介紹Arellano-Valle(2007)[3]、Arellano-Valle(2005)[4]和 Sahu(2003)[5]、Azzalini(1999)[6]、Azzalini(1996)[7]中所介紹的偏正態(tài)分布及其相應(yīng)的特征。

1.1 一元偏正態(tài)分布

假設(shè)X是一個(gè)連續(xù)隨機(jī)變量，服從偏正態(tài)分布則概率密度函數(shù)是：

偏正態(tài)分布的一個(gè)重要性質(zhì)就是各階矩都存在，矩母函數(shù)如下所示：

偏正態(tài)分布可以認(rèn)為是對(duì)正態(tài)分布的擴(kuò)展，一元偏正態(tài)分布的曲線如圖1所示，該圖顯示δ（σ=1，μ=0）取五個(gè)不同值的圖形分布情況，當(dāng)δ=0時(shí)曲線表示對(duì)稱分布標(biāo)準(zhǔn)正態(tài)分布，其他δ值的曲線呈現(xiàn)偏態(tài)分布。

圖1 偏正態(tài)分布曲線圖

從上面公式可以看出，偏正態(tài)分布的表達(dá)式比較復(fù)雜，這就給模型的建模帶來(lái)困難，但是通過(guò)對(duì)模型進(jìn)行轉(zhuǎn)化，可以使該分布通過(guò)兩個(gè)相互獨(dú)立的正態(tài)隨機(jī)變量來(lái)獲得偏正態(tài)分布隨機(jī)變量。Henze(1986)[8]一文對(duì)一元偏正態(tài)分布的隨機(jī)化表述做了一些研究，偏正態(tài)分布形式也可以寫(xiě)成如下等價(jià)形式：

這里U和V是相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)隨機(jī)變量，由和所表示的等價(jià)形式簡(jiǎn)化了偏正態(tài)分布函數(shù)，便于對(duì)偏正態(tài)分布混合模型建模。

1.2 多元偏正態(tài)分布

2 偏正態(tài)線性混合模型

2.1 偏正態(tài)線性混合模型簡(jiǎn)介

本模型主要目的是估計(jì)參數(shù)本文是在貝葉斯框架內(nèi)對(duì)線性混合模型進(jìn)行推斷分析，這種模型的推斷是建立在響應(yīng)變量Y的分布基礎(chǔ)之上，該線性混合模型的隨機(jī)效應(yīng)和誤差項(xiàng)可以服從一元或多元偏正態(tài)分布，這個(gè)模型的一個(gè)重要優(yōu)勢(shì)就是可以用來(lái)評(píng)估偏態(tài)非對(duì)稱數(shù)據(jù)，隨機(jī)效應(yīng)并不局限于服從對(duì)稱的正態(tài)分布。

2.2 貝葉斯偏正態(tài)線性混合模型

要寫(xiě)出各參數(shù)的條件分布。這種算法在開(kāi)始的時(shí)候需要給出所有隨機(jī)變量的初始值，再由條件分布生成樣本，直到樣本收斂。我們可以用Gelman(1992)[11]定義的統(tǒng)計(jì)量檢驗(yàn)所生成的樣本是否收斂。

3 應(yīng)用

3.1 雇主責(zé)任險(xiǎn)的案例

本文引用一個(gè)雇主責(zé)任險(xiǎn)的案例，該險(xiǎn)種主要處理永久性或部分喪失工作能力的索賠，在這個(gè)縱向數(shù)據(jù)中共有121個(gè)職業(yè)分類即風(fēng)險(xiǎn)類別，在7年觀察期內(nèi)發(fā)生的847條數(shù)據(jù)。這些數(shù)據(jù)包括職業(yè)類別、觀察年、收入和賠付額。在Frees(2001)[12]的論文中分析了這些數(shù)據(jù)，他選用純保費(fèi)(PP=Loss/Payroll)為響應(yīng)變量即賠付額與收入的比率，解釋變量有觀察年、職業(yè)類別。Frees(2001)[14]在建立模型時(shí)沒(méi)有直接用純保費(fèi)數(shù)據(jù)，而是使用純保費(fèi)的對(duì)數(shù)數(shù)據(jù)。從圖2的左側(cè)純保費(fèi)頻率分布圖可以看出，純保費(fèi)是一個(gè)右偏的分布圖形。將純保費(fèi)數(shù)據(jù)進(jìn)行對(duì)數(shù)變換，對(duì)數(shù)變換后的數(shù)據(jù)分布如圖2右側(cè)所示，該圖呈現(xiàn)出類似正態(tài)分布的對(duì)稱分布圖形。

圖2 純保費(fèi)直方圖

這些變量的取值在Frees(2001)[12]文章的第四部分中的表3有描述，該表顯示觀察期內(nèi)所有職業(yè)的平均收入為195百萬(wàn)美元，平均純保費(fèi)為0.0203。從表中可以看出平均保費(fèi)沒(méi)有呈現(xiàn)隨觀察年變化的趨勢(shì)，收入表現(xiàn)出隨觀察年增加的趨勢(shì)。另外部分職業(yè)的純保費(fèi)隨觀察年的變化如圖3所示，我們可以看出純保費(fèi)并不隨觀察年變化，但是純保費(fèi)對(duì)數(shù)值的波動(dòng)情況跟純保費(fèi)的波動(dòng)情況不同，純保費(fèi)隨觀察年的變化偏向一側(cè)，純保費(fèi)對(duì)數(shù)值呈現(xiàn)對(duì)稱波動(dòng)。

圖3 部分職業(yè)的純保費(fèi)隨觀察年的變化趨勢(shì)

上面這種表示方式比較重要，因?yàn)橛欣谟肂UGS編寫(xiě)程序。應(yīng)用Gibbs抽樣需要得到各參數(shù)的條件分布，根據(jù)各參數(shù)的條件分布進(jìn)行抽樣。但是本文應(yīng)用Open-BUGS軟件，只需要根據(jù)模型的分層結(jié)構(gòu)進(jìn)行編程，沒(méi)有必

雖然純保費(fèi)跟觀察年之間沒(méi)有明顯的變化趨勢(shì)，但是該數(shù)據(jù)是縱向數(shù)據(jù)，需要考慮同一組數(shù)據(jù)之間的相關(guān)性。在該數(shù)據(jù)中同一職業(yè)類別不同觀察年之間的觀察值具有顯著地相關(guān)性，在建立模型時(shí)需要考慮這些問(wèn)題。

3.2 態(tài)線性混合模型的建立

本文應(yīng)用偏正態(tài)分布定義隨機(jī)效應(yīng)bi和隨機(jī)誤差項(xiàng)eit構(gòu)成偏正態(tài)線性混合模型，本文將通過(guò)三個(gè)模型進(jìn)行比較分析。第一個(gè)模型是一般的線性混合模型，響應(yīng)變量為純保費(fèi)的對(duì)數(shù)值，隨機(jī)效應(yīng)和隨機(jī)誤差項(xiàng)服從正態(tài)分布。第二個(gè)模型引入偏正態(tài)分布，假設(shè)隨機(jī)誤差項(xiàng)服從偏正態(tài)分布，隨機(jī)效應(yīng)服從正態(tài)分布，響應(yīng)變量為純保費(fèi)。第三個(gè)模型的隨機(jī)誤差項(xiàng)和隨機(jī)效應(yīng)都服從偏正態(tài)分布，響應(yīng)變量為純保費(fèi)。

表1 三個(gè)模型參數(shù)的后驗(yàn)分布得出均值與標(biāo)準(zhǔn)差

為了選取合適的模型我們需要對(duì)模型進(jìn)行比較分析，貝葉斯模型的評(píng)估結(jié)果用DIC統(tǒng)計(jì)量Spiegelhalter(2002)[13]進(jìn)行比較分析，模型DIC表達(dá)式如下所示：

另外對(duì)模型的殘差考查如圖4、圖5所示，這兩個(gè)圖只給出模型一和模型二的殘差圖，其中左邊的a圖是模型一的殘差圖，右邊的b圖是模型二的殘差圖。從這兩圖中可以看出模型一的殘差具有一定的異方差性，殘差波動(dòng)范圍比較大，其變化范圍超過(guò)了模型響應(yīng)變量純保費(fèi)的取值范圍。模型二的殘差波動(dòng)范圍比較小，而且殘差主要聚集在零值附近，說(shuō)明該模型二很好的擬合了該數(shù)據(jù)。另外圖6表示純保費(fèi)與其模擬值之間的散點(diǎn)圖，可以看出由b圖所表示的模型二擬合效果明顯好于a圖所表示的模型一，模型二給出的模擬值能夠很好地近似實(shí)際的純保費(fèi)數(shù)據(jù)。從下面的圖形我們可以看出假設(shè)隨機(jī)誤差服從偏正態(tài)分布的模型能夠更好的反映實(shí)際數(shù)據(jù)。

圖4 純保費(fèi)對(duì)殘差的散點(diǎn)圖

圖5 損失對(duì)殘差的散點(diǎn)圖

圖6 純保費(fèi)對(duì)其模擬值的散點(diǎn)圖

4 結(jié)論

本文主要研究如何更好對(duì)具有非對(duì)稱分布特征的縱向數(shù)據(jù)進(jìn)行建模，傳統(tǒng)的方法一般對(duì)數(shù)據(jù)進(jìn)行變換，比如對(duì)數(shù)變換等。但是對(duì)數(shù)據(jù)進(jìn)行變換可能損失一些信息，本文探討引入適合所分析數(shù)據(jù)特征的分布來(lái)處理非對(duì)稱數(shù)據(jù)。

在文中引用一個(gè)具有右偏特征的保險(xiǎn)損失數(shù)據(jù)，在Frees(2001)文中通過(guò)對(duì)數(shù)變換來(lái)處理該數(shù)據(jù)的右偏性，應(yīng)用線性混合模型來(lái)分析，本文根據(jù)所選數(shù)據(jù)的分布特征，選用可以處理偏態(tài)數(shù)據(jù)的偏正態(tài)分布，建立偏正態(tài)線性混合模型，并且與對(duì)數(shù)變換的線性混合模型進(jìn)行比較，評(píng)估結(jié)果得出偏正態(tài)線性混合能夠更好擬合該數(shù)據(jù)。通過(guò)本文的一個(gè)實(shí)例數(shù)據(jù)分析發(fā)現(xiàn)，對(duì)一些非對(duì)稱數(shù)據(jù)的建模，應(yīng)用非對(duì)稱模型更加有效。實(shí)際上一般線性混合模型假設(shè)隨機(jī)效應(yīng)和隨機(jī)誤差服從正態(tài)分布，可以認(rèn)為是對(duì)偏正態(tài)分布線性混合模型的特殊情況。

雖然所涉及的線性混合模型比較復(fù)雜，本文應(yīng)用貝葉斯蒙特卡羅方法Gibbs抽樣對(duì)偏正態(tài)線性混合模型的參數(shù)進(jìn)行估計(jì)，這種方法可以通過(guò)免費(fèi)軟件OpenBUGS非常方便的得出相關(guān)結(jié)果。本模型估計(jì)一個(gè)比較大的優(yōu)勢(shì)就是可以對(duì)模型的隨機(jī)效應(yīng)和隨機(jī)誤差給出更加靈活的假設(shè)。通過(guò)本文的分析可以發(fā)現(xiàn)，基于貝葉斯方法的偏正態(tài)線性混合模型處理縱向數(shù)據(jù)，使模型不再局限于傳統(tǒng)的假設(shè)，可以建立符合實(shí)際數(shù)據(jù)特征的混合模型，從而能夠更好的處理所研究數(shù)據(jù)。

[1]Zhang D,Davidian M.Liner Mixed Models With Flexible Distributions of Random Effects for Longitudinal Data[J].Biometrics,2001,（57）.

[2]Ghidey W,Lesaffre E,Eilers,P.Smooth Random Effects Distribution in A Linear Mixed Model[J].Biometrics,2004,（60）.

[3]Arellano-Valle R B,Bolfarine H,Lachos V H.Bayesian Inference for Skew-Normal Linear Mixed Models[J].Journal of Applied Statistics,2007,（34）.

[4]Arellano-Valle R B,Genton M G.Fundamental Skew Distributions[J].Journal of Multivariate Analysis,2005,（96）.

[5]Sahu S K,Dey D K,Branco M.A New Class of Multivariate Distributions With Applications to Bayesian Regression Models[J].The Canadian Journal of Statistics,2003,（29）.

[6]Azzalini A,Capitanio A.Statistical Applications of The Multivariate Skew-Normal Distributions[J].Journal of The Royal Statistical Society,1999,（61）.

[7]Azzalini A,Dalla-Valle A.The Multivariate Skew-Normal Distribution[J].Biometrika,1996,（3）.

[8]Henze,N.A Probabilistic Representation of The Skew-Normal Dis-Tribution[J].Scand J Statist,1986,（13）.

[9]Arellano-Valle R B,Azzalini,A.On The Unification of Families of Skew-Normal Distributions[J].Scand J Statist,2006,（33）.

[10]Azzalini,A.The Skew-Normal Distribution and Related Multivariate Families[J].Scandinavian Journal of Statistics,2005,32(2).

[11]Gelman A,Rubin D.Inference From Iterative Simulation Using Multiple Sequences[J].Statistical Science,1992,（7）.

[12]Frees E W,Young V R,Luo Y.Case Studies Using Panel Data Models[J].North American Actuarial Journal 2001,5(4).

[13]Spiegelhalter D J,Best N G,Carlin B P,Van Der Linde,A.Bayesian Measures of Model Complexity and Fit(With Discussion)[J].Journal of The Royal Statistical(B),2002,（64）.