第五节 生存分析
生存的直接含意是与死亡结局相对立的结局。但这里“生存”是一个广义的概念,泛指二分类结局中的一种。如可以把发病/不发病这一二分类结局中的不发病看成为“生存”,发病看成为“死亡”。把疾病复发/缓解这一二分类结局中的缓解看成为“生存”,复发看成为“死亡”等。在生存分析中同时考虑两个反应变量:即生存时间的长短和生存状态。以临床治疗效果来说,在同一时间长度情况下,生存概率大者治疗效果为优。
率这个词的中文含义既包含概率的意思,也包含速率的意思。但在概率论中,概率、速率是相关而不相同的两个概念。习惯上人们所称谓的生存率,实际上是生存概率,它与前面介绍的以人年为分母的发病率在数学性质上是不相同的。为避免混淆,在本文中一律用生存概率这一名词,而不用通俗的生存率。
在生存研究中经常有观察对象中途退出研究的现象,即该观察对象在尚未到达观察终点即“死亡”前就终止了观察。对于这类失访的观察对象所提供的时间信息是不完全的信息,因为其未到达研究终点。在统计学上称这类信息为“截尾”数据或“删失”数据(censored data)。如何利用好这类“截尾”数据也是生存分析中所要考虑的问题。
为了便于分析,通常假定截尾为无信息截尾(noninformative censoring)。也就是说,在整个研究期间,个体的截尾原因应与事件的发生无关,不能提供任何与模型参数相关的信息。遗憾的是,在大量的癌症患者生存研究中,无信息截尾通常是不存在的。例如,在实际研究中,截尾的发生往往与肿瘤药物的副作用、患者本身的身体状况等密切相关。因此,需要对是否存在信息截尾进行评估。然而,信息截尾的存在及其对生存分析的影响等评估相当困难。
在各类截尾情形中,失访是最有可能与事件发生存在着关联性的截尾情形。这往往可以从专业的角度进行定性判断。一种较为客观实用的评估信息截尾是否存在的方法是,考察整个研究期间的失访比例(loss to follow-up)。如果失访比例较大,或者,对每一类别分别作Kaplan-Meier生存曲线图,观察每个时段失访比例。如果存在较大差异,那么,就有理由认为失访是信息截尾的。另一种较为简单实用的方法是,通过事后的截尾敏感性分析,如,最佳-最差情形法(best case-worst case scenario),来评估截尾对生存时间的影响大小。例如,第一次分析对截尾个体按照无信息截尾进行处理,第二次分析直接将截尾个体的截尾时间当作其生存时间进行处理,比较这两次分析的结果,如若存在着较大差异,就可以认为截尾是有信息的。这种方法对于截尾较少的情形较为有效。需要注意的是,在这种情形下,信息截尾对生存分析的影响通常也较小。对于存在着较大截尾比例的研究,往往认为该研究的质量不高,其结果的可靠性与意义并不大。因此,在实际研究中,如何确保截尾尽可能少是癌症患者生存研究的关键性问题。
1. Kaplan-Meier乘积极限法[Kaplan-Meier method,product-limit(P-L)method]
当资料不含有截尾观察数据时,生存概率的估计方法十分简单,即生存人数除以期初观察人数。但当资料含有截尾观察数据时,必须将截尾数据所提供的部分信息考虑进来。Kaplan-Meier法是Kaplan和Meier于1958年首先提出,根据概率论中关于总概率是条件概率连乘积的原理来计算生存概率,是生存概率的一种非参数估计方法。适用于观察例数较少的情况。
设有n个观察时间t
i,i =1,…,n,第i个观察若为完整时间则记为t
i,若为截尾时间则记为
,将n个时间分为m组(m≤n),记为:
t1≤t2≤…≤…≤tm
如果
= t
i则分为同一组,并令t
0=0和t
m +1=∞。在第j(j =1,…,m)组中的观察人数为n
j,死亡数为d
j,在已生存到时间t
j -1的基础上再能生存到时间t
j的条件生存概率的计算为:
则从t 0开始一直生存到时间t k的生存概率P(t k)的计算公式为:
又可将 P(t k)看成为生存时间等于或大于t k的概率。
生存率的标准误计算公式为:
生存概率 P(t k)的95%置信区间的计算公式为:
总生存概率P(tk)±1.96×se[P(tk)]
(11-65)
例11-14 48例急性淋巴细胞性白血病病人从治疗开始至死亡的时间(月)如下:
1,1,1,2,2,2,2,3,3,3,3,3,3,5,5,5,5,4,5,5,7,7,7,7,7,7,7,7,8,8,8,8,9,10,10,13 +,13,13,13,14,15,18,18,20,20,21 +,21,23
注:有+号者为截尾时间
其生存概率的计算过程列于表11-27中。
表11-27 48例急性淋巴细胞性白血病病人生存概率的Kaplan-Meier法估计结果
用Kaplan-Meier法估计的总生存概率是一种阶梯形概率,用估计的总生存概率所绘制的图形见图11-13。
图11-13 48例急性淋巴细胞性白血病病人生存概率的阶梯曲线
从图中可见,在10个月之前的曲线下降的坡度较陡,反映死亡风险较大,10个月之后的曲线下降的坡度较平坦,反映死亡风险较小。50%生存概率在开始治疗后5~7个月之间。
2.人寿保险法(actuarial method)
人寿保险法又称精算法,也是一种非参数生存概率估计方法,适用于观察例数较多的情况。
用人寿保险法计算生存概率的第一步是将观察人群按生存时间长短分为数个相等的时间区间。其次是计数在每一个区间开始时的观察人数和在区间内的死亡人数;第三步是计算每一个区间内的条件死亡概率和条件生存概率。最后按概率乘法原理计算总生存概率。
例11-15在5年期间内共随访了411例结肠癌病例,最长观察期为5年。为计算这一组病人的不同时间的生存概率,将总观察时间按每6个月1段共分为10个等长区间。记数每一区间开始时的病人数 n i、在区间内的死亡数 d i和截尾数c i。见表11-28中的第1 至4列。
表11-28 411例结肠癌病例按人寿保险法计算生存概率
续表
表中各列的说明:
第2列,区间(月)(t i- t i +1):t i是月初起点,t i +1是月末终点;
第3列,期初人数n i:n i +1= n i- d i- c i;
第4列,期内死亡人数 d i:在区间( t i- t i +1)内死于结肠癌的病例数;
第5列,期内失访人数c i:由于迁居或死于其他原因等未观察到其死于结肠癌的病例数;
第6列,校正期初人数N i:N i= n i- c i/2,即假定每名失访者在区间中点失访,在此区间内提供了半个区间的的观察时间;
第7列,条件死亡概率q i:q i= d i/N i,表示在区间起点t i生存的一个人死于区间( t i- t i + 1)的概率;
第8列,条件生存概率s i:s i= 1 - q i,表示在区间起点t i-生存的一个人能活过区间( t i- t i +1)的概率;
第9列,生存概率
P i:生存概率
,表示从观察起点开始能一直活过区间i的概率。即通称的生存率。
从生存概率 P i分析,前5年的生存概率下降很快,从0.647下降到0.392;后5年的生存概率比较稳定,从0.361下降到0.319。存活5年的概率为0.392。50%存活时间为:
3. Cox比例风险回归
随着现代医学的不断发展,有关肿瘤的随访研究越来越多。当前肿瘤尚缺少根治方法,这就使得怎样延长病人生存时间和提高生存质量的研究,显得越来越重要。近年来,怎样从众多的风险因素中分辨出对疾病的发生、发展有较大影响的重要因素,已成为肿瘤研究的一个热点。大量医学实例已经表明,肿瘤患者的生存时间与个体的性别、年龄、肝功能、肾功能等有着密切关系。不同的个体具有不同的属性与指标,因而,导致各个个体之间存在着较大的差异。如果已知在某一时刻、某一同龄的肿瘤患者群中发生一个死亡,那么,这一死亡的可能性对这群肿瘤患者中每一个个体而言,是不一样的。换句话说,就是在死亡面前,机会并不均等。这样,利用前面介绍的乘积极限法和人寿保险表法来估计死亡概率等非参数生存分析统计方法,就不能满足多因素分析的要求,必须要有新的模型来分析与研究各个个体背景因素不一致的情形下的生存状况。癌症患者生存资料中,反应变量为生存时间与事件结局,不是单一变量。然而,普通的多元线性回归和logistic回归难于分析此类资料,这是由于这两种模型通常不能全面利用这种存在截尾的不完全数据信息,使得所建立的模型失去统计效能。
目前,生存分析领域最重要的一类分析方法就是时间到事件(time-to-event)的分析方法。其相应提出的模型也被称之为生存模型。生存模型主要分为两大类,即20世纪六七十年代发展起来的参数模型(parametric model)与其同期发展起来的半参数模型(semiparametric model)。生存分析中的参数模型,必须事先指定误差项的基准分布,如指数分布,Gamma分布,Weibull分布等。与参数模型相比较,虽然半参数模型不能估计出各时点的风险率,但是,半参数模型对生存时间分布没有任何先验要求,并且,可以估计出各研究因素对风险率或生存时间的影响,鉴于肿瘤疾病发生、发展和死亡的复杂本质,很难确定分布的形式,因此,在肿瘤研究与应用领域中,半参数模型比参数模型更为可取,应用范围也更为广泛。
经典的半参数生存分析模型,主要包括两类,即一是Cox比例风险回归模型(proportional hazards regression model,Cox model),二是半参数加速失效时间模型(semi-parametric accelerated failure time model,AFT model)。本节先介绍Cox比例风险回归模型,下一节介绍加速失效时间模型。
英国统计学家D. R. Cox于1972年提出了一个巧妙的、用于分析生存时间的半参数模型,即,Cox比例风险模型(Cox's proportional hazard model),现在通常称为Cox模型。有关Cox比例风险回归模型的研究和应用,特别是在肿瘤的预后因素分析中,已取得长足的发展。
肿瘤患者的生存时间与一些外部或内部的因素有着相当强的关联,通常将这些因素称为协变量。协变量可以是多维的,记为X =(x 1,x 2,…,x p)。记每个个体协变量X i= (x i1,x i2,…,x ip),以表示与第i个个体有关的p维协变量,第i个个体的生存时间分布依赖于X i。一般地,将生存分布函数记为 S( t | X),其密度函数记为 f( t | X),风险函数记为h(t | X),以表示对协变量X的依赖关系。
(1)风险函数:
对于Cox比例风险模型而言,风险函数的理解是至关重要的。风险函数基本定义为:
式中的X表示可能对患者生存时间产生影响的各种因素,亦称为协变量(covariate)。通常假定这些变量不随时间的变化而变化。 t表示生存时间, h( t, X)称为具有协变量 X的个体在t时刻的风险函数(hazard function)。
风险函数具有三个基本特性:
1)风险函数不是概率函数,这是因为风险可能会超过1.0。这可以通过风险函数的定义直接得出。因此,风险函数不能按照概率予以解释。另外,虽然风险函数没有上限,但是,它不可能低于0。
2)由于风险函数本身是由条件概率函数的导数来定义的,它是一个不可观测的理论变量,因而,只能通过所获取的数据来估计风险函数,只能从估计的角度来予以解释。
3)风险最好从个体的角度来予以解释,而不能从群体的角度来理解。这是因为每个个体的风险函数是互不相同的。
风险函数 h(t,X)表示生存时间已达t时刻的个体在t时刻的瞬时风险率。它是针对个体而言的。风险函数意味着某个体在某个时间区间内,某事件的发生频数。也就是说,风险函数是有时间单位的,它是一个率的指标。例如,按照以往月统计数据表明,某个体在某个特定时间点患上呼吸道感染的风险为0.01,这就是说,在一个月内该个体患上呼吸道感染的次数为0.01次;如果以年来统计,那么,某个个体在某个特定时间点患上呼吸道感染的风险为1.0,这就意味着,在一年内该个体患上呼吸道感染的次数为1.0次。当然,这里必须假定患上呼吸道感染的风险在一月或一年内是恒定的。如果在一年内患病次数是恒定的,那么,最好采用以年为单位的风险函数。相对而言,一年内患病次数的假定要比一月内患病次数的假定更为严格,这是因为患上呼吸道感染的次数往往因季节而变动。因此,在对风险函数进行解释之前,必须对统计区间予以统一规范,并加以检验,而这通常从医学的专业角度来进行事前考虑。
对于上呼吸道感染这类易患易愈,易于确诊,无后效的疾病来说,基于个体的统计是易于描述与分析的。然而,对于癌症这种罕有的、难以确诊的疾病来说,基于个体的统计几乎是无法实施的,只能假定每个体的基准风险是相似的,也就是说,如果任意两个体基本条件一致的话,那么,他们所经历的风险应是完全一致的。在这种假定下,可以通过基于群体的统计来估计个体的风险,这就使得风险的解释具有群体特性(流行病学特征)。对于癌症之类的疾病,往往通过观察大量的随机样本来确定风险函数。例如,随机抽取10 000人,观察1年,假定未发生截尾,暴露共计10 000人年,在这一年中,共有10人发生某种特定癌症,那么,根据风险相似的假定,风险的最优估计为10/10 000 =0.0001。
风险相似的假定比较难以理解。现在给出一个常识性例子,以帮助理解。例如,当你说“这两辆不同的车都在以每小时80公里的时速在行驶”这句话时,其客观实现是,作为观察者的你,主观地认为如果这两辆车保持当前的行驶状态,那么,无须考虑它们之间是否确实存在着本质上的差异,而导致某个时段这辆车快一点,另外一个时段那辆车快一点,并不影响到你观察它们的整个行驶过程,这两辆不同的车都应在一个小时内行驶完80公里。
条件风险是指各种外在条件会影响事件发生的风险。这也就是说,对于某个特定个体的某个特定事件发生风险是会随着各种外在因素的变化而变化的,并非一成不变。例如,乙肝患者患肝癌的风险较大,而生活规律会降低肝癌患病的风险。大量事实表明,随着外在条件的剧烈改变,风险的改变往往也是阶梯样的,并不是连续性的,而那些大量的只能引起风险轻微改变的外在条件,往往是会相互抵消的。因此,可以建立与常规多元线性回归类似的统计模型,其回归系数解释大致同样类似于多元线性回归。
(2)Cox比例风险回归模型结构:
在比例风险模型中,假设在时点 t时,个体出现观察结局的风险大小可以分解为两个部分。第一部分为一个基准风险函数h 0(t)(baseline hazard function),表示某类个体的共性风险,该风险量是未知的,属于非参数部分,需要注意的是,基准风险函数h 0(t)会随着时间的变化而变化;第二部分为第j个影响因素使得该风险量从h 0(t)增加e βjxj倍,表示某类个体的特殊性质,属于参数部分,其值是可估计的。从而,个体在时点 t的风险量变成 h 0( t) e βjxj。因此,如果在 p个因素同时影响生存过程的情况下,在时点t的风险量(常称为风险比hazard rate,HR),则模型为
h(t,X)= h0(t)eβ1x1eβ2x2…eβpxp
因此,Cox比例风险模型的基本结构如下:
h(t,X)= h0(t)eβ1x1eβ2x2…eβpxp= eβ1x1+β2x2+…+βpxp
(11-66)
将基础风险移到等式左侧,两边取自然对数,得,
ln(Rh(t))= ln(h(t,X)/h0(t))=β1x1+β2x2+…+βpxp
(11-67)
式中R h(t)= h(t,X)/h 0(t)称相对风险。上式和多元线性回归模型非常类似,因此,Cox比例风险模型,也被称为Cox回归。与一般多元线性回归不同的是,Cox比例风险模型的截距项不是恒定的,会随着时间的变化而变化。Cox模型是目前解决多因素对生存过程影响最常用的统计分析方法,它将协变量对生存期的影响表现在风险函数的关系上,从而有效地解决了截尾数据的问题。
(3)Cox比例风险模型的基本假定:
根据Cox比例风险模型基本形式,该模型要求资料事先满足两个假定:一是,对数线性假定;二是,比例风险假定(assumption of proportional hazard)。
1)对数线性假定:
Cox比例风险模型假定协变量的影响为线性模型,风险比 R h(t)= h(t,X)/h 0(t)的自然对数与各影响因素呈线性关系,服从线性模型的一般规则。各风险因素对风险比的影响具有乘积性,而不是可加性。
2)Cox比例风险假定:
比例风险假定,是指假设有两个个体,其协变量的值分别为X 和 X *,其风险比为:
该比值与基础风险量h 0(t)无关,在时间t上为常数。也就是说,风险比是成正比例关系。上式的值被称为具有风险因素X的个体对风险因素为X *的个体的相对危险度(relative risk,RR)或风险比(risk ratio,RR)。这种协变量效应不会随着时间的变化而改变的假定,被称为等比例风险假定,简称PH假定。这也是比例风险模型的由来。该假定暗示各组的风险曲线是成比例的,不能出现交叉的情形。对于Cox模型来说,比例风险假定是至关重要的。如若这一假定不成立,那么,Cox模型就会成为一种统计效能极差的生存分析模型。因此,需要对比例风险假定应作出检验与评估。
(4)参数估计:
除了对数线性假定与比例风险假定,还必须无信息截尾假定,以及,在协变量X给定后,事件发生的时间和截尾时间之间相互独立。为此,Cox比例风险模型提供了一种新的参数估计方法,即偏似然估计(partial maximum likelihood)。通常,似然函数是基于结果变量的分布而来的,而Cox比例风险模型对结果变量(生存时间)没有任何假定分布,所以,该模型不能像参数模型那样建立一个完全基于结果变量分布的似然函数。相反,Cox似然函数的建立是基于事件发生的秩序,而不是事件的联合分布。因此,Cox似然被称为偏似然。偏似然方法基于如下假设:两个生存时间(或事件发生时间)的间隔长度对协变量和风险比之间的关系无任何影响。也就是说,两个连续事件的间隔长度无论是否为0,对偏似然函数不能提供任何有效信息,仅仅考虑事件发生的先后顺序,故而,在Cox比例风险模型中,事件发生时间,也被称为有序事件发生时间(ordered failure times)。由于Cox比例风险模型仅仅使用了部分变量资料,没有估计基准风险函数,Cox的似然函数被称为偏似然函数(partial likelihood function)。偏似然函数不是真正的似然函数,这是因为截尾和未截尾的实际生存时间没有在偏似然函数中得到具体体现。这就使得偏似然法难以克服巨大的理论研究难题。
(5)模型参数的意义及其解释 1)回归系数与相对危险度:
由Cox比例风险模型的公式,可以得到:
两边取自然对数,可得到
因而,回归系数β j不会随着时间的变化而改变,保持恒定。其与风险函数h(t,X)之间有如下关系:当β j>0时,则随着协变量x j的增大,其风险h(t,X)也相应增大,表示患者的风险越大;当β j<0时,则随着协变量x j绝对值的增大,其风险h(t,X)也相应减小,表示患者的风险越小;当β j=0时,则随着协变量x j的增大对其风险h(t,X)没有影响,表示患者的风险无变化。
对于具有协变量X和X *的两个个体,其风险比为:
该比值同样不会随着时间的变化而改变,它表示在任何生存时间上,二者之间的相对危险度。对某个协变量特定取值而言,在保持其他协变量取值不变的情形下,回归系数β j可以解释为该协变量变化一个单位,其相对危险度变化e βj。据此,可以求得各协变量的相对危险度估计值,这就使得Cox比例风险模型具有明确的流行病学意义,即,当协变量x改变一个单位时,引起的死亡风险改变倍数的自然对数值。
2)个体预后指数:
Cox比例风险模型的线性参数部分β 1x 1+β 2x 2+…+β px p与风险函数h(t,X)成正比,即,β 1x 1+β 2x 2+…+β px p越大,风险h(t,X)也相应地越大。由此Cox模型的线性参数部分反映了一个个体的预后,β 1x 1+β 2x 2+…+β px p被称为预后指数(prognostic index,PI)。如果预后指数越大,那么,某特定患者的风险也就越大,预后越差;反之,预后指数越小,预后越好。
如果对各协变量进行标准化变换后,得到的Cox模型的线性参数部分,即为标准化的预后指数sPI。当标准化的预后指数为0时,表明某特定患者的风险达到平均风险水平;当标准化的预后指数大于0时,表示该患者的风险超过平均风险水平;当标准化的预后指数小于0时,表示该患者的风险低于平均风险水平。
(6)比例风险假定的检验:
由于Cox比例风险模型具有一系列优点,如,可以同时分析各种协变量对生存时间的影响,并且无须对基准风险分布进行任何假定,特别是大部分统计软件均可容易地实现Cox比例风险模型的拟合与分析等,从而,导致Cox比例风险模型在大量研究中存在着滥用现象。需要特别强调的是,比例风险假定是Cox比例风险模型的最为根本性的假设。这就意味着,各种组合条件下的生存曲线不能交叉。只有当资料满足比例风险假定时,Cox比例风险模型才是统计有效的。因此,对于Cox比例风险模型而言,PH假设检验至关重要。但是,风险函数是基于个体的,而非群体的,这就导致检验PH假定几乎是不可能的。从实用的角度来说,必须首先假定风险函数是群体相似的。在此种假定下,PH假定的近似检验方法主要有两大类:一类是针对每个协变量,单独进行PH假设检验,另一类是针对Cox比例风险模型,进行残差分析,其中,残差主要有鞅残差,偏差残差,Schoenfeld残差,Score残差等。
第一类PH假设检验方法主要有以下三种方法:一是,如若协变量为分类变量时,可以对每一类别分别作Kaplan-Meier生存曲线图,观察各生存曲线间是否有交叉,若无交叉,可以认为满足比例风险假定;二是,以生存时间t为横轴,对数对数生存率ln(- ln (S(t))为纵轴,绘制分类协变量各个类别的生存曲线,如果这些生存曲线平行,可以认为满足比例风险假定;三是,对于连续型协变量,可以将其与生存时间的对数构建交互作用项x jln(t),纳入到Cox比例风险模型之中,如果该交互作用项无统计学意义,那么,可以认为满足比例风险假定。
第二类PH假设检验方法主要有以下两种方法:一是,以残差为纵坐标,时间为横坐标作残差图,从图中判断残差是否存在着某种非随机性模式,如若存在的,则认为不满足PH假定;二是,直接建立残差关于时间的回归模型,如若回归系数有统计学意义,则可以认为不满足PH假定。
需要注意的是,第一类PH假设检验方法,是将各协变量进行分割处理,即使在各单独的协变量下,满足比例风险假定,也并不必然意味着在多个协变量存在着复杂关系的情形下,仍然能够满足比例风险假定。特别是,当各分层的观察对象个数较少时,更是难以判断比例风险假定是否真正成立。第一类PH假设检验方法的统计效能较低,目前更多地使用第二类PH假设检验方法。虽然第二类PH假设检验方法可以从Cox模型的角度进行整体考虑,但是,Cox模型残差本身包含有过多干扰,其性质复杂,难以研究,并且,存在截尾,某些重要协变量未包含在模型中,这些因素均使得第二类方法也难以对PH假定作出较严格的检验。
当比例风险假定不能满足时,主要采用以下三类方法予以处理。一是,将不成比例关系的协变量作为分层变量,然后,再利用其他协变量构建Cox比例风险模型进行分析;二是,采用参数模型替代Cox比例风险模型;三是,采用无须比例风险假定的其他半参数模型,如,半参数加速失效时间模型,半参数转换模型等,对资料予以分析。需要说明的是,第一类替代方法是对于各协变量存在复杂关系的情形并不适用,尤其是,各分层的观察对象个数较少时,所建立的Cox比例风险模型并不可靠;第二类替代方法,则需要对生存时间进行事先分布假定,丧失了半参数模型的优势;第三类替代方法,既无须进行分布假定,也无须比例风险假定,是较好的Cox比例风险模型的替代方案,特别是,半参数加速失效时间模型,得到研究者与实用者的越来越广泛关注。
例11-16为了解影响乳腺癌患者术后生存状态的因素,对32名手术后的乳腺癌患者进行了随访。其收集的影响因素包括远期生活质量评分(KPS),从确诊到入组的时间(Duration),年龄(Age,岁),是否家庭护理(Nurse,否=0,是=1),治疗方法(Therapy,化学治疗=0,生物制剂治疗=1),结果变量为术后生存时间(Time,月)以及随访结局(Status,截尾=0,未截尾=1),数据列于表11-29试对此资料予以分析。
表11-29 32名乳腺癌患者手术后生存资料
续表
在本例中,采用SAS统计软件的PHREG过程拟合Cox比例风险模型,由于在大量研究中,使用的统计方法为Cox比例风险模型,因而,为了保证研究之间的可比性,本例先使用Cox比例风险模型予以分析,再根据残差分析结果,利用半参数加速失效时间模型对本资料进行分析(见下节)。
第一步,对截尾进行分析。由于本实例中,截尾仅有2例,占全部观察个体的6.25%。截尾比例较小,可以直接作生存分析。
第二步,建立Cox比例风险模型。对每个影响因素单独进行单因素Cox回归分析,得到表11-30。由下表可见,在α= 0.05水准上,有统计学意义的因素仅为远期生活质量评分。
表11-30 32名乳腺癌患者术后生存资料单因素Cox回归分析结果
因此,无须对该资料进行多因素Cox回归分析,仅将KPS纳入Cox回归模型之中。由此,可以认为对乳腺癌患者死亡风险有影响的因素是远期生活质量评分。从回归系数的符号和相对危险度的大小来看,该因素是保护性因素。乳腺癌患者KPS每增加10分,术后死亡风险将下降至0.8088倍,即减少19.12%。此研究表明远期生活质量评分越高,其预后更佳。
本例Cox比例风险模型表达式为:
h(t | X)= h0(t)e-0.02122KPS
表达式右侧指数部分取值越大,则风险函数越大,预后相对越差。本例预后指数 PI = -0.02122KPS。可按适当的预后指数分位将观察个体分成低危组、中危组、高危组。对各组制定更为合理的个体化治疗与康复方案,正确指导乳腺癌患者的治疗,以降低其长期死亡风险。
第三步,Cox回归残差分析。相对于其他残差图而言,Schoenfeld残差图的效能更高。理论上,Schoenfeld残差的期望为0,且近似不相关。PH假定下,Schoenfeld残差散点图应围绕0,呈随机波动。如若残差图存在某种趋势或模式,则有理由认为违背PH假定。本例的残差图见图11-14。
图11-14 32名乳腺癌患者Cox模型的Schoenfeld残差图
由图11-14可见,当生存时间大于100时,残差分布不对称,大于0的残差似乎偏多,且呈上升趋势。这表明本例可能违背PH假定,拟合Cox比例风险模型似乎并不太合适,需要拟合其他模型予以分析。进一步分析见下面5.4节。
4.半参数加速失效时间模型
在大量医学研究与应用中,发现尽管Cox比例风险模型具有众多优点,其适用范围相当广泛,但是,仍有众多资科不适合Cox模型来分析。这主要是因为Cox比例风险回归模型事先要求资料满足两个基本假定:一是比例风险假定;二是对数线性假定。特别是,当资料违背比例风险假定时,如果强行拟合Cox比例风险模型,就极有可能得出与实际问题不相符合的解释与预测等,甚至相反的结论。在肿瘤研究中,存在着大量误用、滥用Cox模型的现象,需要可以替代Cox模型的统计方法,更好地分析医学研究数据。
作为Cox比例风险模型的一种很好的替代模型,半参数加速失效时间模型(semi-parametric accelerated failure time model,semi-parametric AFT model)也是一种线性回归模型,它把生存时间的对数作为反应变量,而且误差项的分布也是未知的。该模型首先是由Pieruschka在1961年提出,并应用于加速寿命试验。与Cox模型相比较,加速失效对间模型研究协变量与对数生存时间的回归关系,模型形式更接近于一般的线性回归方程,回归系数的解释也与一般线性回归更相似,模型结果的解释更为简单、直观,同时,也易于理解,有利于被医学研究人员所接受。特别是,当所研究的因素仅仅是延迟或加快事件的起始时间,而不是对整个生存过程产生影响时,加速失效时间模型具有更好的统计效能,而在此情形下,Cox比例风险模型的统计效能极差。
(1)半参数加速失效时间模型的基本形式:
相关公式符号意义见前一节,加速失效时间模型的基本形式如下:
其中,T 0是不考虑协变量时的基准分布的生存时间,T为在协变量条件下的生存时间。在加速失效时间模型中,将生存时间之比的自然对数,协变量对其存在着线性影响,也就是说,生存时间的延长或缩短仅仅与协变量有关,与时间无关。为此,通常将上式取自然对数,将其线性化为:
进一步变形,可化为:
ln(T)=β1x1+β2x2+…+βpxp+ ln(T0)=β1x1+β2x2+…+βpxp+ε
(11-73)
其中,误差项ε表示独立同分布的随机变量。需要注意是,误差项ε是基准生存函数的对数,其均数不一定为0,如果指定误差项的分布,那么,上述模型即为参数加速失效时间模型,否则,就属于半参数加速失效时间模型。从半参数加速时间模型的形式上看,与一般的线性回归模型极为相似,无须对误差项的分布进行指定,这使得该模型具有良好的可解释性与可适用性。
(2)模型的基本假定:
根据半参数加速失效时间模型基本形式,该模型要求资料事先满足两个假定:一是对数线性假定;二是时间尺度比例假定。
1)对数线性假定:
半参数加速失效时间模型同样假定协变量对失效时间的对数的影响为线性模型,时间之比的
的自然对数与各影响因素呈线性关系,服从线性模型的一般规则。各风险因素对时间之比的影响具有乘积性,而不是可加性。
2)时间尺度比例假定:
假设协变量x j,j = 1,2,…,p,均为0时,基准生存函数S 0(t)= e ε0,即,基准生存函数与误差项ε 0存在着指数关系。如果个体具有协变量x,那么,
其中,参数φ= e -(β1x1+β2x2+…+βpxp),称为加速因子(acceleration factor)。根据上式,协变量X的作用只是改变原来生存时间的尺度,改变的幅度大小由φ= e -(β1x1+β2x2+…+βpxp)来决定。当φ<1时,生存时间以一恒定比例被拉长,反之,生存时间以一恒定的比例缩短。其图11-15如下:
图11-15 不同加速因子的生存概率曲线
由于当φ>1时,其生存曲线存在着一个加速度,下降更快,故而,将该模型称为加速失效时间模型。
(3)与风险比的关系:
根据生存函数,可以得到如下风险函数:
h(t,X)= h0(te-(β1x1+β2x2+…+βpxp))e-(β1x1+β2x2+…+βpxp)
(11-75)
假设有两个个体,其协变量的值分别为X和X *,其风险比为:
根据上式,不难推论出,当生存时间分布服从Weibull分布时,加速失效时间模型与Cox比例风险模型是完全等价。另外,由上式也可知,加速失效时间模型无须满足比例风险假定,可以处理更广泛的生存数据类型。
在回归系数解释上,加速失效时间模型回归系数是针对生存时间之比而言,而Cox比例风险模型是针对风险之比而言的。
(4)参数估计:
本节主要介绍Jin所提出的基于秩的加速失效时间模型估计方法。令
,定义残差:
定义两个计数过程N i(β;t)与Y i(β;t),
Ni(β;t)= Iti≤tI{ ei(β)≤t}
Yi(β;t)= I{ ei(β)≥t}
令
那么,回归系数β的加权log-rank估计函数为
或者,
其中,X(b;t)= S (1)(b;t)/S (0)(b;t),w是加权函数。显然,当φ=1时,该统计量就是log-rank检验统计量;当φ= S (0)时,该统计量就是Wilcoxon(Gehan)检验统计量。因此,可以说,半参数加速失效时间模型是乘极限法与寿命表法的一种自然扩展。然而,Cox比例风险模型却难以视为对一个分组因素时非参数生存分析方法的自然扩展。
加速失效时间模型回归系数的估计为加权log-rank估计函数U
w(β)的根。虽然回归系数的估计较易得到,但是,其方差协方差阵极难计算,使得各回归系数的检验统计量与可信区间的计算难以实现。一种可行的方法就是使用bootstrap法,重复抽样得到回归系数估计的经验分布,从而,得到回归系数
的可信区间及其检验统计量P值。
(5)模型参数的意义及其解释 1)回归系数:
与Cox比例风险模型回归系数的流行病学意义不同的是,加速失效时间模型回归系数具有比较明确的临床意义,它直接反映了各协变量对各个个体生存时间的影响。由加速失效时间模型的公式,可以得到
S(t|X)= S0(te-(β1x1+β2x2+…+βpxp))
(11-80)
回归系数β j的指数e βj称为时间比(time ratio)。其不会随着时间的变化而改变,保持恒定。它与生存函数S(t│X)之间有如下关系:当e βj>1时,则随着协变量x j的增大,生存函数S(t│X)的递减速度也相应减小,表示患者的生存时间越长;当e βj<1时,则随着协变量x j绝对值的增大,生存函数S(t│X)的递减速度也相应增大,表示患者的生存时间越短;当e βj=1时,则协变量x j对生存函数S(t│X)没有影响,表示患者的生存时间无变化。
从上面的公式,难以直观地看出各协变量的意义,采用分位函数(quantile function)的形式,可将其变换如下:
Q(p|X)= Q0(p)e-(β1x1+β2x2+…+βpxp))
(11-81)
由上式可见,对某个协变量特定取值而言,在保持其他协变量取值不变的情形下,回归系数β j可以解释为该协变量变化一个单位,其生存时间变化e βj倍(图11-16)。
图11-16 加速失效时间模型
需要说明的是,根据上图中的左图,两组的生存时间具有明显差异,而通过右图,两组的风险基本一致,这是因为两条生存函数曲线的斜率基本相等之故。由此可见,相对危险度的高低与生存时间的长短并不存在必然联系,相对危险度高不意味着生存时间短,相对危险度为1,也不意味着生存时间相等,相对危险度低也不意味着生存时间就长。与相对危险度不同的是,时间比的解释比较直观,直接反映生存时间的长短,时间比是临床研究的关注焦点。而相对危险度,是流行病学研究的关注焦点,它只能间接地反映个体的生存时间。
2)个体加速因子:
如果加速因子越大,那么,某特定患者的整个生存期越短;反之,加速因子越小,整个生存期越长。当φ>1时,其生存曲线下降较快,整个生存期较短;当φ<1时,其生存曲线下降较慢,整个生存期较长;当φ=1时,其生存曲线与基准生存曲线相同,整个生存期与基准生存期相同。
如果对各协变量进行标准化变换后,得到的加速失效时间模型的线性参数部分,即为标准化的时间比。当标准化的时间比为1时,表明某特定患者的整个生存期达到平均生存水平;当标准化的时间比大于1时,表示该患者的整个生存期短于平均生存水平;当标准化的时间比小于1时,表示该患者的整个生存期长于平均生存水平。
3)时间尺度变化比例:
加速失效时间模型可以描述任意两个个体生存状态之间的关系。对于具有协变量X和X *的两个个体(前者为甲个体,后者为乙个体),其生存函数之间的关系为:
上面公式表明,任意两个个体仅仅是生存时间尺度发生变化,时间尺度变化比例为δ *,表示二者之间的整个生存期相差δ *倍。当时间尺度比例δ *>1时,表示甲个体整个生存期要短于乙个体;当时间尺度比例δ *<1时,表示甲个体整个生存期要长于乙个体;当时间尺度比例δ *=1时,表示甲个体与乙个体的整个生存期相当。例如,若δ *=4,则表示这两个个体在整个生存期均相差4倍。
(6)模型选择:
从理论上讲,Cox比例风险模型的两个基本假定——对数线性假定与比例风险假定难以验证,加速失效时间模型的两个基本假定——对数线性假定与时间尺度比例假定也同样难以验证。因此,从实用的角度来说,选择何种统计模型的较佳途径是依据数据的拟合程度和可解释性来予以判断。
如果加速失效时间模型对数据的拟合程度和可解释性要明显好于Cox比例风险模型,那么,应采用加速失效时间模型。反之亦然。但是,如果这两模型对数据的拟合程度与可解释性相近,那么,就必须要考虑其他的一些因素。例如,对于同一个影响因素采用相对危险度来衡量其对生存时间的影响,并且,其他研究者均采用Cox比例风险模型。此时,为了便于各个研究之间的比较与分析,最好还是采用Cox比例风险模型。反之,当研究目的是为了得到更好的预测效果,并且,需要各影响因素对生存时间影响更为直观的解释时,则采用加速失效时间模型无疑要优于Cox比例风险模型。
到目前为止,半参数加速失效时间模型在参数估计、统计检验与实际效能验证等方面尚不如Cox比例风险模型研究的成熟与全面,并且,缺乏可靠的、易用的统计软件来实现半参数加速失效时间模型,加上,现在的医学研究人员对该模型缺乏深入的了解,这些都导致半参数加速失效时间模型没有得到广泛的应用,较少见于各种医学研究文献之中。
例11-17对例11-16的资料配合半参数加速失效时间模型。资料见表11-29。对该资料配合Cox比例风险模型后的残差分布不对称,大于0的残差似乎偏多,且呈上升趋势。这表明本例可能违背PH假定,利用半参数加速失效时间模型对本资料进行分析。对每个影响因素单独进行单因素半参数AFT回归分析,采用R统计软件中的rankreg程序包拟合半参数加速失效时间模型。得到表11-31。由于半参数加速失效时间模型的参数估计与加权函数的选择有着直接的关系,因此,为了间接排除权重函数对其参数估计的影响,将采用两种常用的加权函数,即,Gehan加权与Logrank加权,对半参数加速失效时间模型的回归系数同时进行估计。如果在两种不同加权函数的情形下,某影响因素的回归系数的估计值相近,可以说明该影响因素对乳腺癌患者存在着较为确切的影响;如若不然,则需要谨慎下结论。
由表11-31可见,无论在Gehan加权的情形下,还是在Logrank加权的情形下,在α=0.05水准上,有统计学意义的因素同样为远期生活质量评分。
表11-31 32名乳腺癌患者术后生存资料单因素半参数AFT模型分析结果
因此,也无须对该资料进行多因素半参数AFT模型分析,仅将KPS纳入到半参数
AFT模型之中。根据上表结果,可以认为对乳腺癌患者生存时间有影响的因素是远期生活质量评分。从回归系数的符号和相对危险度的大小来看,该因素是保护性因素。乳腺癌患者KPS每增加10分,术后整个生存期将延长至1.40倍左右,延长约40%。此研究表明远期生活质量评分越高,其生存时间越长。
本例半参数AFT模型表达式为:
S(t│X)= S0(te-0.011KPS)
或
S(t│X)= S0(te-0.013KPS)
表达式右侧指数部分取值越大,则加速因子越大,生存时间相对越短。本例加速因子为e -0.011KPS或e -0.013KPS。同样,也可按适当的加速因子将观察个体分组,正确指导乳腺癌患者的治疗与康复。
综上所述,从本例资料来看,无论是Cox比例风险模型,还是半参数AFT模型,均表明远期生活质量评分对于乳腺癌患者术后生存起着至关重要的作用,如何有效提高远期生活质量是乳腺癌患者术后治疗与康复的关键性问题。
5.生存分析应用的注意事项
(1)癌症患者的随访期较长,影响因素多,在研究期间,易引入各种混杂因素。因此,不能轻易下结论。生存分析结论的正确性往往不在于方法的选择,更多地在于科学的研究设计与良好的质量控制。如何保证癌症患者生存研究的质量与过程控制才是生存分析的真正关键所在。
(2)在应用各种生存分析回归模型时,必须对各种模型的前提进行检验。例如,Cox比例风险模型的PH假定,半参数AFT模型的时间尺度比例假定等。这些假定都是对客观事物的一种抽象与近似,并不是客观存在的,应更多地从所收集的数据出发,在准确把握数据特性的基础之上,再决定采用何种模型更为合理。这就要求在进行癌症患者生存研究之前,必须要有对研究对象的深入了解与分析和(或)良好的预试验,事先确定各种生存模型的备选路径方案,而不是事后盲目选择。
(3)在建立生存模型的过程中,除了参考统计方面的证据之外,应更多地依赖于专业上的理论与可解释性,来判断研究是否真正满足生存分析的假设前提,删除或增加变量应相当谨慎。另外,为了保证生存分析的严密性,必须进行生存分析模型的敏感性分析,特别是截尾对生存模型的影响,以了解所建立的模型是否具有较强的稳健性。如若模型不太稳健,则需要对整个研究重新进行审视,而不是从统计上作修补的工作。
(陈心广 蒋红卫 余松林)