广义极值回归模型下现状数据的贝叶斯估计

孙 烨， 蒋京京， 王纯杰

(长春工业大学数学与统计学院， 吉林长春130012)

生存分析领域主要研究感兴趣事件的发生时间，比如临床医学中病人的死亡时间或细菌的感染时间、可靠性试验中机器元件的失效时间等[1]。这些事件的发生时间都可归结到生存分析研究领域，由此可见，生存分析应用的领域十分广泛。极值分布[2]是应用研究较为广泛的参数模型，可应用于生存数据分析。作为极值分布的推广，广义极值分布对处理删失数据具有很好分析效果，而且它用统一的表达形式，涵盖了极值I型(Gumbel)分布、极值II型(Frechet)分布和极值III型(Weibull)分布3种类型，其中形状参数符号决定分布类型，因而在寿命试验和可靠性工程中有极其重要应用。同时，在生存分析中，针对不同删失数据类型下的感兴趣变量，研究协变量对感兴趣变量的影响是否显著，就要建立相应生存分析背景下的回归模型。李群等[3]对广义指数分布下的区间删失II型数据进行贝叶斯回归分析，基于广义指数分布的尺度参数引入协变量，建立尺度参数与生存时间的贝叶斯回归模型，关注协变量对于生存时间的影响。

Prescott等[4-5]于1980年和1983年分别基于完整数据和删失数据对广义极值分布的参数进行极大似然估计，并给出三参数广义极值分布信息矩阵的显式表达式；李秀敏等[6]基于完整数据对广义极值分布的参数进行贝叶斯估计，并将其应用到某水文观测站历年的年最高水位观测数据中，得到较好的估计结果；Markose等[7]使用广义极值分布代替对数正态分布研究股票收益率，并开发了一个广义极值分布下的定价模型；鲁帆等[8]基于完整数据下广义极值分布模型和MH抽样算法进行极值洪水的频率分析研究，说明提议分布中的提议方差的选取并不影响计算效果，采用分位数图等多种拟合优度检验方法进行检验，其拟合结果显示贝叶斯方法优于极大似然估计方法和矩法等经典统计方法；樊利利等[9]基于完整数据对海港年最大海平面高度建立广义极值模型，利用R软件进行求解并比较极大似然估计、广义极大似然估计和贝叶斯估计的估计效果，结果表明贝叶斯估计效果相对较好，具有更小的均方根误差；吴云标等[10]在完整数据下对广义极值分布的参数基于MH抽样算法进行贝叶斯估计，给出洪水的后验分布，相较于传统极大似然方法，贝叶斯方法不仅可以得到洪水设计值的点估计，也可得到区间估计，进一步说明了应用贝叶斯方法进行分析结果的可靠性；Naseef等[11]在完整数据下基于广义极值分布和广义帕累托分布研究印度陆架海50年重现期的有效波高的不确定性，并以50年回归值为例，研究热带气旋的影响；Aguirre-Salado等[12]采用惩罚极大似然法拟合广义极值参数的平滑函数，分析2000—2018年墨西哥城市城区紫外辐射的空间极值，结果表明惩罚极大似然法比极大似然估计对模型具有更好的正则化效果；Rypkema等[13]在完整数据下，基于广义极值分布描述气候变化预测，并将广义极值纳入生态人口模型；Sampaio等[14]在完整数据下基于广义极值分布建立贝叶斯层次模型，并将该模型应用到区域洪水频率分析中。上述文献均没涉及广义极值分布的贝叶斯生存回归模型。本文主要研究现状数据下广义极值模型的贝叶斯估计问题，并与极大似然方法进行比较。

1 现状数据下广义极值分布的参数估计

1.1 符号描述和模型介绍

本文主要研究现状数据下广义极值回归模型的建立及贝叶斯估计，首先介绍广义极值分布的一般形式，广义极值分布[2]的分布函数表示为

(1)

特别地，μ=0、σ=1时称为标准广义极值分布，其分布函数为

(2)

当γ=0时，F(t)为极值Ⅰ型(Gumbel)分布

(3)

当γ>0时，取α=1/γ，F(t)为极值Ⅱ型(Frechet)分布

(4)

当γ<0时，取α=-1/γ，F(t)为极值Ⅲ型(Weibull)分布

(5)

由式(1)可得广义极值分布的生存函数为

(6)

广义极值分布的概率密度函数为

(7)

采用广义极值分布对现状数据进行建模，对第i个个体，用Ti表示第i个个体感兴趣的生存时间，Ci表示第i个个体的观测时间点，vi为示性变量，vi=ITi≤Ci，xi表示p×1维协变量向量。因此观测到的现状数据结构为D=Di=Ci,vi,xi,i=1,2,…,n。

在现状数据下，广义极值分布的似然函数可表示为

因此，对数似然函数为

LLγ,σ,βD=lnLγ,σ,βD=

为简化计算，令

用极大似然估计法极大化求解，估计参数θ=(γ,σ,β)T。针对对数似然函数，LL(γ、σ,βD)分别关于参数γ、σ和β求一阶偏导数，并且令其一阶偏导数等于0，得到以下似然方程组

(8)

下面介绍贝叶斯估计方法在现状数据下的研究。

1.2 贝叶斯估计及MCMC算法

本节采用MCMC方法估计模型中的参数，待估参数可表示为θ=(γ,σ,β)T。这里只考虑协变量xi为一维(p=1)的情况，此算法也较易推广到p维。给待估参数一个先验，由贝叶斯定理得到后验分布。选取Gibbs算法和MH算法相结合的方式，从后验分布中抽取随机数，进行后验推断。本文σ和β的先验分布选择参考樊利利等[9]的方法

所以，θ=(γ,σ,β)T的联合后验密度为

π(γ,σ,βD)∝L(γ,σ,βD)π(γa0,b0)π(σc0,d0)π(βμ0,σ0)∝

θ=(γ,σ,β)T的联合对数后验密度为

接下来采用MCMC方法对上述模型中的各个参数进行估计，其中Gibbs采样方法特别适合于目标分布是多元的场合，当参数的条件后验密度函数形式较为复杂时，无法从某个特定分布中抽取，则运用MH算法抽样。所以，为了从条件后验分布中对各个参数进行抽样，需要在Gibbs-sampler中执行Metropolis-Hastings 算法，MH算法具体步骤如下：

① 为待估参数任意选取一个初值θ0=γ0,σ0,β(0)T，并且定义t=0；

② 从提议分布g·θt中抽取θ(new)，提议分布选取正态分布；

③ 从U(0,1)生成随机数u，计算接受概率，如果u≤A，则接受θ(new)，且定义θt+1=θ(new); 否则，θt+1=θt，这里接受概率为αθt,θ(new)=min1,A，式中

④ 令t=t+1，重复进行步骤②、③，直到序列收敛。

2 数值模拟

运用R软件4.0.3版本进行数值模拟，评价基于现状数据建立广义极值回归模型和两类估计方法的有效性，针对极大似然估计和贝叶斯估计方法进行比较研究，讨论不同样本量N=300、500、1 000下循环500次的参数估计效果。

考虑协变量x～U2,10，生存时间T服从广义极值分布，可通过逆变换法生成T的随机数，u服从均匀分布U(0,1)，T从式(9)的模型中产生，即

(9)

设定一组参数为γ1=0.3、σ1=0.5、β1=1.5，超参数设置为a0=8、b0=8、c0=6、d0=9、μ0=0、σ0=10。删失变量C服从均匀分布，右删失比约为30%。最后生成一个示性变量v，当t小于等于C时v取1，否则取0。根据生成的随机数进行贝叶斯估计时，利用Gibbs算法嵌套MH算法进行抽样，模拟中设定的链长为8 000，预烧3 000次，利用剩下的5 000次样本进行参数估计。模拟结果见表1。

表1 x～U(2,10)的模拟结果

从表1可以看出，贝叶斯方法在此设定下拥有较小偏差和标准差，并且随样本量的增大，估计偏差、样本的标准差和理论标准差有逐渐减小的趋势，同时，样本的标准差和理论标准差越来越接近，CP也一直在0.95左右波动，模拟结果表明用贝叶斯方法进行参数估计效果很好。而用极大似然方法时，样本量较大时的估计效果与贝叶斯方法类似，但在样本量较小时估计效果不如贝叶斯方法，贝叶斯估计方法具有更小的样本标准差和理论标准差。因此，基于贝叶斯方法在估计广义极值回归模型参数时是有效的，并且估计效果优于极大似然估计。

敏感性分析通过改变先验中的超参数来进行额外模拟研究，以研究超参数的影响，如表2所示，得到的是敏感性分析结果，参数估计值和真实值都非常接近，参数估计的偏差、样本标准差和理论标准差都很小，经验覆盖概率在95%左右，说明在不同超参数设置下，通过贝叶斯方法得到的参数估计效果很好，结果表明所提出的方法是稳健的。同样，考虑在不同删失比下，对各个参数的估计进行分析，显示出在不同删失比下，各个参数的估计效果没有明显变化，且都显示出较小的偏差和标准差，详细结果见表3。

表2 敏感性分析：超参数为a0=10,b0=10,c0=6,d0=9,μ0=0,σ0=100的模拟结果

表3 X～tN(2,4,3,1)的模拟结果(γ=0.3,σ=0.5,β=1.5)

3 实证分析

本文将广义极值回归模型的贝叶斯估计应用于一个雄性RFM小鼠涉及肺肿瘤的致瘤性实验中，该数据来自1972年Hoel等[16]的研究，共收集144只小鼠的肺肿瘤数据。本文研究目的是比较不同环境(常规环境与无菌环境)对小鼠的肿瘤发病时间是否存在影响，看无菌环境是否可以提高小鼠总的生存率。将小鼠分为2组进行研究，一组为常规环境，另一组为无菌环境。这项研究中小鼠的肿瘤发病时间并不能直接观测到，只有小鼠的死亡时间(以天为单位)以及小鼠死亡时肿瘤是否存在是已知的。表4列出这144只小鼠所涉及的数据名称以及类型。表5中的数据包括小鼠肺肿瘤数据信息，但是没有小鼠肿瘤发病的精确时间，其中82只小鼠在死亡时没有出现肿瘤，所以这部分观测可以认为是右删失数据， 62只小鼠在死亡时肿瘤已经发病了，说明小鼠肿瘤的精确发病时间是在研究开始到小鼠死亡的这段时间内。例如，观测到的小鼠有肿瘤并且死亡时间为381 d，表示小鼠肿瘤的发病时间出现在381 d之前，但精确的时间未知。用区间删失I型时间数据，即现状数据来描述小鼠肿瘤发病时间。

表4 144只雄性RFM小鼠肺肿瘤数据的变量名称

表5 雄性RFM小鼠肺肿瘤数据

在进行数据分析过程中，假定小鼠肿瘤发病时间服从广义极值分布，在进行贝叶斯估计时参数的先验设定与模拟研究中相同，估计结果见表6。

表6 小鼠肺肿瘤数据分析结果

采用贝叶斯估计方法对小鼠肺肿瘤数据进行分析，通过表6可看出，3个参数都表现出对小鼠肿瘤的发作时间有显著效果，且均具有较小的标准差，其中形状参数γ决定了分布的尾部形状：当γ>0时，分布的尾部较长，广义极值分布趋向于Frechet分布；当γ=0时，分布的尾部呈指数状，这是广义极值分布的特殊情况，称为Gumbel分布；当γ<0时，分布具有有限的上端点，广义极值分布趋向于Weibull分布。由表6得到γ的估计值为0.416 0，所以广义极值分布化为极值Ⅱ型(Frechet)分布，并且协变量β对应的P值小于0.05，表明环境因素对生存概率的影响非常显著，同时画出广义极值回归模型的生存曲线如图1。

图1 小鼠的生存曲线

根据图1的生存曲线可以看出：不同环境对小鼠生存率存在显著影响，常规环境(x=1)可以提高小鼠总的生存率，意味着小鼠在常规环境下的生存概率比在无菌环境下的生存概率更高。该分析结果与Sun[17]关于小鼠肺肿瘤例子所得的结论基本一致。由此实例，可以证明贝叶斯估计方法对广义极值分布下区间删失I型寿命时间的回归模型的参数估计是有效果的，且具有较小的标准差。

4 结语

本文主要研究基于现状数据下的广义极值回归模型，利用R软件进行数值模拟，采用MCMC方法和MH算法对后验似然函数进行求解，同时比较极大似然估计方法和贝叶斯估计方法对广义极值回归模型的估计效果，结果表明贝叶斯估计的样本标准差和理论标准差更小，估计效果更精确，证明基于现状数据建立贝叶斯生存回归模型的有效性。在实证分析中，选取144只雄性RFM小鼠的肺肿瘤数据，从分析结果可知，生存环境对小鼠的生存概率存在一定影响，小鼠在常规环境下比在无菌环境下的生存概率更高，即处于常规环境下的小鼠具有更高的生存率。