十八章生存分和COX回归

1、第十八章第十八章生存分析和生存分析和coxcox回归回归上海第二医科大学生物统计教研室第一节第一节 基本概念基本概念 在医学,生物学研究中,常用到生存分析 (survival analysis)方法。例如对于肿瘤等疾病的疗效及预后的考核,通常不用治愈率，有效率等表示，而用将来复发或死亡的时间长短表示,也即生存期来表示。 所谓生存期（survival time）是指从某个标准时刻(如发病,确诊,开始治疗或进行手术的时间)算起至死亡或复发为止的时间。 生存期不同于一般指标的二个特点：1.有截尾数据（censored data) 随访中未能知道病人的确切生存时间，只知道病人的生存时间大于某时间。（1

2、）病人失访或因其他原因而死亡-失访（2）到了研究的终止期病人尚未死亡-终访截尾数据可记为t+,如： 4+ = 生存时间大于4年。 虽然截尾数据提供的信息是不完全的，但不能删去，因为这不仅损失了资料，而且会造成偏性。2. 生存期的资料一般不服从正态分布。 由于上述原因，常用的统计方法不适用，而要用特殊的统计方法。 生存分析是指对于生存期这一指标进行分析的一系列特殊的统计方法。 生存时间不一定专用于死与活的情况,生存时间(存活时间)可定义为从某种起始事件到达某终点事件所经历的时间跨度。例如急性白血病病人从治疗开始到复发为止之间的缓解期;冠心病病人在两次发作之间的时间间隔;已作输卵管结扎的妇女从施行

3、输卵管吻合手术后至受孕的时间间隔;在流行病学研究中,从开始接触危险因素到发病所经历的时间等都可作为生存时间用作生存分析。 有时还收集一些有关因素(称为自变量或协变量), 以分析这些协变量是否对生存时间有影响,影响的大小，是缩短或延长生存时间。这可以通过cox回归进行分析,因此,cox 回归可看成带有协变量的生存分析。 包括： (1)开始观察日期，终止观察日期-生存时间 (2)结局（最终的观察到的是死亡还是存活） 死于该病-完全数据 存活或死于其他原因-截尾数据 每个生存期数据要用2个变量表示：观察到的生存时间和是否截尾（如：用1表示截尾，用0表示死亡；4+ 用4，1表示；4用4，0表示）。 (

4、3)协变量-各种影响生存期长短的因素。随访资料的记录：第二节 描述生存时间分布规律的函数 一. 生存率(survival rate) 又称为生存概率或生存函数,它表示一个病人的生存时间长于时间t的概率,用s(t) 表示： s(t)=p（tt） 如5年生存率： s(5)=p（t5） 以时间t为横坐标,s(t)为纵坐标所作的曲线称为生存率曲线, 它是一条下降的曲线,下降的坡度越陡,表示生存率越低或生存时间越短,其斜率表示死亡速率。 1.2 概率密度函数 (probability density function) 简称为密度函数,记为f(t),其定义为: f(t)=lim (一个病人在区间(t,t

5、+t)内死亡概率/t) 它表示死亡速率的大小。如以t为横坐,f(t) 为纵坐标作出的曲线称为密度曲线,由曲线上可看出不同时间的死亡速率及死亡高峰时间。纵坐标越大,其死亡速率越高,如曲线呈现单调下降,则死亡速率越来越小,如呈现峰值,则为死亡高峰。 1.3 风险函数(hazard function) 用h(t)表示,其定义为: h(t)=lim(在时间t生存的病人死于区间(t,t)的概率/t) 由于计算h(t)时,用到了生存到时间t,这一条件,故上式极限式中分子部分是一个条件概率。可将h(t)称为生存到时间t的病人在时间t的瞬时死亡率或条件死亡速率或年龄别死亡速率。当用t作横坐标,h(t)为纵坐标

6、所绘的曲线,如递增,则表示条件死亡速率随时间而增加,如平行于横轴,则表示没有随时间而加速(或减少)死亡的情况。 风险函数的不同情况：常数， 如：死于飞机失事。下降， 如：急性损伤。上升， 如：持续接触危险因素。澡盆样，如：人的一生。 生存分析目的：（1）估计生存函数。（2）比较各组的生存函数。（3）研究影响生存期长短的因素。第三节 生存率的估计方法 生存率s(t)的估计方法有参数法和非参数法。常用非参数法，非参数法主要有二个,即,乘积极限法与寿命表法,前者主要用于观察例数较少而未分组的生存资料,后者适用于观察例数较多而分组的资料,不同的分组寿命表法的计算结果亦会不同,当分组资料中每一个分组区间

7、中最多只有 1个观察值时,寿命表法的计算结果与乘积极限法完全相同。t(年)s(t)0110.6720.4530.340.250.1460.0970.06 参数法可求出一个方程表示生存函数s(t)和时间t的关系，画出的生存曲线是光滑的下降曲线。 非参数法只能得到某几个时间点上的生存函数，再用直线联起来，画出的生存曲线是呈梯型的。一. 乘积极限法(product-limit method) 简称为积限法或pl法,它是由统计学家kaplan和meier于1958年首先提出的, 因此又称为kaplan-meier法, 是利用条件概率及概率的乘法原理计算生存率及其标准误的。 设s(t)表示t年的生存率,

8、s(ti/ti-1)表示活过ti-1年又活过 ti年的条件概率,例如s(1),s(2)分别表示一年,二年的生存率,而s(2/1)表示活过一年者,再活一年的条件概率,据概率的乘法定律有: s(2)=s(1)s(2/1),一般地有 s(ti)=s(ti-1)s(ti/ti-1) 例22.1 用某中药加化疗(中药组)和化疗(对照组)两种疗法治疗白血病后, 随访记录各患者的生存时间,不带+号者表示已死亡,即完全数据,带+ 号者表示尚存活,即截尾数据,试作生存分析。时间单位为月。 中药组 10,2+,12+,13,18,6+,19+,26,9+,8+,6+,43+,9,4,31,24 对照组 2+,13

9、,7+,11+,6,1,11,3,17,7资料中药组积限法计算生存率时间 状态 期初人数 死亡人数 条件生存率 累积生 di di/ni(ni-di)累积生存 ti si ni di (ni-di)/ni 存率s(ti)ni(ni-di) 率标准误 = 2 活 4 死 15 1 0.9333 0.9333 0.004762 0.004762 0.0644 6 活 6 活 8 活 9 死 11 1 0.9090 0.8485 0.009091 0.013853 0.0999 9 活 10 死 9 1 0.8889 0.7542 0.013889 0.027742 0.1256 12 活 13 死

10、 7 1 0.8571 0.6465 0.023810 0.051551 0.1468 18 死 6 1 0.8333 0.5387 0.033333 0.084885 0.1570 19 活 24 死 4 1 0.7500 0.4040 0.083333 0.168218 0.1657 26 死 3 1 0.6667 0.2694 0.166667 0.334885 0.1559 31 死 2 1 0.5000 0.1347 0.500000 0.834885 0.1231 43 活积限法的计算步骤为: 1.将n个生存数据ti,按小到大排列,当截尾数据与完全数据(非截尾值)相同时,截尾数据排

11、列在后,并写出每个生存数据的状态si(即死或活),见表22.1的1,2列 2.写出各个完全数据(即死亡状态)的期初人数ni和死亡人数di,见表22.1的3,4列。 3.计算条件生存率的估计值,见表中第5列,s(ti/ti-1)=(ni-di)/ni 4.计算累积生存率,即时间ti的生存率估计值(见表中第6列) s(ti)=s(ti-1)s(ti/ti-1) 5.计算s(ti)的标准误(见表中第9列) i se(sti)=s(ti)dj/nj(nj-dj) j=1 表中已列出了积限法的全部结果,各个时间点的生存率和标准误分别在6,9两列,例如二年生存率(即24个月)为0.40400.1657二.

12、 寿命表法(life table method) 适用于随访的病例数较多, 将资料按生存期进行分组,在分组的基础上计算生存率 ,本法也能用于不分组的资料,此时计算结果与积限法相同。 某医院1946年1月1日到1951年12月31日收治的126例胃癌病例,生存情况如表22.2,试用寿命表法估计生存率。 表22.2 126例胃癌患者寿命表法估计生存率 时间(年) 期初例数 死亡例数 失访例数 截尾例数 有效例数 条件生存率 累积生存率 di di/ni(ni-di)累积生存 ti ni di ui wi ni s(ti/ti-1) s(ti) ni(ni-di) 率标准误 = 0- 126 47

13、4 15 116.5 0.5966 0.5966 5.80510-3 5.80510-3 0.0455 1- 60 5 6 11 51.5 0.9029 0.5386 2.08810-3 7.89310-3 0.0479 2- 38 2 0 15 30.5 0.9344 0.5033 2.30110-3 0.0102 0.0508 3- 21 2 2 7 16.5 0.8788 0.4423 8.35910-3 0.0186 0.0602 4- 10 0 0 6 7.0 1.0000 0.4423 0 0.0186 0.0602 5- 4 0 0 4 2.0 1.0000 0.4423 0 0

14、.0186 0.0602 寿命表法估计生存率步骤如下: 1.将观察例数按时间段(年)0-,1-,2-,划分,分别计数期初例数,死亡,失访, 截尾例数列入表22.2的1-5列。事实上,从第二个时间段开始,期初人数ni 系由下式算得: ni=ni-1-di-ui-wi 例如第二行,即时间段1-,有 n2=126-47-4-15=60 2.计算各时间段期初实际观察例数,(亦称有效例数)ni ni=ni-ui/2-wi/2 上式表明该时间段期初例数中的失访,及截尾例数只计其半时,即得有效例数。 如第一行,n1=126-4/2-15/2=116.5 3.分别用(22.5)(22.6)(22.7)式计算条

15、件生存率s(ti/ti-1),累积生存率s(ti)及其标准误。 计算结果已列于表22.2中,第7,8,11列,表中9,10二列系用于第11列的计算。 例如时间段0-中 s(ti/ti-1)=(116.5-47)/116.5=0.5966 s(ti)=10.5966=0.5966 se(s(ti)=0.59665.80510-3=0.0455 故一年生存率的估计为0.59660.0455 同样二年生存率的估计为0.53860.0479 由于寿命表法与积限法的累积生存率及其标准误的计算公式完全相同,所以,当分组资料中每一个分组区间中最多只有1个观察值时,寿命表法就是积限法。第四节 生存率的比较 当

16、有两个或两个以上的生存分布时,我们常需比较它们是否来自同一生存分布,此时的假设检验为: h0:样本所来自的总体生存分布相同。 h1:样本所来自的总体生存分布不相同。 可选用的检验方法有:logrank法, 广义wilcoxon法,和cox-mantel法等。当拒绝h0时,认为几个生存分布不相同。 当不需要整体比较，而只要比较个别时间点上几组生存率时可用下面方法：（1）两个生存率比较 生存率s1和s2，其方差为v1和v2 用卡方检验： 2 =（s1-s2）2 / （v1+v2） df=1（2）两个以上两个生存率比较 生存率s1，s2和s3，方差为v1，v2和v3 用卡方检验： 权重w1=1/v1

17、, w2=1/v2，w3=1/v3 加权平均生存率： s=（w1*s1+w2*s2+w3*s3） /（w1+w2+w3） 2 =w1*(s1-s)2 +w2*(s2-s)2 +w3*(s3-s)2 df=3-1表22.3 例22.1资料两疗法生存比较的logrank检验时间 秩次 中药组序号 死亡例序号 处于危险 logrank (非截尾数据) 状态例数 变换值wti(年) i1 i2 i3 r 1 1 1 26 -0.96 2+ 2 1 0.04 2+ 3 0.04 3 4 2 23 -0.92 4 5 2 3 22 -0.87 6 6 4 21 -0.82 6+ 7 3 0.18 6+ 8

18、 4 0.18 7 9 5 18 -0.77 7+ 10 0.23 8+ 11 5 0.23 9 12 6 6 15 -0.70 9+ 13 7 0.30 10 14 8 7 13 -0.63 11 15 8 12 -0.54 11+ 16 0.46 12+ 17 9 0.46 13 18 10 9 9 -0.37 13 19 10 8 -0.37 17 20 11 7 -0.16 18 21 11 12 6 0 19+ 22 12 1.00 24 23 13 13 4 0.25 26 24 14 14 3 0.59 31 25 15 15 2 1.09 43+ 26 16 2.093.1 l

19、ogrank检验(log rank test) 当比较的几个样本生存分布,全部为完全数据时,本检验又称为savage检验。 logrank检验的计算步骤如下:1.将两样本的生存数据混合,由小到大排列,并给以秩次i1, 当截尾数据与完全数据数值相同时,截尾数据排列在后。并设两样本含量分别为m1,m2,总例数n=m1+m2。 例22.1中药组与对照组生存数据排列结果见表22.3中第1,2列。2.列出所比较的两组中任一个组的序号i2(本处选用中药组),记入表22.3中第3列。3.列出死亡例的序号i3(见表22.3中第4列)。4.计算非截尾数据(完全数据)各时间点处于危险状态的例数r,它表示该时刻时还

20、剩下多少例数。r系由与i3相应的i1值计算而得。 r=n-i1+1 (22.10)例如与生存期7(月)相应的r值系由 r=26-9+1=18 算得,见表中第5列。5.对秩次i1作logrank变换,即计算logrank变换值w,其算法为秩次为i1 的序号为i3非截尾数据的w值为 i3 w= 1/rj-1 (22.11) j=1 秩次为i1的截尾数据,首先判断它在那二个非截尾数据之间,如果它在序号i3与i3+1之间则w为 i3 w= 1/rj (22.12) j=1特别地,当截尾数据在第一个非截尾数据之前时,取w=0,几个截尾数据落在同样序号的非截尾数据之间时,它们具有相同的w值。例如表22.3

21、中第1个数据为非截尾(已死亡)则由(22.11)式得 w=1/26-1=-0.96第2,3数据都是截尾数据(存活),它处于序号i3=1与2之间据(22.12)式,有 w=1/26=0.04第4个数据为非截尾,则由(22.11)式 w=1/26+1/23-1=-0.92余类推(见表22.3第6列)6.计算所指定的组别(本例为中药组,序号为i2)的logrank变换值之和t t=w i2上式的连加系在指定的i2范围内相加。其均数与方差分别为 e(t)=m1/n w v(t)=m1m2/n(n-1)(w-e(t)/m1)2(22.14)式中连加系在全部观察值上完成,m1系指所指定的组别的例数,(本例

22、为中药组m1=16),n为总例数。 z=t-e(t)/v(t) z服从标准正态分布,故可由z0.05=1.96,z0.01=2.58作出统计推断。 本例资料有t=3.822,e(t)=0.440210-6,v(t)=3.1755,z=2.145,故p0.05, 拒绝h0,认为两种疗法生存分布不相同。3.2 cox-mantel检验(cox-mantel test) 又称广义savage检验(generalized savage test),可用于两个或多个生存分布的比较。仍用例22.1的资料说明本检验的计算过程,为叙述方便现将中药组称为a组,对照组称为b组。本检验的h0,h1同前。表22.4

23、例22.1资料两疗法生存分布比较的cox-mantel检验时间 状态 组别 a组 b组 合并死亡率 期望死亡数 ti si 期初人数 死亡数 期初人数 死亡数 (b组) n1i d1i n2i d2i pi =(+)/(+) = 1 死 b 16 10 1 0.038462 0.38462 2 活 b 16 9 2 活 a 16 9 3 死 b 15 8 1 0.043478 0.347826 4 死 a 15 1 7 0.045455 0.318181 6 死 b 14 7 1 0.047619 0.333333 6 活 a 14 6 6 活 a 14 6 7 死 b 12 6 1 0.05

24、5556 0.333333 7 活 b 12 6 8 活 a 12 4 9 死 a 11 1 4 0.066667 0.266667 9 活 a 11 410 死 a 9 1 4 0.076923 0.30769211 死 b 8 4 1 0.083333 0.33333311 活 b 8 412 活 a 8 213 死 a,b 7 1 2 1 0.222222 0.44444417 死 b 6 1 1 0.142857 0.14285718 死 a 6 1 0 0.166667 019 活 a 5 024 死 a 4 1 0 0.250000 026 死 a 3 1 0 0.333333 0

25、31 死 a 2 1 0 0.500000 043 活 a 1 0 小计 3.212284计算步骤为1.将两组生存数据混合由小到大排列,当截尾数据与非截尾数据数值相同时,截尾数据排列在后。并指明各生存数据的状态(死或活)及所属组别(见表22.4中1-3列)2.列出a,b两组各生存时间上的期初人数及死亡人数分别以n1i,d1i,n2i,d 2i表示(见表22.4中第4-7列)。3.在完全数据的相应行中计算合并死亡率pi pi=(d1i+d2i)/(n1i+n2i) (22.16)4.在两组中任选一组(本处用b组)计算 各生存时间点上的期望死亡人数,它由该组期初人数乘以合并死亡率而得 e(d2i)

26、=n2ipi (22.17)参见表22.4中8.9两列。5.所指定的组别(本处可b组)死亡人数的期望值与方差为 e(d2i)=n2ipi (22.18) v(d2i)=n1in2i/(n1i+n2i-1)pi(1-pi) (22.19)计算服从标准正态分布的统计量z z=d2i-e(d2i)/v(d2i) (22.20)可据z0.05=1.96,z0.01=2.58,作出统计推断。代入本例资料有d2i=7, e(d2i)=3.212284,v(d2i)=1.916190, z=2.7363 故p0.01,拒绝h0,认为两种疗法的生存期不相同。3.3 广义wilcoxon检验(generaliz

27、ed wilcoxon test) 又称为breslow检验法(breslow test),可用于两个或两个以上生存分布的比较,其h0,h1同前。本处仍用例22.1资料说明其计算过程。表22.5 例22.1资料两疗法生存分布比较的广义wilcoxon检验时间 状态 组别 累积生存率 计分值 a组计分值 ti si s(ti) ui ui(a) 1 死 b 0.9615 0.9615 2 活 b -0.0385 2 活 a -0.0385 -0.0385 3 死 b 0.9197 0.8812 4 死 a 0.8779 0.7976 0.7976 6 死 b 0.8361 0.714 6 活 a

28、 -0.1639 -0.1639 6 活 a -0.1639 -0.1639 7 死 b 0.7897 0.6258 7 活 b -0.2103 8 活 a -0.2103 -0.2103 9 死 a 0.7370 0.5267 0.5267 9 活 a -0.2630 -0.263010 死 a 0.6803 0.4173 0.417311 死 b 0.6236 0.303911 活 b -0.376412 活 a -0.3764 -0.376413 死 a 0.4851 0.1087 0.108713 死 b 0.4851 0.108717 死 b 0.4158 -0.099118 死 a

29、 0.3465 -0.2377 -0.237719 活 a -0.6535 -0.653524 死 a 0.2598 -0.3937 -0.393726 死 a 0.1732 -0.567 -0.567031 死 a 0.0866 -0.7402 -0.740243 活 a -0.9134 -0.91341.将两组生存数据混合由小到大排列,当截尾数据与完全数据数值相同时,截尾数据排列在后。并写出每个生存数据的状态(死或活)及所属组别(见表22.5中第1-3列)。2.用积限估计法对两组合并资料估计生存率(即累积生存率)s(ti),列于表中第4列3.计算各生存时间点的计分值ui。 ui= s(ti

30、-1)+s(ti)-1 观察值为完全数据 s(ti)-1 观察值为截尾数据 (22.21)其中s(0)=1例如第一个时间点ti=1是完全数据(死亡),故u1=1+0.9615-1=0.9615,第2,3个时间点为截尾数据 u2=u3=0.9615-1=-0.0385第4个时间点为完全数据 u4=0.9615+0.9197-1=0.8812 余类推。4.计算任一组的计分值之和的绝对值,本处选择a组,其计分值已记第6列中, t=ui(a) (22.22)t的期望值为0,方差为 v(t)=m1m2u2/(m1+m2)(m1+m2-1) (22.23)(22.23)式中u2系指全部生存时间点的u值平方和, z=t/v(t) (22.24)z服从标准正态分布,故可据z0.05=1.96,z0.01=2.58作出统计推断。本例资料有t=2.8712,u2=6.6559,v(t)=1.6384,z=2.243,p0表示该协变量是危险因素，越大使生存时间越短 i 0表示该协变量是保护因素，越大使生存时间越长 h(t,x)=h0(t)exp(1x1 + 2x2 + mxm ） h0(t)为基础风险函数，它