基于模型平均的右删失数据下加速失效时间模型的统计推断研究

王淑影,张亚男,程云飞

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

一、引言

在经济学、生物医学、临床试验等研究领域中,由于各种原因经常会出现数据删失的情况。例如,在临床试验过程中,由于患者的异质性或者费用的限制等原因导致在试验结束时,某些个体还是没有发生感兴趣的事件。此时研究人员不能精确观测到某个个体的确切失效(死亡)时间,只能观测到生存时间大于删失时间,即右删失数据。针对右删失数据构建模型,最常用的就是比例风险(Proportional Hazards,PH)模型和加速失效时间(Accelerated Failure Time,AFT)模型。两者的不同之处在于PH模型假设对数危险率函数与协变量之间存在线性关系,而AFT模型直接对生存时间构建模型,研究协变量与对数生存时间的线性关系。但是,在不同实际问题中,PH模型的比例风险假设可能很难验证,因此本文直接考虑更具有解释性的AFT模型来对右删失数据进行建模。

在传统的AFT模型构建过程中,往往都是在假设协变量给定的情况下进行统计推断研究,忽略了此过程中的不确定性,导致估计量的方差被低估,进而得到不准确的预测结果。在这种情况下,有必要将传统的AFT模型与模型平均方法结合起来,避免忽略模型选择过程引入的不确定性。因此,本文考虑在右删失数据下,构建基于信息准则的AFT模型平均方法,并且基于估计和预测结合的双重视角来分析模型平均后AFT模型的拟合效果,以此来避免选择单一模型潜在的风险,减少估计的均方误差,提高覆盖概率,为以后更好地对AFT模型进行研究提供新的思路。

二、文献述评

加速失效时间(AFT)模型最早是由Pieruschka在1961年提出的,后来国内外专家学者对加速失效时间模型进行了大量的研究[1]。例如Buckley等将加速失效时间模型应用到了右删失数据下,并提出了经典的Buckley-James估计方法[2]。Tian等在现状和区间删失数据下,对加速失效时间模型提出了通过Wald-type来构造的一种新的参数估计方法[3]。Arnošt等针对任意删失数据,提出了一种半参数方法来对加速失效时间模型进行估计[4]。Li等考虑在带治愈组的右删失数据下,提出了一种半参数加速失效时间治愈模型,并针对该模型开发了一种EM算法[5]。Scolas等针对带治愈组的区间删失数据,提出了一个十分灵活的扩展广义Gamma加速失效时间模型[6]。现有文献大多考虑具有固定协变量的加速失效时间模型,但是这样做可能会丢失候选模型中包含的一些有用信息。克服这种问题的有效途径就是使用模型平均方法,通过对候选模型的估计量进行加权,可以有效地减少有用信息的遗失。

模型平均方法分为贝叶斯模型平均(Bayesian Model Averaging,BMA)和频率模型平均(Frequentist Model Averaging,FMA)两大类。许多早期工作从贝叶斯的角度讨论了模型平均方法,比如高磊等对BMA方法在一些难点问题上的理论进展进行了综述,并介绍了BMA方法在大数据背景下的应用前景[7]。Draper采用BMA方法进行预测研究并通过油价预测等两个实例证明了该方法的有效性[8]。但是,考虑到BMA方法存在模型假设十分复杂且难以直观解释的问题,FMA方法开始受到越来越多的关注和研究。例如,Zhang等提出了一个基于函数线性回归模型的交叉验证模型平均估计方法来预测响应变量[9]。Bai等针对期望回归模型进行了最优模型平均估计,并证明了所得模型平均估计量是渐近最优的[10]。Zhang等在固定参数设置的条件下,研究了Mallows模型平均估计和Jackknife模型平均估计下模型权重的渐近性质等[11]。这些工作为频率模型平均研究提供了深刻的理论支持和有效的方法。

文献梳理表明,国内外关于模型平均方法的研究主要聚焦于完整数据类型,对复杂删失数据的研究不多,且很少有文章研究加速失效时间模型的模型平均问题。因此无论在候选模型集合构建以及如何评价最终模型的好坏等方面都有待进一步深入研究。

三、数据、模型及似然函数

分布函数为:

密度函数为:

则对数似然函数为:

四、基于信息准则的AFT模型平均方法

在生物统计领域中,已有研究大多是先选定一个模型,然后对选定模型进行统计推断,这样可能会造成估计结果不稳健。而模型平均是解决这个问题的有效方法,模型平均方法通过考虑模型的不确定性和优化模型权重来最小化误差。

在本文中,我们考虑p+q个协变量,其中xi中包括p个确定被包含在模型中的协变量,zi中包括q个可能被包含在模型中的协变量,那么就会产生2q个候选模型。对第m个候选模型,假设为如下AFT模型:

均方根误差是一个能够较好地衡量估计精度的评价值,一般来说均方根误差值越小则估计精度越高[12]。

在实际情况中,受试者可能更关心的指标是生存概率,也就是在单位时间内生存的可能性大小。

对应累积生存率为0.5时的生存时间叫做中位生存时间。这是一个非常重要的指标,有时比平均生存时间重要[13]。

其中,r为T分布的50%分位点。

如果给定了一组数据,但是产生数据的真实模型未知时,需要找到所有可能的模型,这时就要通过比较研究哪个模型的拟合效果较好来决定,这个过程就是模型选择。在模型选择中,最常用的是基于Akaike信息准则(AIC)的模型选择和基于Bayesian信息准则(BIC)的模型选择,表达式分别如下所示:

AICm=-2logm+2lm,BICm=-2logm+lmlogn

通过模型选择可能得到存在有局部误差的模型,进而导致估计结果不准确,这时模型平均方法就会被考虑。模型平均方法最核心的问题就在于权重的选择,本文采用基于Smoothed AIC(S-AIC)和Smoothed BIC(S-BIC)的模型平均方法来计算组合权重,权重的表达式如下:

其中,m代表第m个模型,ωm是各个候选模型的权重,xIC表示AIC或BIC。设ω=(ω1,ω2,…,ωm)T是m个模型的权重向量,则:

对每个候选模型进行极大似然估计后会得到各个参数对应的估计值,则模型平均后的未知参数估计的表达形式如下:

模型平均后的生存概率的具体形式如下:

模型平均后的中位生存时间表达形式如下:

其中,r为T分布的50%分位点。

五、模拟研究

其中,i=1,2,…,n,β1、β2、γ1、γ2、γ3、γ4是回归系数,σ是尺度参数,εi是服从标准正态分布的误差项,μ是截距项。xi1、xi2是确定包含在模型中的协变量,zi1、zi2、zi3、zi4是可能包含在模型中的协变量,那么候选模型为24=16个。各个协变量由以下分布生成:xi1～N(0,1),xi2～N(0,1),zi1～N(0.5,1),zi2～U(0,1),zi3～P(1),zi4～B(1,0.5)。C服从均匀分布U(0,cl),其中cl控制删失率。本文考虑两组参数真值,设置A为窄模型(即只包含确定协变量的模型):β1=1,β2=1,γ1=0,γ2=0,γ3=0,γ4=0,μ=0.5,σ=0.4;设置B为宽模型(即包含确定协变量和所有可能协变量的模型):β1=-0.2,β2=0.2,γ1=0.2,γ2=-0.3,γ3=0.2,γ4=-0.2,μ=0.5,σ=0.4。删失比例为30%,模拟循环1 000次。

表1 设置A两组样本量下三个感兴趣指标的均方根误差

表2 设置B两组样本量下三个感兴趣指标的均方根误差

在表1和表2中,AIC代表基于AIC准则的模型选择方法,BIC代表基于BIC的模型选择方法,S-AIC代表基于Smoothed AIC准则的模型平均方法,S-BIC代表基于Smoothed BIC准则的模型平均方法(下同)。在设置A中,窄模型是真实模型。从理论上来说,BIC是一个一致的模型选择器,通常适用于具有少量参数的模型,因此在真实模型为窄模型的情况下,BIC优于AIC,S-BIC优于S-AIC。在指标(a)和指标(c)下,基于模型平均方法得到的结果与基于模型选择方法得到的结果相差较大,基于信息准则的AFT模型平均方法在参数数量较少时的估计精度更高。对于设置B,如果真实模型中的非零参数较多,S-BIC还是最好的方法。这可能是因为AIC假设真实模型是高维的,需要在有限的情况下尽可能多的参数来描述它。Sakamoto等指出,AIC不是估计真实模型的标准,而是最佳模型拟合的标准[15]。也就是说,AIC力求找到最佳的近似模型,如果数据较少模型的维数将会降低,随着可用信息的增多,模型的维数将会增加。相比之下,BIC是一致的。它假设真实模型存在且是低维的,可对真实模型进行一致性估计。

总的来说,无论真实模型为窄模型还是宽模型,基于Smoothed BIC准则的AFT模型平均方法的表现都是最好的。这主要是因为模型选择方法在使用过程中是将最终选择的模型视为提前选择,忽略了模型选择过程带来的额外不确定性,进而导致了较大的估计误差。随着样本量大小的增加,在局部错误指定的假设下,参数将会更接近真实值,所有候选模型都更接近真实模型,因此模型平均方法和模型选择方法的估计精度都在逐渐提高。

六、肺癌数据分析

使用本文提出的模型平均方法对一组肺癌数据进行研究,以此来进一步证明所提出方法的有效性。这组完整的数据集来自于R软件包survival中,是为比较接受两种不同治疗方案的肺癌患者生存时间而设计的一个临床试验。在这项研究中,137名肺癌患者被随机分为两组,分别接受两种不同的治疗方案,感兴趣的是肺癌患者的生存时间。数据集中有6个协变量,包括治疗方案(x1)、细胞类型(x2)、卡诺斯基评分(z1)、从诊断到治疗的时间(z2)、年龄(z3)以及是否有过既往治疗(z4)。我们感兴趣的问题是基于上述6个协变量,对数据进行预测。Kalbfleisch等人的研究结果表明,不同的治疗方案和细胞类型会对肺癌患者的生存时间产生显著影响[16]。所以,选定治疗方案和细胞类型为必选协变量,其他4个协变量则为候选协变量,这样共有24=16个不同的候选模型。

为了比较AIC、BIC、S-AIC和S-BIC四种方法的效果,随机地将数据分为两部分,分别记为训练集D(0)和测试集D(1)[17]。训练集的样本量n0依次选取了80,90,100,110和120,而测试集的样本量大小为n-n0,其中n=137是整体的样本量。首先对n0个样本进行极大似然估计,其次对n-n0个样本进行预测,并计算样本的均方预测误差(Mean Squared Prediction Error,MSPE),即:

从表3中MSPE的均值和中位数的角度来看,最优的方法是S-AIC,它比其他三种方法的MSPE的均值和中位数都要小。S-AIC和S-BIC方法的MSPE的均值和中位数都比相应的AIC和BIC方法要小,这说明基于模型平均方法的AFT模型预测结果比模型选择的预测结果更准确。此外S-AIC方法略优于S-BIC方法,AIC方法略优于BIC方法。尽管表3中四个方法的MSPE的均值和中位数有较大的不同,但是可以清楚地分辨出各个方法的表现。为了检验所提方法得出的结论,进一步考虑了Diebold-Mariano检验,通过R软件包forecast中的dm.test函数计算出两种方法构造的模型之间的DM统计量和p值,以此来观察两种模型预测能力的优劣性[18]。DM检验的原假设为H0:DM统计量均值为0,即两种方法构造的模型预测能力相同;备择假设为H1:DM统计量的均值不为0,也就是说两种方法构造的模型具有不同的预测能力。正的DM值表明分母对应的方法优于分子对应的方法,若DM检验的p值小于0.05,拒绝原假设,那么就说明两种方法构造的模型具备不同的预测能力。

表3 肺癌患者生存时间的MSPE

表4 肺癌患者生存时间MSPE的DM统计量

在表4的前3列中,我们可以看到DM的值都是负的,且p值几乎为0,这表明S-AIC方法显著优于其他三种方法。第4列和第6列中的DM值都是正的,且p值几乎为0,说明S-BIC方法显著优于AIC和BIC方法。第5列中的DM值都是负的,且p值几乎为0,说明AIC方法优于BIC方法。这些结论与表3中的结论是一致的。

七、总结及讨论

本文主要考虑右删失数据下,提出基于AFT模型的频率模型平均方法,使用极大似然方法进行参数估计及预测,其中误差项服从标准正态分布。模拟研究表明,所提方法的表现更好。模型平均方法考虑了所有协变量对响应变量的影响,而模型选择是只用选择出来的部分协变量来解释响应变量。由于模型平均方法是根据协变量的个数来构建候选模型,所以当协变量比较多的时候,具体实现计算就会变得十分复杂[19]。近些年来频率模型平均方法快速发展,提出了各种最优模型平均的方法,但本文仅考虑了基于信息准则的模型平均方法。如果使用渐近最优的模型平均方法,如Mallows模型平均方法、Jackknife模型平均方法、OPT模型平均方法等,也许会得到更好的估计结果。这些问题值得我们进一步研究。