貝葉斯非線性混合效應模型及其應用研究

王明高，孟生旺

(1.山東工商學院 統(tǒng)計學院，山東 煙臺 264005； 2.中國人民大學 應用統(tǒng)計科學研究中心，北京100872)

?

【統(tǒng)計理論與方法】

貝葉斯非線性混合效應模型及其應用研究

王明高1,2，孟生旺2

(1.山東工商學院 統(tǒng)計學院，山東 煙臺 264005； 2.中國人民大學 應用統(tǒng)計科學研究中心，北京100872)

由于常用的線性混合效應模型對具有非線性關系的縱向數據建模具有一定的局限性，因此對線性混合效應模型進行擴展，根據變量間的非線性關系建立不同的非線性混合效應模型，并根據因變量的分布特征建立混合分布模型?；谝唤M實際的保險損失數據，建立多項式混合效應模型、截斷多項式混合效應模型和B樣條混合效應模型。研究結果表明，非線性混合效應模型能夠顯著改進對保險損失數據的建模效果，對非壽險費率厘定具有重要參考價值。

混合效應模型；非線性模型；保險損失；貝葉斯

一、引 言

混合效應模型在教育、醫(yī)學、社會和金融保險等領域具有廣泛的應用價值。很多保險數據需要進行分層分析，比如同一地區(qū)的車輛和房屋等保險標的都具有相似的風險特征，而不同地區(qū)的這些保險標的具有一定的差異性，這些保險標的的損失數據都不滿足相互獨立的假設條件。對于不滿足獨立性假設并具有層次性的縱向數據，普通的回歸模型將不再適合，而考慮數據間的相關性和層次性的混合效應模型既含有固定效應參數，又含有隨機效應參數，能夠更好地分析縱向數據?；旌闲Ｐ偷暮诵乃枷胧峭ㄟ^在模型中加入隨機效應，來體現(xiàn)組內數據的相關性和不同組數據的異質性。對于觀察數據比較少的個體，混合效應模型的推斷結果會產生向平均值聚集的傾向，從而得到更加穩(wěn)健的估計結果[1-2]。

線性混合效應模型假設變量之間存在簡單的線性關系。為了描述變量之間的非線性關系，需要借助于非線性混合效應模型，該模型能夠很好地處理同一個體重復觀測數據的相關性以及變量之間的非線性關系。多項式在非線性回歸模型中較常使用，但其缺陷是，當多項式的階數較低時模型的擬合效果不佳，而增加多項式的階數又會產生不穩(wěn)定的結果[3]。Hastie 等提出了廣義可加模型[4]112-113；Wood在廣義可加模型中考慮了連續(xù)變量的平滑效應，使模型的擬合效果更加符合實際[5]45-78。對于分層數據的處理，Ruppert 等探討了廣義可加混合效應模型的應用[6]67-86[7]；Brezger 等應用MCMC方法估計該類模型的參數[8-10]；Silva 等應用貝葉斯混合效應B樣條來分析健康指標在各地區(qū)間的非線性變化[11]；Rigby 等提出了基于位置、尺度和形狀參數的廣義可加模型(GAMLSS)[12]，該模型可以進一步擴展為對因變量的多個參數通過連接函數建立解釋變量的可加模型，并在回歸方程中包含解釋變量的線性效應、連續(xù)變量的非線性平滑效應等[13]223-226。

為了探索非線性混合模型在非壽險精算中的應用，分析一組具有非線性關系的保險損失數據，嘗試應用一般多項式、截斷多項式和B樣條三種非線性表達式建立非線性混合效應模型，這些非線性形式在上面的文獻中都有涉及。但是，根據響應變量的分布特征，本文進一步探討因變量服從混合分布條件下的建模問題，并對混合分布中的每一個分布進行非線性回歸分析。通過對各種模型的比較分析，可以發(fā)現(xiàn)非線性混合效應模型可以有效改善對保險損失的預測效果。

二、非線性混合效應模型

線性混合效應模型假設均值參數的回歸方程中含有線性可加的解釋變量[14]254-256，雖然線性混合效應模型的適用范圍較為廣泛，但是變量之間的非線性關系也比較常見。最簡單的非線性形式是多項式，該形式能夠給出數據整體性的趨勢變化，但是缺乏對局部異常波動的擬合能力。在一些實際應用中，非線性表達式的形式是未知的，需要通過數據進行估計，這種非線性關系稱為非參數形式。在回歸模型中，如果只有部分變量具有非線性效應，這就生成了半參回歸模型[15]。Wu 等討論了如何應用似然方法處理非參數混合效應模型[16]341-342，本文將主要應用貝葉斯方法建立非線性混合效應模型。

假設觀察數據中共有m個不同的觀察對象，每個觀察對象具有n個觀察值{yij,i=1,2,…,m,j=1,2,…,n}，每個觀察值服從下述的指數分布族：

(1)

上述分布的均值為E(yij)=μij=a′(θi)，均值的回歸方程為g(μij)=ηij，其中g(·) 為連接函數。

假設數據中有K個解釋變量(xij1,xij2,…,xijK)，應用非參數的方式進行建模，則均值回歸方程可以表示為：

g(μij)=ηij=Xijα+Zijβi+S1(xij1)+…+SK(xijK)

(2)

其中Xij和Zij表示解釋變量向量，參數α和βi分別表示固定效應和隨機效應，回歸方程的前半部分Xijα+Zijβi為線性方程，后半部分為解釋變量的平滑函數S(xijk)。本文主要研究的平滑函數為多項式、截斷多項式和B樣條。

(一)多項式模型

多項式非線性模型比較簡單，假設某解釋變量xijk的平滑函數S(xijk)為多項式，即S(xijk)=γi0+γi1xijk+…+γilxijkl，這是一個l階多項式，其中回歸參數γ是隨機效應參數，也可以定義成固定效應參數。多項式回歸方程對于波動趨勢不顯著的曲線擬合效果較好，但是對于比較復雜的曲線關系，需要在變量xijk的定義區(qū)間內選取不同的節(jié)點。

(二)截斷多項式模型

(3)

(三)B樣條多項式模型

l階的B樣條多項式有r個節(jié)點，可以表示成r+l+1=K個基函數的線性組合。變量xijk的B樣條平滑函數S(xijk)可以表示為：

(4)

其中B樣條基函數的定義如下：

(5)

當l=1時的基函數為：

(6)

每一個基函數由區(qū)間[κj-1,κj)和[κj,κj+1)上的兩個線性部分組成，在節(jié)點κj處連接。

更高階的B樣條基函數可以通過下面的遞推公式求得：

(7)

(四)懲罰樣條(P樣條)

在截斷多項式回歸模型中，最小化下述的懲罰殘差平方和可以避免出現(xiàn)較大的參數估計值以及過擬合現(xiàn)象：

(8)

平滑系數λ≥0用于控制懲罰項的影響。當λ→0時，懲罰項的影響消失，得到參數γ的最小二乘估計值；當λ→時，參數估計值為，平滑函數f(yij)即為l階多項式。

三、貝葉斯非線性混合效應模型的參數估計與模型評價

(一)貝葉斯非線性混合效應模型的參數估計

非線性回歸模型的參數估計比較復雜，對于上面介紹的懲罰樣條模型除了可以應用懲罰最小二乘方法進行估計外，也可以應用貝葉斯方法對懲罰樣條參數進行估計。下面主要對懲罰樣條多項式回歸模型進行分析。

在貝葉斯方法中，無需單獨列出懲罰項，而是給出樣條參數γ的先驗分布，應用隨機游走的動態(tài)表達式表示差異性懲罰項，一階隨機游走可以表示為：

γj=γj-1+uj,uj～N(0,τ2),j=2,3,…,K

(9)

對于上面的一階差分需要給出起點值γ1的先驗分布，如果沒有任何信息，可以假設γ1的先驗分布為p(γ1)∝const，其中const表示任意常數值。其他樣條參數的先驗分布以條件分布的形式出現(xiàn)：

γj|γj-1,…,γ1～N(γj-1,τ2)

(10)

參數γj的期望值為其滯后值γj-1，誤差項uj的方差為τ2，該方差值越大，說明實際值對期望均值的偏離越大。

(11)

可見，該聯(lián)合分布是多元正態(tài)分布。對于更高階的隨機游走過程，可以得出相似的結果。參數γ二階隨機游走先驗分布定義為：

γj=2γj-1-γj-2+uj,uj～N(0,τ2) j=3,4,…,K

(12)

無信息初始參數的先驗分布為p(γ1)∝const和p(γ2)∝const。根據二階隨機游走的假設，條件分布為：

γj|γj-1,γj-2,…,γ1～N(2γj-1-γj-2,τ2)

(13)

應用貝葉斯方法估計模型參數還需要給出其他未知參數的先驗分布，比如方差參數可以選用逆伽瑪分布，模型中的一般回歸參數可以選用均值為零的正態(tài)分布等。假設模型中各參數的先驗分布相互獨立，可以得出模型各參數的聯(lián)合后驗分布，通過該后驗分布可以得出模型參數的估計值及估計區(qū)間。本文應用蒙特卡羅方法中的Gibbs抽樣，通過模型參數的全條件分布進行重復抽樣，得到所估參數的后驗分布樣本，進而可以得到參數的均值、方差和置信區(qū)間等統(tǒng)計量。

(二)貝葉斯非線性混合效應模型的評價

應用Gibbs抽樣可以得到模型參數后驗分布的一系列樣本，這些樣本可能存在自相關，因而需要檢驗樣本是否收斂。對于MCMC模擬樣本的收斂性檢驗方法有：模擬樣本軌跡圖、自相關檢驗和Gelman-Rubin統(tǒng)計量等[17]。為了得到Gelman-Rubin統(tǒng)計量的估計值，需要同時給出不同初始值生成幾個樣本鏈，比較分析樣本之間和樣本內部的可變性，如果該統(tǒng)計量的結果趨向于1，說明模型參數估計值達到了平穩(wěn)的收斂狀態(tài)。

貝葉斯方法可以給出各參數的后驗分布，從而求得各個參數的估計值及其信度區(qū)間。應用參數估計值的信度區(qū)間就可以對參數進行顯著性檢驗，如果該區(qū)間內含有零值，說明該參數值是不顯著的，具體可參見Kruschke的研究[18]113-117。

一般而言，參數較多的模型對觀察數據的擬合較好，但其預測能力較差。在比較和選擇不同的模型時，通常使用AIC或BIC統(tǒng)計量，它們都含有一個懲罰項來抵消由于參數增加而增加的模型擬合效果。對于固定效應模型，模型的參數個數比較容易確定；對于混合效應模型，模型的有效參數個數并不明確。如果隨機效應的方差較大，該隨機效應參數就比較顯著；反之，如果隨機效應的方差較小，隨機效應參數在模型中的作用就比較小。由于混合效應模型含有隨機效應參數，所以本文應用貝葉斯模型常用的檢驗標準DIC(Deviance information criterion)統(tǒng)計量來對不同模型進行比較和評價，它能夠估計貝葉斯混合效應模型的“有效”參數個數[19-20]。DIC統(tǒng)計量表達形式如下所示：

(14)

四、案例分析

歐洲交通安全理事會每年統(tǒng)計并發(fā)布歐洲31個國家的交通事故傷亡人數數據*參見Ranking EU Progress on Road Safety-8th Road Safety Performance Index Report, http://etsc.eu/ 2014(6)。。本文選用死亡人數數據進行建模分析。由于不同國家的人口數量不同，各國交通死亡人數差別較大，為了便于比較，本文選用百萬人口交通事故死亡率(以下簡稱死亡人數)進行分析。各國死亡人數的觀察期為2001年至2013年，圖1和圖2分別給出了歐洲31個國家13年內交通事故死亡人數分布和各年死亡人數趨勢。

圖1 交通事故死亡人數分布圖

圖2 交通事故死亡人數趨勢圖

從分布圖可以看出，死亡人數是多峰結構，且在觀察期限內具有非線性變化趨勢。下面以死亡人數為響應變量建立回歸模型，可以選擇的解釋變量有人均GDP和觀察年度。在建模過程中發(fā)現(xiàn)人均GDP的作用并不明顯，所以只考慮年度作為解釋變量。

假設死亡人數yij服從泊松分布，解釋變量各觀察年表示為xij，兩個下標i=1,2,…,31和j=1,2,…,13分別表示31個國家和13個觀察年。

基于本文的數據，可以建立下述10個模型進行比較分析。

模型一：假設死亡人數服從泊松分布，對均值參數建立具有對數連接函數的線性混合效應模型：

yij～Poisson(μij)

log(μij)=γi1+γi2xij

模型二：假設死亡人數服從混合泊松分布，其均值參數的回歸方程為線性混合效應模型：

log(μijk)=γik1+γik2xij

模型三：圖2表明死亡人數隨各年度呈現(xiàn)出非線性變化趨勢，所以考慮建立多項式混合效應模型：

yij～Poisson(μij)

模型四：在模型三的基礎上假設死亡人數服從混合泊松分布，同時對不同的均值參數建立回歸方程：

模型五：圖2 顯示出死亡人數隨年度呈現(xiàn)出遞減的變化趨勢，且各條曲線波動明顯，考慮在均值參數的回歸方程中加入截斷多項式，假設截斷多項式參數為隨機效應參數，其他參數為固定效應參數，其中κ1,κ2,…,κm表示節(jié)點。

yij～Poisson(μij)

模型六：在模型五的基礎上假設死亡人數服從混合泊松分布：

模型七：在模型五的基礎上假設回歸方程中所有參數為隨機效應參數：

yij～Poisson(μij)

模型八：在模型七的基礎上假設死亡人數服從混合泊松分布:

模型九：針對死亡人數隨年度呈現(xiàn)出遞減的變化趨勢，考慮在均值參數的回歸方程中加入B樣條，模型中B樣條參數為隨機效應參數，其他參數為固定效應參數：

yij～Poisson(μij)

模型十：在模型九的基礎上，假設死亡人數服從混合泊松分布：

以上10個模型可以分為四類，模型一和模型二是線性混合效應模型；模型三和模型四是多項式混合效應模型；模型五、六、七、八是截斷多項式混合效應模型；模型九、十可以稱為B樣條混合效應模型。

由于這些模型比較復雜，本文使用貝葉斯方法對模型進行估計。在上述模型中，假設混合效應參數(γi1,γi2)的先驗分布為均值為零的二元正態(tài)分布，即(γi1,γi2)～N2(0,Σ)，其中Σ為超參數，其先驗分布為自由度為2的逆沙維特分布(Inv-Wishart)，即Σ～Inv-Wishart2(I)，I表示二維的單位矩陣。模型三和模型四的混合效應參數服從三元正態(tài)分布。

當死亡人數服從混合泊松分布時，模型中含有比例參數p，該參數在(0，1)區(qū)間內取值，將其定義為pk=eλk/(eλ1+eλ2+…+eλC)，其中λ1=0，當k>1時λk～N(0,1)。經過模型比較發(fā)現(xiàn)，當兩個泊松分布進行混合時，模型的擬合效果較好，所以在本文中C=2。

在上述模型中分別采用了1階截斷多項式和3階B樣條，并對13個觀察年度進行了中心化處理。選取的5個節(jié)點分別是-4、-2、0、2、4，所以在截斷多項式模型中m=5，即含有5個截斷項，在B樣條模型中K=9。

根據表1給出部分參數的估計值和模型擬合統(tǒng)計量，可以對上述模型進行比較分析。在模型三和模型四中，參數γ服從三元正態(tài)分布，表1只給出了部分結果，在其他模型中，該參數服從二元正態(tài)分布。每個混合效應參數γ有31個估計值，表1中只給出它的協(xié)方差矩陣Σ。在表1中給出參數估計值的標準誤，進而得出估計值的信度區(qū)間，通過信度區(qū)間確定估計值的統(tǒng)計顯著性，比如大部分參數γ的協(xié)方差在其95%的信度區(qū)間內含有零值，這說明參數γ之間近似相互獨立。參數p給出了兩個泊松分布的混合比例。Deviance(負二倍的對數似然函數)和DIC統(tǒng)計量用于模型比較，這兩個統(tǒng)計量的值越小表示模型越好?？梢钥闯觯旌喜此煞植寄Ｐ偷腄eviance統(tǒng)計量小于同類模型，但是DIC統(tǒng)計量相對較大。這說明雖然混合泊松分布模型的擬合效果較好，但是涉及參數較多，優(yōu)勢并不明顯。從上述十個模型的擬合效果來看，非線性模型的擬合效果優(yōu)于線性模型，樣條多項式優(yōu)于一般多項式，B樣條優(yōu)于截斷樣條?？傮w來看，模型九的Deviance統(tǒng)計量為2 791.28，DIC統(tǒng)計量為2 911，在上述所有模型中是最小的，所以是相對最優(yōu)的模型。另外模型九涉及參數相對較少，便于操作，更具有實際應用價值。

表1 參數估計結果

注：括號中的數字為估計值的標準誤。

五、結 論

由于實際的保險損失縱向數據中，變量之間關系比較復雜，變量間的非線性關系更具有一般性。很多研究為了分析的簡便，假設變量之間存在簡單的線性關系，這并不適合具體的實際應用。本文基于一組實際保險損失數據建立了截斷多項式混合效應模型和B樣條混合效應模型，并與一般的線性混合效應模型進行了比較，結果表明，非線性混合效應模型能夠更好地擬合該組實際數據。在非壽險損失預測和費率厘定實踐中，解釋變量與響應變量之間往往存在非線性變化關系，甚至存在波動性的趨勢變化，在這種情況下應該考慮建立非線性混合效應模型。非線性混合效應模型是對一般線性混合效應模型的擴展，能夠更好地滿足實際應用中對復雜保險損失數據建模的需要，對非壽險費率厘定具有重要價值。

[1] 孟生旺，邱子真. 混合效應模型及其在非壽險費率厘定中的應用[J]. 數理統(tǒng)計與管理, 2016(1).

[2] 王明高，孟生旺. 基于貝葉斯偏態(tài)線性混合模型的費率厘定[J]. 統(tǒng)計與信息論壇，2015(1).

[3] Klein N, Denuit M. Stefan Lang, Thomas Kneib. Nonlife Ratemaking and Risk Management with Bayesian Generalized Additive Models for Location, Scale, and Shape[J].Insurance: Mathematics and Economics, 2014(5).

[4] Hastie T, Tibshirani R. Generalized Additive Models[M]. London: Chapman and Hall, 1990.

[5] Wood S N. Generalized Additive Models: An Introduction with R[M]. London: Chapmann & Hall, 2006.

[6] Ruppert D, Wand M P, Carroll R J. Semiparametric Regression[M]. Combridge: Cambridge University Press,2003.

[7] Wood S N. Fast Stable Direct Fitting and Smoothnesss Election for Generalized Additive Models[J]. Journal of the Royal Statistical Society (B), 2008(4).

[8] Brezger A, Lang S. Generalized Structured Additive Regression Based on Bayesian P-splines[J]. Computational Statistics & Data Analysis, 2006(3).

[9] Jullion A, Lambert P. RobustSpecification of the Roughness Penalty Prior Distribution in Spatially Adaptive Bayesian P-splines Models[J]. Computational Statistics & Data Analysis, 2007(2).

[10]Lang S, Umlauf N, Wechselberger P, et al. Multilevel Structured Additive Regression[J]. Statistics and Computing, 2014(6).

[11]Silva G, Dean C, Niyonsenga T, Vanasse A. Hierarchical Bayesian Spatiotemporal Analysis of Revascularization Odds Using Smoothing Splines[J]. Statistics in Medicine, 2008(2).

[12]Rigby R A, Stasinopoulos D M. Generalized Additive Models for Location, Scale and Shape (with Discussion)[J]. Applied Statistics, 2005(3).

[13]Fahrmeir L, Kneib T, Lang S, Marx B. Regression-Models, Methods and Applications[M]. New York:Springer, 2013.

[14]Goldstein H. Multilevel Statistical Models[M]. 4rd ed. Chichester: John Wiley, 2011.

[15]Beck N, Jackman S. Beyond Linearity by Default: Generalized Additive Models[J]. American Journal of Political Science, 1998(5).

[16]Wu H, Zhang J T. Nonparametric Regression Methods for Longitudinal Data Analysis[M]. Chichester: John Wiley, 2006.

[17]Gelman A,Rubin D. Inference from Iterative Simulation Using Multiple Sequences[J]. Statistical Science, 1992(7).

[18]Kruschke J K. Doing Bayesian Data Analysis[M].London: Academic Press, 2011.

[19]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 Society(B), 2002(2).

[20]Spiegelhalter D J, Best N G, Carlin B P, et al. The Deviance Information Criterion: 12 Years on [J]. Journal of the Royal Statistical Society(B), 2014(1).

(責任編輯：崔國平)

Bayesian Nonlinear Mixed Effects Models and Their Applications

WANG Ming-gao1,2，MENG Sheng-wang2

(1. School of Statistics, Shandong Institute of Business and Technology, Yantai 264005, China;2. Center for Applied Statistics, Renmin University of China, Beijing 100872, China)

Linear mixed effects models have some limitations to model longitudinal data with nonlinear relationship. The paper extends the linear mixed effects models, and based on the nonlinear relationship of variables and the distribution of dependent variable, different nonlinear mixed-effects regression models are established. Using a set of insurance loss data, the paper establishes polynomial mixed effects models, truncated polynomial mixed effects models and B-spline mixed effects models. The result shows that the nonlinear moixed effects models can significantly improve the prediction of the insurance loss. Models established in this paper extends the application of mixed-effects model, and has important value for non-life insurance ratemaking.

mixed effects models; nonlinear models; insured loss; Bayesian

2016-03-25

國家自然科學基金項目《考慮風險相依的非壽險精算模型研究》(71171193)；山東省社會科學基金項目《基于貝葉斯分層MCMC對山東保險市場區(qū)域差異及協(xié)同發(fā)展對策性研究》(15CTJJ01)；山東工商學院博士科研基金項目《貝葉斯分層模型的應用研究》(BS201508)

王明高，男，山東日照人，經濟學博士，講師，研究方向: 風險管理與精算； 孟生旺，男，甘肅秦安人，經濟學博士，教授，博士生導師，研究方向: 風險管理與精算，應用統(tǒng)計。

O212∶F840.3

A

1007-3116(2016)12-0010-07