王明高

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

0 引言

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

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

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

基于以上原因，本文探討對線性混合模型隨機效應(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是一個連續(xù)隨機變量，服從偏正態(tài)分布則概率密度函數(shù)是：

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

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

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

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

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

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

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

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

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

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

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

3 應(yīng)用

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

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

圖2 純保費直方圖

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

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

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

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

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

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

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

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

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

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

圖5 損失對殘差的散點圖

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

4 結(jié)論

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

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

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