半參數雙重Tweedie復合泊松回歸模型的貝葉斯分析

段星德,伍震寰,張鐘妮,張文專

(貴州財經大學數學與統計學院,貴州 貴陽 550025)

1.引言

在對健康保險行業進行研究時,人們常常分析衛生保健利用費用數據.而部分被保人在保險期間沒有到醫院進行醫學治療,因此這部分個體沒有產生衛生保健利用費用.上述的衛生保健費用數據就是典型的半連續數據,即由零和正的連續數據所構成.近年來,對衛生保健利用費用數據進行統計建模,已取得了大量的研究成果.首先,Mihaylova等[1]綜述了分析衛生保健資源及費用數據的各種統計方法,這些數據具有偏態、零過多、多峰、重右尾等特點;Smith等[2]利用邊際兩部分模型分析半連續衛生保健服務數據;Neelon等[3-4]綜述了在衛生保健服務領域中零調整計數數據和半連續數據的建模方法及其應用;Merlo等[5]利用兩部分分位數回歸模型來分析半連續衛生保健費用縱向數據.上述文獻中使用的兩部分模型分別對零數據和連續數據進行建模,這樣的分割處理給半連續數據整體屬性的解釋帶來困難.其次,Kurz[6]利用Tweedie回歸模型對半連續衛生保健費用數據進行建模,并與Tobit模型、泊松回歸模型及兩部分模型進行比較分析.

眾所周知,Tweedie復合泊松分布是分析半連續數據的一個重要工具并且具有可解釋半連續數據整體屬性的優勢,因此對Tweedie復合泊松回歸模型的研究引起眾多統計工作者的青睞.一方面,Smyth和J?rgensen[7]以及Andersen和Bonat[8]分別研究了雙重Tweedie復合泊松回歸模型（即對Tweedie復合泊松分布的均值和散度參數聯合建模）統計推斷問題并用這類模型分析半連續保險數據;Halder等[9]在雙重Tweedie復合泊松回歸模型引入空間效應并用它分析半連續保險費率制定數據.另一方面,在貝葉斯框架下,利用Markov Chain Monte Carlo(簡稱MCMC)技術對各類Tweedie復合泊松回歸模型進行統計推斷.比如: Peters等[10]利用Dunn和Smyth[11]給出的的數值方法去逼近Tweedie復合泊松分布的密度函數,并給出這類模型的貝葉斯分析;ZHANG[12],Swallow等[13]以及YE等[14]研究了Tweedie復合泊松隨機效應模型的貝葉斯估計問題;段星德等[15]研究了Tweedie復合泊松回歸模型的貝葉斯數據刪除影響問題.在本文中,基于上述研究工作提出一類半參數雙重Tweedie復合泊松回歸模型,進一步對這類模型進行貝葉斯估計,最后利用這類模型分析衛生保健費用數據以及影響因素.

2.統計模型

本節將首先介紹Tweedie復合泊松分布以及逼近它的密度函數的數值方法,其次引入它們所對應的雙重廣義線性模型: 帶有異質結構的Tweedie復合泊松回歸模型.

Ⅰ Tweedie復合泊松分布

指數分布族是一類常見的分布族,在某些條件下Tweedie復合泊松分布是它的特殊情形.指數分布族的概率函數具有如下的一般形式:

其中,a(·)和k(·)的形式是已知的;θ常被稱作自然參數,?常被稱作離散參數且?>0.另外,指數分布族的均值和方差分別為:μ=E(Y)=k′(θ),var(Y)=?k′′(θ),其中k′(θ)和k′′(θ)表示k(θ)關于未知參數θ的一階導數和二階導數;特別地,函數k′′(θ)稱為方差函數.進一步,如果方差和均值有如下關系var(Y)=?μp,其中參數p是取值范圍為(1,2)的冪指標參數,則有k′′(θ)=μp,θ=μ1-p/(1-p)和k(θ)=μ2-p/(2-p).因此,(2.1)式可以表示為[16]:

如果一個隨機變量Y的概率密度函數具有(2.2)的形式且1

其中V(y)=ypI(y>0)+(y+v0)pI(y=0),這里v0是一個給定的較小的正數[9,17].

Ⅱ 半參數雙重Tweedie復合泊松回歸模型

在本文中考慮以下半參數雙重Tweedie復合泊松回歸模型,即對Tweedie復合泊松分布的均值參數和散度參數進行聯合建模:

其中,Y=(Y1,Y2,···,Ym)T是m維響應變量且相互獨立,xi=(xi1,xi2,···,xik)T表示均值模型中的協變量,zi=(zi1,zi2,···,ziq)T表示散度模型中的協變量.β=(β1,β2,···,βk)T,ζ=(ζ1,ζ2,···,ζq)T分別是k×1和q×1維未知待估參數向量,且k

3.貝葉斯分析

Ⅰ 先驗分布與后驗分布

其中μβ,Σβ,μζ,Σζ,aτ,bτ,aδ,bδ為已知的超參數,IG表示逆Gamma分布,Γ(a,b) 表示服從參數為a和b的Gamma分布.此外,本文在抽樣過程中使用的條件分布、Gibbs抽樣、MH算法如附錄所示.

Ⅱ 貝葉斯估計

為了得到平穩的隨機序列,我們舍棄序列的前D個值,并保留來自聯合后驗分布p(θ|Y,X,Z,T)的隨機樣本{θ(n) :n=D+1,···,N},則有

4.數值例子

Ⅰ 模擬研究

在模擬研究中,假設從如下模型結構中產生半連續響應數據{yi,i=1,2,···,m},即:

其中樣本量m=200,協變量xi ～N(0,1),zi ～N(0,1),ti～U(0,1).令參數和非參數真實值為β=(β0,β1)T=(0.5,1.5)T,ζ=(ζ0,ζ1)T=(-1,0.4)T,p=1.6,g(ti)=1.5×cos(2πti).

在貝葉斯框架下進行的模擬研究中,我們通常研究下面三種不同先驗信息對貝葉斯估計的影響,即:

類型Ⅰ(良好的先驗信息) 設超參數μβ的取值為β的真值,即μβ=β=(0.5,1.5)T,協方差陣Σβ=0.25I2,I2表示2階單位陣;μζ=ζ=(-1,0.4)T,Σζ=0.25I2;aτ=1,bτ=0.005,aδ=0.5,bδ=0.5.

類型Ⅱ(不準確的先驗信息) 設μβ=1.5×β=1.5×(0.5,1.5)T,協方差陣Σβ=0.75I2,I2表示2階單位陣;μζ=1.5×ζ=1.5×(-1,0.4)T,Σζ=0.75I2;其它超參數的取值和類型Ⅰ一致.

類型Ⅲ(無先驗信息) 設μβ=(0,0)T,協方差陣Σβ=100I2,I2表示2階單位陣;μζ=(0,0)T,Σζ=100I2;其它超參數的取值和類型Ⅰ一致.

基于上述三種類型的先驗信息,我們分別做100次實驗的模擬研究,且每次實驗都迭代10000次,為了避免最初產生的非平穩樣本序列對后驗推斷的影響,我們舍棄前面產生的5000次迭代值,利用后面的5000 次迭代值來進行貝葉斯估計.另外,在實施MH算法時,我們選擇方差調節參數=4,=0.7,=2,=0.6使得在抽樣過程中所有參數的平均接收率在區間[0.26,0.35]上.表1給出了所有參數的Bayes估計、標準差和RMS(表示參數的Bayes估計與真值的偏差的平方的平均值的算術平方根).從表1中發現: 在三類不同先驗信息下,所有參數的貝葉斯估計與真值的偏差都很小,說明我們模擬研究中所得到的貝葉斯估計都具有較高的精度且對先驗信息不敏感;另外,參數的標準差和RMS 值也比較接近.我們在圖1中列出類型Ⅰ、類型Ⅱ和類型Ⅲ先驗下非參數g(t)的估計值與真實值的擬合圖形,從圖1中發現,非參數部分關于真實函數的擬合是比較好的,說明我們所使用的貝葉斯P-樣條方法是有效的.

圖1 類型I (左圖),類型II (中圖),類型III(右圖)時非參數函數g(t)擬合圖

表1 隨機模擬研究中未知參數的Bayes估計

Ⅱ 實證分析

研究的數據來源于蘭德健康保險實驗(RAND HIE),該實驗是對美國醫療成本、衛生保健利用率及相關結果的一個綜合研究.[6]為了設計可靠的實驗和得到精準的數據,該項研究跟蹤了隨機分配到不同計劃的人群并記錄了他們的醫療費用及個人信息.這里,我們選擇了第五年觀察期的1713個個體作為樣本,并用ID標識不同個體.數據集中,衛生保健費用包括如下5種: 門診費用(‘outpdol’)、藥物費用(‘drugdol’)、供應費用(‘suppdol’)、心理治療費用(‘mentdol’)和住院費用(‘indol’),我們把每個個體的5種衛生保健費用之和作為響應變量,并記為yi.另外,把個體信息: 性別(‘female’:1=女性，0=男性)、種族(‘black’:1=戶主是黑人，0=戶主不是黑人)、家庭收入對數(‘linc’)、身體缺陷數(‘physlm’)、慢性病數(‘disea’)、家庭規模對數(‘lfam’)、戶主受教育年限(‘educdec’)和表示自評健康狀況良好的虛擬變量(‘hlthg’)、保險特定變量包括對數共同保險率加1(‘logc’)、個人免賠額計劃(‘idp’)的虛擬值、參與激勵支付(‘lpi’)的對數和最大支出函數(‘fmde’)作為協變量,并把每個個體對應的12個協變量表示為xi1,xi2,xi3,xi4,xi5,xi6,xi7,xi8,xi9,xi10,xi11,xi12,把年齡(‘age’)作為非參數函數里的時間變量進行考慮并記為ti.在建模過程中,假定zij=xij,j=1,2,···,12,則用如下的模型擬合上述數據集:

同樣,利用上述提出的MH算法和Gibbs抽樣從聯合后驗分布中產生隨機序列,在10000次迭代值中剔除前5000次,把后5000次迭代值作為后驗樣本值進行貝葉斯估計,計算結果見表2.

表2 實例分析中參數的估計值和標準差

由表2中的均值模型可以看出,對衛生保健費用支出有顯著影響的協變量有家庭收入對數、身體缺陷數、慢性病數、保險特定變量、最大支出函數,而且前四個協變量對應的回歸系數的符號都為正,這與定性分析的結果是一致的.另一方面,從散度模型中可以看出,性別、種族、家庭規模對數、個人免賠額計劃的虛擬值、參與激勵支付的對數這些協變量都是顯著的;因此我們對散度參數進行建模也是合理的.非參數光滑函數g(age)的貝葉斯估計結果見圖2,由函數g(age)的P樣條估計可見其具有明顯的非線性效應.

圖2 實證分析中非參數函數g(age)的貝葉斯估計

5.結論

本文利用Gibbs抽樣技術和Metropolis-Hastings(MH)算法的混合算法去研究半參數雙重Tweedie復合泊松回歸模型的聯合貝葉斯估計問題.最后利用提出的方法分析蘭德健康保險實驗中的半連續衛生保健費用數據.

6.附錄

6.1 條件分布和Gibbs抽樣

為了利用Gibbs抽樣技術從后驗分布進行抽樣,我們首先基于(3.1)-(3.2)式推導出如下步驟的滿條件分布:

步1 給定Y,X,Z,T,p,ζ,ξ下β的滿條件分布為

步2 給定Y,X,Z,T,β,p,ξ下ζ的條件分布為

步3 給定Y,X,Z,T,ζ,β,ξ下p的條件分布為

步4 給定Y,X,Z,T,β,ζ,,δ下ξ的條件分布為

步5 給定ξ,δ下的條件分布為

6.2 MH算法

為了利用MH算法對(6.1)-(6.4)給出的滿條件分布進行抽樣,步驟如下:

第二,從均勻分布U(0,1)中產生隨機數u,若u ≤ a(β(n),β?),令β(n+1)=β?,否則令β(n+1)=β(n),其中接受概率為:

同樣,參數ζ,p,ξ可用上述方法產生隨機樣本.