1、生存分析概述及其应用实例 侯笛 摘 要 : 本文对概括性地介绍了生存分析的概念和主要的研究内容。对生存分析中常用的术语进行了描述,包括 生 存数据,生存函数,风险函数等。并阐述了生存 数据估计中 常用的三种模型(非参数模型,参数模型,半参数模型。)的原理和方法。最后借助 SPSS软件,用不同的生存分析模型对白鼠的治疗实验数据进行了统计处理,最后得出了影响白鼠生存时间的主要因素。 关键词 :生存分析; SPSS;生存数据 生存分析是指将终点事件和出现此事件所经历的时间结合起来分析的一种统计分析方法,研究生存现象和现象的响应时间数据以及其规律,是处理以生存时间( survival time)为反应
2、变量、含有删失数据一类资料的统计方法。此类资料的生存时间变量大多不服从正态分布,且常含有删失值,故不适于用传统的数据分析方法如 t 检验或线性回归进行分析。所谓删失值,就是因各种原因对随访对象的随访可能失访或终检。能处理删失数据也是生存分析的一个优点。作为统计科学的一个重要分支,生存分析的应用已经由最初的医学领域扩展到了社会科学、经济学等诸 多领域。 “事件 ”和 “寿命 ”是生存分析研究中的两个重要变元。生存分析研究的主要内容包括以下两个方面:( 1)描述生存过程;( 2)分析生存过程的影响因素并对生存的结局加以预测。通常将生存分析所用到的数据称为生存数据,生存数据又分为完整数据和删失数据。
3、生存分析中常用的分析方法很多,按照是否使用参数可以分为非参数方法,参数方法和半参数方法 1。下面将从基本概念,生存分析方法,实例解析三方面作具体的介绍。 1.生存分析中的基本概念 1.1 事件和寿命 生存生存分析中定义的事件有死亡、损坏、失败、解雇、病发等等。例如病人的死亡,产品的失效,疾病的发生,职员被解雇。而寿命则是指试验或记录开始到事件发生所经历的时间。 1.2.生存数据 生存数据可以分为完全数据和删失数据。完全数据是指提供了完整信息的数据。例如,在研究产品的失效时间时,某个样品从进入研究直到失效都在我们的观察中,可以得到该样品的具体失效时间,这就是一个完全数据。删失数据则是指在观测期内
4、,我们并没有看见个体的状态发生改变,无法确定个体具体的生存时间。又分为左删失数据,右删失数据,区间删失数据。 生存分析中最常见的是右删失数据,是指试验个体在观测终止前(或者离开实验前)一直存活或工作,其失效或死亡发生在观测终止(或离开实验)之后。左删失数据是指我们知道某个体在某时刻前已经失效或死亡,但不知道具体的时间。区间删失数据是指个体在某个时间区间内死亡,但不能确定具体的死亡时间。 1.3 生存函数和风险函数 生存分析中有两个重要的函数:生存函数和风险函数 2。生存函数( survival function) ttint()iii n dS又称为累积生存率,记作 S( t),是指个体生存时
5、间超过 t 的概率。风险函数( hazard function)又称瞬时死亡率,记作 h( t),是指在 t 时刻存活的个体,在 t+ t 时刻死亡的概率。 用 T 表示个体生存时间, t 表示观测时间,可以用下式描述两个函数: S( t) =P( T t) ( 1.1) ( 1.2) 2.生存分析方法分类 2.1 非参数方法 非参数方法是生存分析中最常用的一种方法。当研究的事件不能被参数模型很好地拟合时,通常可以采用非参数方法研究它的生存特征。常用的非参数方法包括生命表分析和Kaplan-Meier 分析。 生命表法可以用来描述群体的生存现象,适用于大样本的情况,它的主要优点在于对生存 时间
6、的分布没有限制。生命表分析将观测时间分成时间段,按时 间段逐个统计事件发生的情况,以此估计生存函数。假设共有 k个 时间段 t 0 , t 1) , t 1 , t 2) , , t k-1 , t k ) , 每个区间中事件发生的次数分别为 d 1 , d 2 , , d k , 每个区间中的个体总数分别为 n 1 , n 2 , , n k ,所以在第 i 个区间个体存活的概率为( n i - d i ) / n i ,而个体可以从第一个区间存活到第 i 个区间的概率 (累积生存率)为: (2.1) 易知 S( t)为递减函数。 3 Kaplan-Meier 分析,也 称 为乘积极限估计,
7、是 Kaplan 和 Meier 二人在 1958 年提出的一种估计生存函数的非参数方法。与生命表分析不同, K-M 分析以事件发生的时间点将观测区间分段,用来估计生存函数。 2.2 参数方法 若已经证明某事件的发展可以用某个参数模型很好地拟合,就可以用参数方法做该事件的生存分析。常用的参数模型有指数分布模型、对数分布模型、正态分布模型,威泊分布模型, Gompertz 分布模型等 2。 指数分布是一种应用极广泛的单参数分布,它最大的特点是风险函数为常数,即: h( t) = (2.2) Gompertz 分布是 Gompertz 在 1825 年提出的,常用于刻画人的生存分布,其风险函数为:
8、 t ttTtPlim =)t( 0t )(h( t) =exp( t), t 0 ( 2.3) Gompertz 分布能够比较好地拟合出生物出生、成长、衰老死亡的全过程,所以目前在人口,保险精算,生物医疗等领域都得到了广泛应用。 威泊分布对应的风险函数为: H(t) = pt p-1, t 0 ( 2.4) 当 p 1 时,风险函数是时间的增函数; p=1 时,风险函数是常数,威泊分布退化为指数分布; p 1 时,风险函数是时间的减函数。 2.3 半参数方法 生存分析中我们常常遇到个体的生存状况受到多种因素影响的 情况。这些对生存时间有影响的变量称为协变量。在分析生存数据时要将协变量的影响考
9、虑进去。 Cox 半参数模型就很好地解决了这个问题。它假定风险函数由两部分构成:基准风险函数和协变量线性组合的指数。 Cox 半参数模型又分为独立协变量比例风险模型和时间相依性协变量比例风险模型两种。二者的区别在于协变量的取值是否和时间有关。 Cox 独立协变量比例风险模型可以写成如下形式 2: 122) (2.5) 式中, Z1, Z2, Z m 为协变量,这里的协变量与时间无关, 1, 2, m 为对应协变量的未知参数。 h 0( t)是基准风险函数,是当所有协变量的值为 0 时在 t 时刻风险函数的值。 实际应用中常常计算两个不同个体风险函数的比率,称为危险率。假定给定的两个个体的协变量
10、观察值分别为( Z1, Z2, Z m)和( Z1*, Z2*, Z m*),那么对应的危险率为: miiimmmmmmzzzzzzzzzzztzzzt)(e x p )e x p ()e x p (),(),(h),(),(h221122112121( 2.6) 因此危险率是与时间无关的常数,式 (2.5)因此称为风险比例模型。 实际应用中协变量常常是随时间变化的,此时危险率不再是常数,协变量随时间变化的 Cox 模型称为时间相依性协变量比例风险模型。 3.实例分析 生存分析常用的软件有 SAS 和 SPSS 两种统计分析软件,本文采用 SPSS进行统计分析。 两组小白鼠用来检验癌症的治疗状
11、况。一组使用传统治疗方法,另一组使用试验方法,试验人员记录了小白鼠的存活时间及状态: Days 为存活时间或观测时间(天); Status 表示生存状态,取值 1 表示死亡, 0 表示存活; Group 表示治疗方法,取 0 表示传统疗法,取 1)exp()(),(),(h 2211021 mmm zzzthzzzt 表示试验疗法,共有 64 组数据。 表 3. 1 不同治疗方法下白鼠的存活时间 Days Status Group Days Status Group 144 1 0 196 0 1 164 1 0 198 1 1 188 1 0 204 0 1 188 1 0 205 1 1
12、190 1 0 208 0 1 192 1 0 220 0 1 206 1 0 240 1 1 209 1 0 242 1 1 210 0 0 242 1 1 214 1 0 244 1 1 216 0 0 244 1 1 216 1 0 244 1 1 220 1 0 244 1 1 220 1 0 249 1 1 227 1 0 254 0 1 227 1 0 261 1 1 240 1 0 266 0 1 244 0 0 266 0 1 244 0 0 269 1 1 244 0 0 280 1 1 244 1 0 280 1 1 246 1 0 284 0 1 246 1 0 296
13、1 1 256 0 0 296 1 1 265 1 0 330 0 1 265 1 0 340 1 1 266 1 0 344 0 1 304 1 0 348 0 1 344 1 0 142 1 1 404 1 0 414 0 1 156 1 1 142 1 1 164 1 1 424 1 1 对于一批生存数据 , 在事先不知道寿命分布的总体趋势 , 且又不好判断应该用何种模型 最合适时 , 多数学者一般直接采用非参数方法或半参数法。但是 , 由于非参数方法的精度一 般低于参数方法。常用的非参数法为生命表分析法和 K-M方法,生命表方法对生存函数和生存率的估计依赖于生命表中所有的区间,应用假定
14、是总体在每个区间内各处有近似相等的生存概率,乘积限方法是基于一个个数据的估计,而非生命表中按区间分组的数据,其假设为事件发生的概率仅依赖于时间 4。本例中并不知道个体在各个区间的生存概率是否相等,依次用生命表法和 K-M法进行分析。 3.1 生命表法 图 3.1 生命表法得到的生存函数 图 3.2 生命表法得到的风险函数 观察生存函数的分布可以看出,大约在 200 天时两种治疗方法的生存函数相交,在 200天以前传统治疗方法的存活率较高,而在 200 天以后试验方法的治疗效果明显优于传统治疗方法。 观察风险函数的分布,试验方法的风险函数基本分布在传统方法一下,风险较低。试验方法的分布趋势整体呈
15、现较为稳定的低水平,而传统方法的风险则整体呈现增长趋势。 表 2 的显著性检验结果 p 值为 0.088 0.05,说明两种治疗方法并不存在显著性差异。 3.2 K-M分析 表 3.3 K-M分析 生存表 治疗方法 时间 状态 此时生存的累积比例 累积事件数 剩余个案数 估计 标准误 0 1 144.000 1 .967 .033 1 29 2 164.000 1 .933 .046 2 28 3 188.000 1 . . 3 27 4 188.000 1 .867 .062 4 26 5 190.000 1 .833 .068 5 25 6 192.000 1 .800 .073 6 24
16、 7 206.000 1 .767 .077 7 23 8 209.000 1 .733 .081 8 22 9 210.000 0 . . 8 21 10 214.000 1 .698 .084 9 20 11 216.000 1 .663 .087 10 19 12 216.000 0 . . 10 18 13 220.000 1 . . 11 17 14 220.000 1 .590 .092 12 16 15 227.000 1 . . 13 15 16 227.000 1 .516 .094 14 14 17 240.000 1 .479 .094 15 13 18 244.000
17、1 .442 .094 16 12 19 244.000 0 . . 16 11 表 3.2 检验结果 Wilcoxon (Gehan) 统计量 df Sig. 2.920 1 .088 20 244.000 0 . . 16 10 21 244.000 0 . . 16 9 22 246.000 1 . . 17 8 23 246.000 1 .344 .095 18 7 24 256.000 0 . . 18 6 25 265.000 1 . . 19 5 26 265.000 1 .229 .092 20 4 27 266.000 1 .172 .085 21 3 28 304.000
18、1 .115 .073 22 2 29 344.000 1 .057 .055 23 1 30 404.000 1 .000 .000 24 0 1 1 142.000 1 . . 1 33 2 142.000 1 .941 .040 2 32 生存表 治疗方法 时间 状态 此时生存的累积比例 累积事件数 剩余个案数 估计 标准误 1 3 156.000 1 .912 .049 3 31 4 164.000 1 .882 .055 4 30 5 196.000 0 . . 4 29 6 198.000 1 .852 .061 5 28 7 204.000 0 . . 5 27 8 205.00
19、0 1 .820 .067 6 26 9 208.000 0 . . 6 25 10 220.000 0 . . 6 24 11 240.000 1 .786 .072 7 23 12 242.000 1 . . 8 22 13 242.000 1 .718 .080 9 21 14 244.000 1 . . 10 20 15 244.000 1 . . 11 19 16 244.000 1 . . 12 18 17 244.000 1 .581 .090 13 17 18 249.000 1 .547 .091 14 16 19 254.000 0 . . 14 15 20 261.000
20、 1 .510 .092 15 14 21 266.000 0 . . 15 13 22 266.000 0 . . 15 12 23 269.000 1 .468 .093 16 11 24 280.000 1 . . 17 10 25 280.000 1 .383 .094 18 9 26 284.000 0 . . 18 8 27 296.000 1 . . 19 7 28 296.000 1 .287 .092 20 6 29 330.000 0 . . 20 5 30 340.000 1 .230 .089 21 4 31 344.000 0 . . 21 3 32 348.000
21、0 . . 21 2 33 414.000 0 . . 21 1 34 424.000 1 .000 .000 22 0 K-M分析方法与生命表法相比记录了删失数据,由于生命表的结果不够直观,最好观察生存函数分布,可以较为清楚的看到试验疗法与传统疗法的差异。 图 3.3 K-M 分析 累积 生存函数 表 3.4个案处理摘要 治疗方法 总数 事件数 删失 N 百分比 0 30 24 6 20.0% 1 34 22 12 35.3% 整体 64 46 18 28.1% 图 3.4 K-M 分析 累积 风险函数 K-M分析的生存函数分布与生命表法的生存函数类似, 200 天左右是两种治疗方法生存率的
22、分界线。但由于是按照事件的发生分段,区段较多,整体呈现密集的锯齿,而生命表分析的分布则较为平缓。 观察 累积 风险函数的分布, 发现 传统治疗方法的 累积风险率高于试验方法 , 这种趋势在观测后期更加明显。 表 3.5不同检验方法的 比较 卡方 df Sig. Log Rank (Mantel-Cox) 4.795 1 .029 Breslow(Generalized Wilcoxon) 2.967 1 .085 Tarone-Ware 3.762 1 .052 由表 3.5 知,不同检验方法呈现不同的结果,其中 Log Rank 检验的 p 值小于 0.05,该检验表明两种治疗方法有显著性差
23、异。 除了治疗方法对小白鼠的生存状况有影响,其他因素如性别,年龄,体重等都可能对其生存时间造成影响。考虑到这些协变量的影响,结合相关数据,用 Cox 独立协变量比例风险模型重新分析。其中年龄的单位是天,体重的单位是盎司, F 代表雌性, M代表雄性。 Days Status Group Sex Age Weight 144 1 0 F 94 10 164 1 0 M 50 10 188 1 0 M 57 11 188 1 0 F 81 11 190 1 0 F 89 10 192 1 0 M 71 11 206 1 0 M 69 11 209 1 0 F 81 12 210 0 0 F 99
24、13 214 1 0 M 52 12 216 0 0 M 69 12 216 1 0 M 67 13 220 1 0 M 51 13 220 1 0 F 78 12 227 1 0 M 58 12 227 1 0 F 79 13 240 1 0 F 92 13 244 0 0 F 98 13 244 0 0 F 79 12 244 0 0 M 68 13 244 1 0 M 59 13 246 1 0 F 88 13 246 1 0 M 73 13 256 0 0 F 83 12 265 1 0 F 79 12 265 1 0 M 66 12 266 1 0 M 59 13 304 1 0 M 67 13 344 1 0 M 69 13 404 1 0 F 86 12 156 1 1 F 97 11 164 1 1 M 72 10 196 0 1 M 61 10 198 1 1 M 60 11 204 0 1 M 62 11 205 1 1 F 81 10 208 0 1 M 54 12 220 0 1 F 78 12