貝葉斯偏正態線型混合模型的應用

王明高

（1.山東工商學院 統計學院，山東 煙臺 264005；2.中國人民大學 統計學院，北京 100872）

0 引言

線性混合模型經常用于分析重復觀測的縱向數據，在很多領域如農業、生物、衛生、經濟、保險精算等都有廣泛應用。線性混合模型在線性解釋變量中包含隨機效應，隨機效應不但可以考查同類觀測值之間的相關性結構，還可以考慮不同類別的異質性。該模型具有廣泛的應用，主要是他們建模的靈活性，可以處理對稱數據和非對稱數據，而且存在比較有效的軟件來擬合線性混合模型。

線性混合模型是線性模型的擴展，假設一個具有m個類別的數據，每類數據具有ni個個體i=1，…，m，假設Yi是一個對樣本i連續性測量的ni維向量。則Yi所具有的線性混合模型一般表達式為：

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

基于以上原因，本文探討對線性混合模型隨機效應和誤差效應的非正態假設。介紹更加靈活的分布處理非正態分布的這種情況，而不去進行特別的數據變換。

1 偏正態分布

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

1.1 一元偏正態分布

假設X是一個連續隨機變量，服從偏正態分布則概率密度函數是：

偏正態分布的一個重要性質就是各階矩都存在，矩母函數如下所示：

偏正態分布可以認為是對正態分布的擴展，一元偏正態分布的曲線如圖1所示，該圖顯示δ（σ=1，μ=0）取五個不同值的圖形分布情況，當δ=0時曲線表示對稱分布標準正態分布，其他δ值的曲線呈現偏態分布。

圖1 偏正態分布曲線圖

從上面公式可以看出，偏正態分布的表達式比較復雜，這就給模型的建模帶來困難，但是通過對模型進行轉化，可以使該分布通過兩個相互獨立的正態隨機變量來獲得偏正態分布隨機變量。Henze(1986)[8]一文對一元偏正態分布的隨機化表述做了一些研究，偏正態分布形式也可以寫成如下等價形式：

這里U和V是相互獨立的標準正態隨機變量，由和所表示的等價形式簡化了偏正態分布函數，便于對偏正態分布混合模型建模。

1.2 多元偏正態分布

2 偏正態線性混合模型

2.1 偏正態線性混合模型簡介

本模型主要目的是估計參數本文是在貝葉斯框架內對線性混合模型進行推斷分析，這種模型的推斷是建立在響應變量Y的分布基礎之上，該線性混合模型的隨機效應和誤差項可以服從一元或多元偏正態分布，這個模型的一個重要優勢就是可以用來評估偏態非對稱數據，隨機效應并不局限于服從對稱的正態分布。

2.2 貝葉斯偏正態線性混合模型

要寫出各參數的條件分布。這種算法在開始的時候需要給出所有隨機變量的初始值，再由條件分布生成樣本，直到樣本收斂。我們可以用Gelman(1992)[11]定義的統計量檢驗所生成的樣本是否收斂。

3 應用

3.1 雇主責任險的案例

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

圖2 純保費直方圖

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

圖3 部分職業的純保費隨觀察年的變化趨勢

上面這種表示方式比較重要，因為有利于用BUGS編寫程序。應用Gibbs抽樣需要得到各參數的條件分布，根據各參數的條件分布進行抽樣。但是本文應用Open-BUGS軟件，只需要根據模型的分層結構進行編程，沒有必

雖然純保費跟觀察年之間沒有明顯的變化趨勢，但是該數據是縱向數據，需要考慮同一組數據之間的相關性。在該數據中同一職業類別不同觀察年之間的觀察值具有顯著地相關性，在建立模型時需要考慮這些問題。

3.2 態線性混合模型的建立

本文應用偏正態分布定義隨機效應bi和隨機誤差項eit構成偏正態線性混合模型，本文將通過三個模型進行比較分析。第一個模型是一般的線性混合模型，響應變量為純保費的對數值，隨機效應和隨機誤差項服從正態分布。第二個模型引入偏正態分布，假設隨機誤差項服從偏正態分布，隨機效應服從正態分布，響應變量為純保費。第三個模型的隨機誤差項和隨機效應都服從偏正態分布，響應變量為純保費。

表1 三個模型參數的后驗分布得出均值與標準差

為了選取合適的模型我們需要對模型進行比較分析，貝葉斯模型的評估結果用DIC統計量Spiegelhalter(2002)[13]進行比較分析，模型DIC表達式如下所示：

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

圖4 純保費對殘差的散點圖

圖5 損失對殘差的散點圖

圖6 純保費對其模擬值的散點圖

4 結論

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

在文中引用一個具有右偏特征的保險損失數據，在Frees(2001)文中通過對數變換來處理該數據的右偏性，應用線性混合模型來分析，本文根據所選數據的分布特征，選用可以處理偏態數據的偏正態分布，建立偏正態線性混合模型，并且與對數變換的線性混合模型進行比較，評估結果得出偏正態線性混合能夠更好擬合該數據。通過本文的一個實例數據分析發現，對一些非對稱數據的建模，應用非對稱模型更加有效。實際上一般線性混合模型假設隨機效應和隨機誤差服從正態分布，可以認為是對偏正態分布線性混合模型的特殊情況。

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

[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）.