半参数双重Tweedie复合泊松回归模型的贝叶斯分析.pdf

应用数学MATHEMATICA APPLICATA2024,37(1):272-279半参数双重Tweedie复合泊松回归模型的贝叶斯分析段星德,伍震寰,张钟妮,张文专(贵州财经大学数学与统计学院,贵州 贵阳 550025)摘要：为了分析健康保险行业中出现的半连续卫生保健费用数据,本文提出一类半参数双重Tweedie复合泊松回归模型.在分析中,首先采用修正鞍点逼近的数值方法去近似Tweedie复合泊松分布的密度函数;其次,利用Gibbs抽样技术和Metropolis-Hastings(MH)算法的混合算法获得了模型参数的联合贝叶斯估计;最后,给出了几个模拟研究以及把这些方法用来分析兰德健康保险实验中的卫生保健费用数据.关键词：卫生保健利用;复合泊松分布;半连续数据;MH算法;Gibbs抽样;贝叶斯P-样条中图分类号：O212.8AMS(2010)主题分类：62F15;62J12文献标识码：A文章编号：1001-9847(2024)01-0272-081.引言在对健康保险行业进行研究时,人们常常分析卫生保健利用费用数据.而部分被保人在保险期间没有到医院进行医学治疗,因此这部分个体没有产生卫生保健利用费用.上述的卫生保健费用数据就是典型的半连续数据,即由零和正的连续数据所构成.近年来,对卫生保健利用费用数据进行统计建模,已取得了大量的研究成果.首先,Mihaylova等1综述了分析卫生保健资源及费用数据的各种统计方法,这些数据具有偏态、零过多、多峰、重右尾等特点;Smith等2利用边际两部分模型分析半连续卫生保健服务数据;Neelon等3-4综述了在卫生保健服务领域中零调整计数数据和半连续数据的建模方法及其应用;Merlo等5利用两部分分位数回归模型来分析半连续卫生保健费用纵向数据.上述文献中使用的两部分模型分别对零数据和连续数据进行建模,这样的分割处理给半连续数据整体属性的解释带来困难.其次,Kurz6利用Tweedie回归模型对半连续卫生保健费用数据进行建模,并与Tobit模型、泊松回归模型及两部分模型进行比较分析.众所周知,Tweedie复合泊松分布是分析半连续数据的一个重要工具并且具有可解释半连续数据整体属性的优势,因此对Tweedie复合泊松回归模型的研究引起众多统计工作者的青睐.一方面,Smyth和Jrgensen7以及Andersen和Bonat8分别研究了双重Tweedie复合泊松回归模型（即对Tweedie复合泊松分布的均值和散度参数联合建模）统计推断问题并用这类模型分析半连续保险数据;Halder等9在双重Tweedie复合泊松回归模型引入空间效应并收稿日期：2023-02-26基金项目：国家自然科学基金项目(12161014);全国统计科学研究项目(2021LY011);贵州省省级科技计划项目资助(黔科合基础20201Y009);贵州省教育厅自然科学基金资助项目（KY2021134）通讯作者：段星德,男,彝族,云南人,教授,研究方向:贝叶斯统计计算.第 1 期段星德等：半参数双重Tweedie复合泊松回归模型的贝叶斯分析273用它分析半连续保险费率制定数据.另一方面,在贝叶斯框架下,利用Markov Chain MonteCarlo(简称MCMC)技术对各类Tweedie复合泊松回归模型进行统计推断.比如:Peters等10利用Dunn和Smyth11给出的的数值方法去逼近Tweedie复合泊松分布的密度函数,并给出这类模型的贝叶斯分析;ZHANG12,Swallow等13以及YE等14研究了Tweedie复合泊松随机效应模型的贝叶斯估计问题;段星德等15研究了Tweedie复合泊松回归模型的贝叶斯数据删除影响问题.在本文中,基于上述研究工作提出一类半参数双重Tweedie复合泊松回归模型,进一步对这类模型进行贝叶斯估计,最后利用这类模型分析卫生保健费用数据以及影响因素.2.统计模型本节将首先介绍Tweedie复合泊松分布以及逼近它的密度函数的数值方法,其次引入它们所对应的双重广义线性模型:带有异质结构的Tweedie复合泊松回归模型.Tweedie复合泊松分布指数分布族是一类常见的分布族,在某些条件下Tweedie复合泊松分布是它的特殊情形.指数分布族的概率函数具有如下的一般形式:p(y;,)=a(y,)expy k(),(2.1)其中,a()和k()的形式是已知的;常被称作自然参数,常被称作离散参数且 0.另外,指数分布族的均值和方差分别为:=E(Y)=k(),var(Y)=k(),其中k()和k()表示k()关于未知参数的一阶导数和二阶导数;特别地,函数k()称为方差函数.进一步,如果方差和均值有如下关系var(Y)=p,其中参数p是取值范围为(1,2)的幂指标参数,则有k()=p,=1p/(1 p)和k()=2p/(2 p).因此,(2.1)式可以表示为16:p(y;,p,)=a(y,p)exp1(y1p1 p2p2 p),(2.2)如果一个随机变量Y 的概率密度函数具有(2.2)的形式且1 p 0)+(y+v0)pI(y=0),这里v0是一个给定的较小的正数9,17.半参数双重Tweedie复合泊松回归模型在本文中考虑以下半参数双重Tweedie复合泊松回归模型,即对Tweedie复合泊松分布的均值参数和散度参数进行联合建模:Yi Twp(i,i)log(i)=xiT+g(ti)log(i)=ziT(2.5)其中,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 m,q m,g()是未知的光滑可导的连续函数.另外,E(Yi)=i,var(Yi)=pi,p (1,2).为了估计参数 和,首先需要估计光滑函数g().根据Lang等18的结论,我们用如下的贝叶斯P-样条去逼近光滑函数:g(t)=Jj=1Bj(t)j=BT(t),(2.6)274应用数学2024其中,Bj(t)表示第j个B样条基函数,J表示B样条的条数(等于变量t的节点数与B样条的阶数之和);且有B(t)=(B1(t),BJ(t)T,=(1,J)T表示待估的样条系数.为了估计样条系数,我们对j定义如下的一阶随机游动模型j=j1+j,其中j独立且服从N(0,2/j),j=2,J.这里,我们称方差参数2/j为自适应光滑参数来控制非参函数的光滑程度.令=(2,J),可把样条系数的先验分布表示为|2,exp(122TQ),其中,Q表示惩罚矩阵.3.贝叶斯分析先验分布与后验分布令Y=(Y1,Y2,Ym)T,X=(x1,x2,xm)T,Z=(z1,z2,zm)T,T=(t1,t2,tm)T,=(,p,2,j:j=2,J).在贝叶斯框架下,参数的后验分布可表示为:p(|Y,X,Z,T)mi=1p(yi|xi,zi,ti,)p()mi=1(2iV(yi)1/2expd(yi,i)2ip(),(3.1)其中p()为参数的先验分布且给定为:p()N(,),p()N(,),log(p12p)N(0,100),|2 exp(122TQ),p(2)IG(a,b),p(j)(a,b),(3.2)其中,a,b,a,b为已知的超参数,IG表示逆Gamma分布,(a,b)表示服从参数为a和b的Gamma分布.此外,本文在抽样过程中使用的条件分布、Gibbs抽样、MH算法如附录所示.贝叶斯估计为了得到平稳的随机序列,我们舍弃序列的前D个值,并保留来自联合后验分布p(|Y,X,Z,T)的随机样本(n):n=D+1,N,则有b=1N DNn=D+1(n),cvar(|Y,X,Z,T)=1N D 1Nn=D+1(n)b)(n)b)T.(3.3)上式中b可作为参数的联合Bayes估计,而用样本协方差矩阵cvar(|Y,X,Z,T)的对角线元素来估计参数的标准差.4.数值例子模拟研究在模拟研究中,假设从如下模型结构中产生半连续响应数据yi,i=1,2,m,即:Yi Twp(i,i),log(i)=0+1xi+g(ti),log(i)=0+1zi,其中样本量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(2ti).在贝叶斯框架下进行的模拟研究中,我们通常研究下面三种不同先验信息对贝叶斯估计的影响,即:类型(良好的先验信息)设超参数的取值为的真值,即=(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 期段星德等：半参数双重Tweedie复合泊松回归模型的贝叶斯分析275类型(不准确的先验信息)设=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算法时,我们选择方差调节参数2=4,2p=0.7,2=2,2=0.6使得在抽样过程中所有参数的平均接收率在区间0.26,0.35上.表1给出了所有参数的Bayes估计、标准差和RMS(表示参数的Bayes估计与真值的偏差的平方的平均值的算术平方根).从表1中发现:在三类不同先验信息下,所有参数的贝叶斯估计与真值的偏差都很小,说明我们模拟研究中所得到的贝叶斯估计都具有较高的精度且对先验信息不敏感;另外,参数的标准差和RMS 值也比较接近.我们在图1中列出类型、类型和类型先验下非参数g(t)的估计值与真实值的拟合图形,从图1中发现,非参数部分关于真实函数的拟合是比较好的,说明我们所使用的贝叶斯P-样条方法是有效的.表1随机模拟研究中未知参数的Bayes估计Type IType IIType IIIPar.BiasSDRMSBiasSDRMSBiasSDRMS0-0.00250.08750.0871-0.00290.08810.0877-0.00280.08840.08801-0.00160.03750.0373-0.00040.03760.0374-0.00170.03760.037500.07380.10650.12910.06840.10990.12900.07730.11040.134410.01840.11530.11620.02200.11820.11970.01920.12110.1220p-0.05210.04590.0692-0.05180.04620.0692-0.05210.04590.0693RMS0.43890.44300.451200.10.20.30.40.50.60.70.80.91tij-1.5-1-0.500.511.5g(tij)true valueestimate00.10.20.30.40.50.60.70.80.91ti-1.5-1-0.500.511.5g(ti)true valueestimate00.10.20.30.40.50.60.70.80.91ti-1.5-1-0.500.511.5g(ti)true valueestimate图1类型I(左图),类型II(中图),类型III(右图)时非参数函数g(t)拟合图实证分析研究的数据来源于兰德健康保险实验(RAND HIE),该实验是对美国医疗成本、卫生保健利用率及相关结果的一个综合研究.6为了设计可靠的实验和得到精准的数据,该项研究跟踪了随机分配到不同计划的人群并记录了他们的医疗费用及个人信息.这里,我们选择了第五年观察期的1713个个体作为样本,并用ID标识不同个体.数据集中,卫生保健费用包括如下5种:门诊费用(outpdol)、药物费用(drugdol)、供应费用(suppdol)、心理治疗费用(mentdol)和住院费用(indol),我们把每个个体的5种276应用数学2024卫生保健费用之和作为响应变量,并记为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,则用如下的模型拟合上述数据集:yi Twp(i,i),log(i)=0+1xi1+2xi2+3xi3+4xi4+5xi5+6xi6+7xi7+8xi8+9xi9+10 xi10+11xi11+12xi12+g(ti),log(i)=0+1zi1+2zi2+3zi3+4zi4+5zi5+6zi6+7zi7+8zi8+9zi9+10zi10+11zi11+12zi12.同样,利用上述提出的MH算法和Gibbs抽样从联合后验分布中产生随机序列,在10000次迭代值中剔除前5000次,把后5000次迭代值作为后验样本值进行贝叶斯估计,计算结果见表2.由表2中的均值模型可以看出,对卫生保健费用支出有显著影响的协变量有家庭收入对数、身体缺陷数、慢性病数、保险特定变量、最大支出函数,而且前四个协变量对应的回归系数的符号都为正,这与定性分析的结果是一致的.另一方面,从散度模型中可以看出,性别、种族、家庭规模对数、个人免赔额计划的虚拟值、参与激励支付的对数这些协变量都是显著的;因此我们对散度参数进行建模也是合理的.非参数光滑函数g(age)的贝叶斯估计结果见图2,由函数g(age)的P样条估计可见其具有明显的非线性效应.表2实例分析中参数的估计值和标准差均值模型散度模型参数估计值标准差参数估计值标准差04.34340.354302.12320.296210.10800.08451-0.28600.06802-0.12260.157720.62520.092530.09920.037530.04760.031440.57410.14554-0.04190.112350.02400.00685-0.00830.00596-0.14130.10616-0.17110.06777-0.00270.02027-0.00860.01438-0.05260.090580.03530.074890.18590.068990.06300.047110-0.22570.1316100.38010.094611-0.03610.021711-0.03010.014812-0.13010.0400120.01310.0280第 1 期段星德等：半参数双重Tweedie复合泊松回归模型的贝叶斯分析277010203040506070age-1-0.8-0.6-0.4-0.200.20.40.60.81nonparameter estimate图2实证分析中非参数函数g(age)的贝叶斯估计5.结论本文利用Gibbs抽样技术和Metropolis-Hastings(MH)算法的混合算法去研究半参数双重Tweedie复合泊松回归模型的联合贝叶斯估计问题.最后利用提出的方法分析兰德健康保险实验中的半连续卫生保健费用数据.6.附录6.1 条件分布和Gibbs抽样为了利用Gibbs抽样技术从后验分布进行抽样,我们首先基于(3.1)-(3.2)式推导出如下步骤的满条件分布:步1给定Y,X,Z,T,p,下的满条件分布为p(|Y,X,Z,T,p,)expmi=11i(2pi2 pyi1pi1 p)12()T1()(6.1)步2给定Y,X,Z,T,p,下的条件分布为p(|Y,X,Z,T,p,)expmi=11/2imi=11i(2pi2 pyi1pi1 p)12()T1()(6.2)步3给定Y,X,Z,T,下p的条件分布为p(p|Y,X,Z,T,)expmi=1log(V(yi)mi=11i(2pi2 pyi1pi1 p)log(p 1)(2 p)12 100(log(p 12 p)2(6.3)步4给定Y,X,Z,T,2,下的条件分布为p(|Y,X,Z,T,2,)expmi=11i(2pi2 pyi1pi1 p)122TQ.(6.4)步5给定,下2的条件分布为p(2|,)(2)rk(Q)2a1exp122TQ b2,(6.5)278应用数学2024于是有p(2|,)IG(rk(Q)2+a,12TQ+b),其中IG为逆伽马分布.步6给定,2下l的条件分布为p(j|,2)(j)a12expj(122j2+b),(6.6)于是有p(j|,2)(a+12,122j2+b),其中为伽马分布,j=j j1.由(6.1)-(6.6)可以看出,p(|Y,X,Z,T,p,),p(|Y,X,Z,T,p,),p(p|Y,X,Z,T,),p(|Y,X,Z,T,2,)并不是标准分布,因此,我们采用下述介绍的MH算法来对这些非标准后验分布进行随机抽样.6.2 MH算法为了利用MH算法对(6.1)-(6.4)给出的满条件分布进行抽样,步骤如下:第一,分别从提案分布:N(n),2),N(n),2),N(log(p(n)12p(n),2p),N(n);2)中产生随机备选值,p,其中2,2,2p,和2称为方差调节参数使得MH算法的平均接收率在区间0.26,0.35上;另外,上述提案多元正态分布的协方差阵可计算得:=(mi=12pi2pixTixi+mi=1p1iyi1pixTixi+1)1,=(12pmi=11ii2pzTizi+1p1mi=11iyi1pizTizi+1)1,=(2 p)mi=11i2piB(ti)TB(ti)+(p 1)mi=11iyi1piB(ti)TB(ti)+12Q)1.第二,从均匀分布U(0,1)中产生随机数u,若u a(n),),令(n+1)=,否则令(n+1)=(n),其中接受概率为:a(n),)=min1,p(|Y,X,Z,T,p,)p(n)|Y,X,Z,T,p,),同样,参数,p,可用上述方法产生随机样本.参考文献：1 MIHAYLOVA B,BRIGGS A,OHAGAN A,THOMPSON S G.Review of statistical methods foranalysing healthcare resources and costsJ.Health Economics,2011,20(8)：897-916.2 SMITH V,PREISSER J,NEELON B,MACIEJEWSKI M.A marginalized two-part model for semi-continuous dataJ.Statistics in Medicine,2014,33(28)：4891-4903.3 NEELON B,OMALLEY A,SMITH V.Modeling zero-modified count and semicontinuous datain health services research Part 1:background and overviewJ.Statistics in Medicine,2016,35(27)：5070-5093.4 NEELON B,OMALLEY A,SMITH V.Modeling zero-modified count and semicontinuous data inhealth services research Part 2:case studiesJ.Statistics in Medicine,2016,35(27)：5094-5112.5 MERLO L,MARUOTTI A,PETRELLA L.Two-part quantile regression models for semi-continuouslongitudinal data:A finite mixture approachJ.Statistical Modelling,2022,22(6):485-508.6 KURZ C F.Tweedie distributions for fitting semicontinuous health care utilization cost dataJ.BMCMedical Research Methodology,2017,17(1):171.7 SMYTH G K,JRGENSEN B.Fitting Tweedies compound Poisson model to insurance claims data:dispersion modellingJ.ASTIN Bulletin:The Journal of the IAA,2002,32(1):143-157.8 ANDERSEN D A,BONAT W H.Double generalized linear compound poisson models to insuranceclaims dataJ.Electronic Journal of Applied Statistical Analysis,2017,10(2):384-407.第 1 期段星德等：半参数双重Tweedie复合泊松回归模型的贝叶斯分析2799 HALDER A,MOHAMMED S,CHEN K,DEY D K.Spatial Tweedie exponential dispersion models:an application to insurance rate-makingJ.Scandinavian Actuarial Journal,2021,2021(10):1017-1036.10 PETERS G W,SHEVCHENKO P V,WUTHRICH M V.Model uncertainty in claims reservingwithin Tweedies compound Poisson modelsJ.ASTIN Bulletin:The Journal of the IAA,2009,39(1):1-33.11 DUNN P K,SMYTH G K.Series evaluation of Tweedie exponential dispersion model densitiesJ.Statistics and Computing,2005,15(4):267-280.12 ZHANG Y.Likelihood-based and Bayesian methods for Tweedie compound Poisson linear mixedmodelsJ.Statistics and Computing,2013,23(6):743-757.13 SWALLOW B,BUCKLAND S T,KING R,TOMS MP.Bayesian hierarchical modelling of continuousnon-negative longitudinal data with a spike at zero:An application to a study of birds visiting gardensin winterJ.Biometrical Journal,2016,58(2):357-371.14 YE T,LACHOS V H,WANG X,DEY D K.Comparisons of zero-augmented continuous regressionmodels from a Bayesian perspectiveJ.Statistics in Medicine,2021,40(5):1073-1100.15 段星德,张实,罗露璐,张文专.Tweedie复合泊松回归模型的Bayes估计和影响分析J.高校应用数学学报,2020,35(4):393-404.16 JRGENSEN,B.The Theory of Dispersion ModelsM.London:Chapman and Hall,1997.17 DUNN P K,SMYTH G K.Generalized Linear Models With Examples in RM.New York:Springer,2018.18 LANG S,BREZGER A.Bayesian P-splinesJ.Journal of Computational and Graphical Statistics,2004,13(1):183-212.A Semiparametric Bayesian Approach to Double TweedieCompound Poisson Regression ModelDUAN Xingde,WU Zhenhuan,ZHANG Zhongni,ZHANG Wenzhuan(School of Mathematics and Statistics,Guizhou University of Finance and Economics,Guiyang 550025,China)Abstract:We introduce a new class of semiparametric double generalized linear models basedon Tweedie Compound Poisson distribution,accommodating semicontinuous health care cost data inthe health actuarial field.Firstly,a modified saddle-point approximation approach is used to evaluatethe density function of the Tweedies Compound Poisson distribution.In addition,a hybrid algorithmcombining the block Gibbs sampler and the Metropolis-Hastings algorithm is proposed for producing thejoint Bayesian estimates of unknown parameters and nonparametric function.Finally,several simulationstudies and a real example with the RAND Health Insurance Experiment are presented to illustrate theproposed methodologies.Key words:Health care utilization;Compound Poisson distribution;Semicontinuous data;Metropolis-Hastings algorithm;Gibbs sampler;Bayesian P-spline