唐璐薇,吕 超,文秋月,韦程东

(1.广西科技师范学院 数学与计算机科学学院,广西 来宾 546100; 2.柳州市妇幼保健院,广西 柳州 545000; 3.南宁师范大学 数学与统计学院,广西 南宁 530000)

0 引 言

生存分析是一种既考虑结果又考虑生存时间的统计方法,其可充分利用截尾数据所提供的不完全信息,对生存时间的分布特征进行描述,以及对影响生存时间的主要因素进行分析.比例风险模型是研究生存数据最广泛的统计半参数模型之一[1-2],但其对数据的拟合效果并不理想.因此,需要我们寻求一个比例风险模型的替代模型,即加性风险模型.与比例风险模型不同的是,加性风险模型是假设基底风险函数与协变量之间的一个加性结构.在实际应用中,加性风险模型对数据的拟合效果更好,正是由于它这一特性,加性风险模型中的回归参数更容易解释实际意义[3-4].

目前,运用不同模型对数据进行拟合的研究有很多.顾刘金首先利用Cox模型对数据进行拟合,但由于Cox模型具有比例风险性,故基于SPSS采用Logistic回归模型对数据进行拟合得出结果[5].Jardim等利用微分方程模型对葡萄牙的COVID-19相关数据进行拟合分析,同时指出该方法可针对不同的流行病学数据进行拟合分析[6].此外,很多学者也运用B样条拟合方法进行研究.曾卓等提出了基于三次B样条小波的曲线拟合方法,该方法行之有效,为曲线拟合提供了一种兼顾拟合精度、光顺性与加工精度的方法[7].张永华等针对传统几何轨迹跟踪算法切向角获取依赖高精度惯导设备的问题,提出了基于三次B样条曲线拟合的轨迹跟踪算法,该算法有效解决了传统算法的问题[8].Bi等基于全局拟合法的缺点,提出一种通用、快速的B样条误差拟合方案,该方法与传统的拟合方法相比,显著地提高了工作效率[9].庞宇等提出一种结合幅度系数法和波形特征法对脉搏波采用三次B样条曲线拟合的血压检测算法,该算法能够提高特征点的准确率[10].芦穗豪等提出一种基于改进麻雀搜索算法的B样条曲线拟合方法,旨在利用最少控制点高效地达到曲线拟合的目标精度,进而提升传统建模方法的精度和效率[11].徐超清等基于神经纤维走向信息考量的问题,提出一种基于B样条拟合与回归模型的脑神经纤维聚类方法,该方法在功能区层面的聚类可以更有效地分割出具有解剖学结构的脑神经纤维[12].

上述研究大多是运用传统的模型对生存数据进行分析,且运用三次B样条平滑函数对生存数据进行拟合,而利用加性风险模型进行相关研究的较少.本文运用三次B样条拟合方法,研究经药物治疗后白血病患者的生存周期与实际情况的拟合程度,从而判断药物治疗的有效性,并运用单变量加性风险模型和多变量加性风险模型找出对生存时间影响显著的协变量,作出生存时间的预测图,从而为决策者提供一定的参考依据.

1 模型介绍

1.1 加性风险模型

假设研究队列数据包含N个独立的样本,设T为失效时间,C为相应的截尾时间,Z(t)(0≤t≤τ)为协变量过程的向量,其中τ<∞表示后续时间.由于数据右删失,所以只能观察X和δ,其中X=min(T,C′)和δ=I(T≤C).用Yi(t)=I(Ti≥t)表示历险过程的示性函数,Ni(t)=ΔiI(Ti≤t)表示计数过程,这里I(·)是示性函数.假定T和C是条件独立的,可考虑以下加性风险模型.当给定协变量Zi(t)时,失效时间T的风险函数为:

(1)

其中,λ0(t)为一个未指定的基底风险函数,β=(β1,β2,…,βp)为未知回归参数的p维向量.

1.2 三次B样条插值函数

B样条曲线是在Bezier曲线基础上发展起来的,其克服了Bezier曲线整体控制性所带来的不方便.它是通过逼近一组控制点生成的曲线,计算公式为:

(2)

其中,pk为输入的一组数据中的第k个数据点;Nk,d为B样条混合函数,这里为加性风险模型,k为第k个混合函数,d为次数.本文采用三次B样条插值函数进行平滑,即d=3.

B样条曲线保留了Bezier曲线的优势,同时曲线在拼接时又比Bezier曲线方便,在修改时可做局部修改,而且在调整某一控制点时不会影响整条曲线的趋势.其性质主要有以下几点:

(1)局部性.K阶B样条曲线上的一点至多与k个控制顶点有关,与其他控制顶点无关.

(2)几何不变形.B样条曲线的形状和位置与坐标系的选择无关,不管坐标系如何变化,B样条曲线的形状仍保持原样.

(3)凸包性.与Bezier曲线一样,B样条曲线落在Pi构成的凸包中,其凸包区域小于或等于同一组控制顶点定义的Bezier曲线凸包区域.

2 白血病数据的三次B样条平滑

常见的生存数据处理模型有比例风险模型(Cox模型)、加性风险模型等,二者在处理生存数据上都有一定的优势.下面给出Cox模型和加性风险模型对白血病数据分析的结果.

(1)Cox模型数据拟合.根据已知白血病数据中的自变量,选取6个自变量作为协变量:age(年龄)、sex(性别)、 ph.ecog(ECOG评分)、ph.karno(医师的Karnofsky评分)、meal.cal(用餐时消耗的卡路里)、wt.loss(最近6个月的体重减轻),并运用Cox模型对其进行分析,结果见表1.

表1 Cox模型对6个协变量的分析

由表1可知,性别和ECOG评分与生存时间显著相关,且P值远小于0.05,模型拟合效果较好.

(2)加性风险模型数据拟合.为与Cox模型拟合效果进行对比,下面运用加性风险模型对上述6个协变量进行分析,查看加性风险模型对数据的拟合效果及相关性,结果见表2.

表2 加性风险模型对6个协变量的分析

由表2可知,加性风险模型对数据的拟合效果与Cox模型一样,性别和ECOG评分与生存时间显著相关,且P值远小于0.05,模型的拟合效果较好.但加性风险模型还给出了对生存时间的整体方差解释率,更进一步说明了该模型对数据的拟合效果.因此,本文在三次B样条平滑函数下,运用加性风险模型来研究白血病的治疗数据,从而给出相关结论.

2.1 三次B样条平滑

为减轻协变量对术后患者生存时间的影响,本文利用三次B样条平滑法对受干扰的协变量进行平滑,以消除协变量对响应变量即生存时间的影响.下面首先讨论利用三次B样条平滑法对受干扰的加性风险模型中的协变量进行平滑,然后讨论单变量加性风险模型和多变量加性风险模型中协变量与生存时间的关系.三次B样条平滑过程分以下几点讨论:

(1)59岁以下的男性,其生存时间与用餐时消耗的卡路里、最近6个月体重减轻的关系见图1.其中,第一行指定的平滑系数分别为k=3,k=5;第二行指定的平滑系数分别为k=5,k=7.由图1可知,从整体看,白血病患者的生存时间随着用餐时消耗卡路里的增加而增加,但不同平滑程度的选择影响了局部趋势的解释,特别是不同年龄生存时间的变化趋势是不同的,且最近6个月的体重减轻对生存时间的影响也较明显.

图1 59岁以下男性两协变量与生存时间的关系

(2)59岁以上的男性,其余条件与(1)相同,得到的结果见图2.由图2可知,59岁以上的男性,随着年龄的增加,其消耗的卡路里逐渐降低,而相对(1)来说,近6个月的体重减轻呈细微的增加趋势.

(3)59岁以下的女性,其余条件与(1)相同,结果见图3.由图3可知,59岁以下的女性,随着年龄的增加,其消耗的卡路里趋于平稳,而随着近6个月体重减轻的越来越多,生存时间呈先下降后上增的趋势.

图3 59岁以下女性两协变量与生存时间的关系

(4)59岁以上的女性,其余条件与(1)相同,结果见图4.由图4可知,59岁以上的女性,从整体看与图3中59岁以下的女性分析结果相似,但在平滑系数k=7时,脂肪消耗量发生细微变化,随着脂肪消耗的增多,生存时间呈细微上增的趋势.

通过三次B样条平滑对年龄(age)、性别(sex)、ECOG评分(ph.ecog,0=好,5=死)、医师的Karnofsky评分(ph.karno,0=差,100=好)、用餐时消耗的卡路里(meal.cal)和最近6个月的体重减轻(wt.loss)等6个协变量中受干扰的协变量进行平滑后,暂时不能确定哪些协变量对生存时间的影响显著.下面研究生存时间的单变量加性风险模型,并根据研究结果选择协变量进行多变量的加性风险模型研究,以探索这些协变量对生存时间的影响.

2.2 单变量加性风险模型研究

选择协变量:年龄(age)、性别(sex)、ECOG评分(ph.ecog,0=好,5=死)、医师的Karnofsky评分(ph.karno,0=差,100=好)、用餐时消耗的卡路里(meal.cal)和最近6个月的体重减轻(wt.loss),并用这些协变量拟合加性风险模型,结果见表3.其中,s()为运用三次B样条平滑后的符号表示.

表3 单变量加性风险模型拟合结果

由表3可知,性别、年龄、ph.karno和ph.ecog变量具有较好的统计学意义(P<0.05),说明它们对生存时间的影响是显著的.此外,年龄和ph.ecog具有正β系数,而性别和ph.karno具有负β系数.因此,年龄较大和ph.ecog较高与事件发生率呈正相关,而女性(sex=2)和ph.karno则与事件发生率呈负相关,即年龄和ph.ecog是死亡的危险因素,女性性别和医师的Karnofsky评分是死亡的保护因素.

2.3 多变量加性风险模型研究

为研究性别、年龄、ph.ecog和ph.karno如何共同影响生存时间,本文利用多变量加性风险模型进行分析,结果见表4.

表4 多变量加性风险模型拟合结果

由表4可知,加性风险模型对生存时间的整体方差解释率达78.5%,拟合效果较好.在多元加性风险模型分析中,协变量性别和ph.ecog保持着显著性(P<0.05),但协变量年龄和ph.karno不显著(年龄:P=0.171 226,ph.karno:P=0.186 82均大于0.05).

性别的P值为0.000 712,风险比HR=exp(coef)=0.56,表明患者的性别与死亡风险降低之间有很强的关系.在保持其他协变量不变的前提条件下,女性(sex=2)相比男性,其死亡风险低44%.可见,性别为女性与良好的预后相关.

同样,ph.ecog的P值为0.000 323,风险比HR=1.88,表明ph.ecog值与死亡风险增加之间有很强的关系.相比之下,年龄的P值为0.171 226,风险比HR=exp(coef)=1.01,95%置信区间为0.99～1.03.由于HR的置信区间为1,因此该结果表明,在调整ph.ecog值和患者性别后,年龄对HR差异的贡献较小,且不显著.

为研究白血病患者治疗的生存数据与治疗的关系,以及白血病治疗的有效性,针对上述对协变量的研究,本文选取性别和ph.ecog作为协变量来拟合加性风险模型.使用K-M估计方法,以存活天数为横坐标,以生存率为纵坐标,给出估计的生存率和生存率的置信区间,见图5.

图5 存活天数与生存率的关系

如图5所示,随着存活天数的增多,术后存活率也逐渐下降,这与自然规律相符.此外,由上述对多变量加性风险模型协变量的分析可知,性别和ph.ecog对模型具有显著影响.性别和ph.ecog对生存时间和生存率的影响见图6.

图6 性别和ph.ecog对生存时间的影响

由图6可知,女性的整体生存率高于男性,在ph.ecog评分表中对应的评分为5,且生存率随着生存时间的变化逐渐下降.但这种差异是否显著并不确定,还需要进行统计检验,检验结果见表5.

表5 性别和ph.ecog差异显著性检验

由表5可知,将性别和ph.ecog综合考虑,得到的P值远小于0.05,说明男女之间以及不同的ECOG评分,其生存率是有差异的.

与K-M生存曲线不同的是,加性风险模型拟合曲线得到的生存率是在运用三次B样条平滑其他协变量后所预测的生存率,并不是实际观察到的生存率.加性风险模型拟合生存率曲线见图7.

图7 加性风险模型拟合生存率曲线

由图7可知,预测患者的生存率随着生存时间的增多而逐渐下降,与使用K-M估计方法对白血病患者的生存时间和生存率所作出的曲线趋势一致.由于本文研究的是白血病治疗数据,涉及对医疗效果的评价,因此选择精确度较高,且对解决高度非线性预测问题有突出能力的加性风险模型来对数据进行拟合是可行的.

3 结 论

本文研究的白血病数据涉及分析白血病药物治疗的效果,精度要求较高,需要一个能够对生存数据进行分析且拟合效果较好的模型.本文利用加性风险模型对其进行拟合研究,并通过图形展示患者生存时间的变化.为避免数据缺失给加性模型研究带来误差,利用三次B样条平滑函数对其进行平滑,并对删失的数据进行处理,以避免造成对已有数据的影响.通过对数据分析可知,不同平滑程度对拟合效果影响不大,不能直接表达哪些协变量对生存时间影响显著.本文不足之处在于运用三次B样条平滑函数对加性风险模型进行平滑后,只对白血病治疗数据进行研究,没有考虑到其他疾病数据的拟合效果是否一样.此外,针对白血病的治疗数据也只运用了加性风险模型对其进行研究,没有运用其他模型进行拟合.今后的研究将针对这两个问题进行深入探讨,从多方面进行生存分析研究.