半?yún)?shù)雙重Tweedie復(fù)合泊松回歸模型的貝葉斯分析

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

(貴州財(cái)經(jīng)大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,貴州 貴陽 550025)

1.引言

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

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

2.統(tǒng)計(jì)模型

本節(jié)將首先介紹Tweedie復(fù)合泊松分布以及逼近它的密度函數(shù)的數(shù)值方法,其次引入它們所對應(yīng)的雙重廣義線性模型: 帶有異質(zhì)結(jié)構(gòu)的Tweedie復(fù)合泊松回歸模型.

Ⅰ Tweedie復(fù)合泊松分布

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

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

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

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

Ⅱ 半?yún)?shù)雙重Tweedie復(fù)合泊松回歸模型

在本文中考慮以下半?yún)?shù)雙重Tweedie復(fù)合泊松回歸模型,即對Tweedie復(fù)合泊松分布的均值參數(shù)和散度參數(shù)進(jìn)行聯(lián)合建模:

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

3.貝葉斯分析

Ⅰ 先驗(yàn)分布與后驗(yàn)分布

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

Ⅱ 貝葉斯估計(jì)

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

4.數(shù)值例子

Ⅰ 模擬研究

在模擬研究中,假設(shè)從如下模型結(jié)構(gòu)中產(chǎn)生半連續(xù)響應(yīng)數(shù)據(jù){yi,i=1,2,···,m},即:

其中樣本量m=200,協(xié)變量xi ～N(0,1),zi ～N(0,1),ti～U(0,1).令參數(shù)和非參數(shù)真實(shí)值為β=(β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).

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

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

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

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

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

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

表1 隨機(jī)模擬研究中未知參數(shù)的Bayes估計(jì)

Ⅱ 實(shí)證分析

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

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

表2 實(shí)例分析中參數(shù)的估計(jì)值和標(biāo)準(zhǔn)差

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

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

5.結(jié)論

本文利用Gibbs抽樣技術(shù)和Metropolis-Hastings(MH)算法的混合算法去研究半?yún)?shù)雙重Tweedie復(fù)合泊松回歸模型的聯(lián)合貝葉斯估計(jì)問題.最后利用提出的方法分析蘭德健康保險(xiǎn)實(shí)驗(yàn)中的半連續(xù)衛(wèi)生保健費(fèi)用數(shù)據(jù).

6.附錄

6.1 條件分布和Gibbs抽樣

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

步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)給出的滿條件分布進(jìn)行抽樣,步驟如下:

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

同樣,參數(shù)ζ,p,ξ可用上述方法產(chǎn)生隨機(jī)樣本.