11.生存分析与SAS程序.doc

第十一章 生存分析与 SAS 程序 在医学研究中 考察研究因素对研究对象的效应 经典的研究设计方法只观察其所出现 的结局或结果 但是在一些诸如癌症等险恶疾病和慢性病的随访研究中 只观察其结局 就不够全面 还需要考察研究对象出现某种结局所经历的时间 这段时间称为生存时间 survival time 生存分析是用来分析生存时间资料的统计方法 是近几十年来产生且 发展甚为迅速的一门应用统计的分支 到目前为止 生存分析已形成一套完整的体系 包括描述生存规律和进行组间比较的参数和非参数方法以及分析影响生存期因素的回归 模型等 本章主要介绍描述生存规律的寿命表和乘积限估计法 用于组间比较的 logrank 检验和 作预后因素分析的 Cox 回归模型以及相应的 SAS 过程 11 1 生存分析的基本概念 11 1 1 生存时间生存时间 在医学研究中 对于肿瘤 心血管等慢性疾病 要考察其治疗方法优劣 疾病预后的好 坏以及影响疾病预后的因素 通常采用随访研究的方法 对某一疾病作随访研究时 一 般是从某一时间开始 观察到某一规定时间截止 而研究对象是始点以后陆续进入观察 随访中要规定一个事件作为随访结局 例如 病人死于研究疾病 如果病人的随访结果 是规定的结局 则称为失效或死亡 Failure 那么病人从进入观察随访到规定的结局出 现 其间所经历的这段时间称为生存时间或失效时间 Survival time Failure time Waiting time 通常用 Ti 表示 i 1 2 n n 为观察个体数 在随访中 由于客观条件限制 不能也不可能将全部观察对象都观察到规定的结局 Failure 在观察截止时 一组观察对象除了出现规定的结局外 还有三种结果 1 死于其他疾病 2 由于迁移等原因失去联系 失访 3 随访截止时尚未出 现规定的结局 这三种结果虽原因不一 但提供的信息是一致的 这类病人不能获得确 切的生存时间 而只知道其生存期比随访观察到的时间长的信息 我们将此类病例数据 称为截尾 Censored 数据 其观察到的时间用 Ci表示 在获得生存时间时要注意以下几点 1 随访结局可以不是死亡 随访的结局根据研究问题的性质和目的的不同可以是死亡 复发 恢复等 但要明确时间界线 例如 肿瘤病人已作了手术切除治疗 我们随访观 察他有否复发 复发就是我们规定的结局 2 确定进入随访观察的起点时间 随访研究中 无论是同时进入观察还是陆续进入观察 对所有对象都要规定某一事件作为进入观察起点事件 例如 病人确诊 治疗开始或手 术时间 动物试验中染毒等 作为起点事件要有明确的时间界线 3 时间尺度 医学研究一般用日历时间作为时间尺度 但在某些医学问题研究中 可能 不以日历时间作为时间尺度更能反映所研究的问题 例如 对儿童龋齿出现的随访研究 中常以儿童出到第几颗牙才出现龋齿作为时间尺度 11 1 2 生存数据的特征生存数据的特征 生存数据在结构上有它自己的特点 主要表现为 1 生存时间数据的分布与常见的统计数据分布有明显不同 常呈指数分布 Weibull 分 布 对数正态分布 对数 Logistic 分布 Gamma 分布或更为复杂的分布 生物医学中的 生存时间分布有时呈现不规则状态 因而难以用传统的统计方法对这类数据进行处理 2 生存数据中往往包含有截尾数据 在随访研究中 若能观察到研究所规定的结局就能 获得确切的生存时间 这类数据称为完全数据 但在实际工作中 由于时间和客观条件 限制 部分患者难以观察到终点事件 以致不能获得确切的生存时间 而只能获得进入 观察至失访时这段时间 在分析总结时若剔除这部分对象往往导致样本有偏性 并且这 部分数据仍然提供了其实际生存时间大于观察到的截尾时间的信息 应充分利用这部分 资料的信息 因此 一般随访研究的样本数据中不可避免地包含有截尾数据 11 1 3 生存时间函数生存时间函数 描述生存时间分布规律的函数主要有生存函数 死亡概率函数 概率密度函数和危险率 函数 为了后文叙述方便 这里主要介绍生存函数和危险率函数 1 生存函数 在描述生存规律的数量指标中 以往常用的指标是某个特定时间的生存率 例如 3 年生存率 5 年生存率 这一指标的主要缺陷为不能反映整个生存规律 一个理想的指 标应该是任意时间的生存率 即生存率是任意时刻 t 的函数 其意义是研究个体生存时 间长于 t 的概率 若令 T 为生存期 s t 为任意时刻 t 的生存率 得 s t p T t 0 t 20 05 1 3 84 P 0 05 两组样本生存率间差异有统计学意义 可以认 为用猪苓提取物合用化疗治疗急性白血病的生存期大于对照疗法 表表 11 3 对数秩和检验计算表对数秩和检验计算表 编 号 时间 用猪苓组 第 1 组 对照组 第 2 组 合计 itin1id1iw1ie1iV1in2id2iw2ie2iV2inidi 1 2 3 4 5 6 7 8 9 10 11 12 13 1 1 516000 6150 236710100 3850 2367261 22 16010 00 0 9010 00 0250 33 515000 6520 2268 8100 3480 2268231 4415100 6820 2169 7000 3180 2169221 5614000 6670 2381 7100 3330 2381211 66 14020 00 0 6000 00 0200 76 512000 6670 2222 6100 3330 2222181 86 5 12000 00 0 5010 00 0170 97 5 12010 00 0 4000 00 0160 108 511100 7330 1956 4000 2670 1956151 119 10010 00 0 4000 00 0140 12109 100 6920 2130 4000 3080 2130131 13118 000 6670 2222 4100 3330 2222121 1411 8 000 00 0 3010 00 0110 1512 8 010 00 3025 2000 00 0100 16137 101 5560 1224 2100 4440 3025 92 17176 000 8570 0 1100 1430 1224 71 18186 101 0000 0 0000 0 61 1919 5 010 00 0 0000 0 50 20244 101 0000 0 0 000 0 41 21263 101 0000 0 0000 0 31 22312 101 0000 0 0000 0 21 2343 1 010 0 0000 0 10 合计811 7882 196473 2122 1964 精确法 2统计量的计算公式为 2 Dg Eg 2 vg 11 16 式中 vg为第 g 组的 Fg的方差估计值 其计算公式为 vg 11 17 1n n dn d nn n i 2 i iiigiigi 所计算的 v1 v2 精确法 2值服从自由度为 1 的 2分布 本例先求 vgi 再求 vg 得 v1 v2 2 1964 2统计量为 l 1i gi 2 6 533 专业结论同近似法 1964 2 788 118 2 对于多组生存率比较 其近似法 2统计量的计算方法与两组相同 设有 m 组 g 1 2 m 这时自由度 m 1 关于多组 Logrank 检验的精确法 2统计量计算较为复 杂 其计算公式为 2 s v 1s m 1 11 18 其中 s s1 s2 sg sm 1 s 为向量 s 的转置 Sg的计算公式为 sg Og Eg g 1 m 1 11 19 l i igiigi nndd 1 V 为 m 1 m 1 矩阵 记为 V Vgh m 1 m 1 其中 Vgh为第 g 组与第 k 组之间的方差与 协方差 计算公式为 Vgh 11 20 l i iiiiiigihghihi nndndnnnn 1 2 1 当 g h 时 gh 1 当 g h 时 gh 0 从多组生存率比较的计算公式可知 当 m 2 时 即为两组生存率的计算公式 2 Wilcoxon 检验法及与 Logrank 检验法比较 1 Wilcoxon 检验法 当 g 1 2 g m 时 Wilcoxon 检验法 2统计量计算公式仍可表示为 2 s v 1s m 1 其中 s s1 s2 sg sm 1 s 为向量 s 的转置 Sg的计算公式为 sg 11 21 l i iigiigi nnddw 1 V 为 m 1 m 1 矩阵 记为 V Vgh m 1 m 1 Vgh的计算公式为 Vgh 11 22 l i iiiiiigihghihii nndndnnnnw 1 22 1 上面 sg和 vgh计算公式中 wi为权重 这里 wi nI 对于例 11 3 资料 若用 Wilcoxon 检验法 得 2 5 2822 结论同 Logrank 检验 2 两种检验方法的比较 Logrank 检验法和 Wilcoxon 检验法实际上可以用统一公式来表示 即 Wilcoxon 检验法 的公式 公式中的权重 wi 1 时为 Logrank 检验法 wi ni时为 Wilcoxon 检验 因而可以 发现 Logrank 检验对生存时间较长的个体在检验中权重较大 对生存时间较短的个体 在检验中权重较小 在生存率 曲线 比较中 这种方法对尾部较为敏感 而 Wilcoxon 检验则与 Logrank 检验相反 对生存时间较短的个体在检验中权重较大 比较中对数据 的头部的差别较为敏感 从理论和实践中均发现 当生存资料的各死亡点的危险率在两 组或多组间成比例时 Logrank 检验的效率高于 Wilcoxon 检验法 宜选用 Logrank 检验 当生存资料各时点的危险率服从其他状态时 Wilcoxon 检验法效率高于 Logrank 检验 宜选用 Wilcoxon 检验法 11 2 3 应用实例应用实例 例例 11 4 对例例 11 3 的资料 试估计两组的生存函数 并用 Logrank 检验和 Wilcoxon 检 验法作比较 1 SAS 程序 DATA SA INPUT T IF T16 THEN GROUP 2 ELSE GROUP 1 T ABS T CARDS 2 4 6 6 7 5 8 5 9 10 12 13 18 19 24 26 31 43 1 5 2 3 5 6 6 5 6 5 11 11 13 17 PROC LIFETEST METHOD PL PLOT S TIME T CENSOR 1 STRATA GROUP RUN 2 程序说明 数据输入时截尾数据用负值表示 第一个 if 语句产生指示变量 censor 其取值为 1 时为 截尾数据 取值为 0 时为完全数据 第 2 个 if 语句产生分组变量 group 前 16 个数值属 于用猪苓组 1 组 后 10 个数据属于对照组 2 组 对 t 取绝对值是为了保证参与计 算的生存时间都是正值 PROC LIFETEST 过程可选择寿命表法或乘积限估计法对生存资料作生存率估计 并可 对两组或多组间的生存率作 Logrank 检验和 Wilcoxon 检验 method 选项有两种选择 即 PL 和 life PL 选项为隐含值 表示用乘积限估计法估计生存率 Life 表示用寿命表 法估计生存率 当用寿命表法分析时 程序自动形成生存时间的分组区间 也可人为指 定分组区间 在 proc 语句中加上 intervals a to b by c a b c 分别表示初值 终值和 步长 Plots 要求绘图 其中 s 表示生存率曲线 L 表示取对数 H 表示危险率函数 图形的横纵坐标分别为 s t s LS t log s LLs log t log Log log s H t h time 语为 lifetest 过程中必须语句 为设置生存时间和截尾指示 当有分组变量时用 strate 表 示 如果资料中还包含有协变量 可在 num 语句前加 test 语句 如 test X1 X2 以便检验协 变量与生存时间联系的密切程度 如果用下一节介绍的 phreg 过程来揭示协变量与生存 时间的关系则更好 3 输出结果及解释 结果首先输出两组的乘积限法估计的生存率 survival 及死亡率 Failure 生存率的标 准误 Survival Standard Error 死亡数 Number Failed 和存活数 Number left 表示 有 号者为截尾观察值 输出中有生存时间的四分位数 第 50 位数为中位生存期 第 1 组为 24 个月 第 2 组 为 11 个月 两组平均生存期分别为 20 35 和 10 76 个月 由于第 1 组最后一个值为截尾 数据 所以均数的估计是有偏性的 The LIFETEST Procedure Product Limit Survival Estimates GROUP 1 Survival Standard Number Number T Survival Failure Error Failed Left 0 0000 1 0000 0 0 0 16 2 0000 0 15 4 0000 0 9333 0 0667 0 0644 1 14 6 0000 1 13 6 0000 1 12 7 5000 1 11 8 5000 0 8485 0 1515 0 0999 2 10 9 0000 2 9 10 0000 0 7542 0 2458 0 1256 3 8 12 0000 3 7 13 0000 0 6465 0 3535 0 1468 4 6 18 0000 0 5387 0 4613 0 1570 5 5 19 0000 5 4 24 0000 0 4040 0 5960 0 1657 6 3 26 0000 0 2694 0 7306 0 1559 7 2 31 0000 0 1347 0 8653 0 1231 8 1 43 0000 8 0 Censored Observation Summary Statistics for Time Variable T Point 95 Confidence Interval Quantile Estimate Lower Upper 75 31 0000 18 0000 50 24 0000 13 0000 31 0000 25 13 0000 8 5000 24 0000 Mean 20 3549 Standard Error 2 9640 NOTE The last observation was censored so the estimate of the mean is biased The LIFETEST Procedure Product Limit Survival Estimates GROUP 2 Survival Standard Number Number T Survival Failure Error Failed Left 0 0000 1 0000 0 0 0 10 1 5000 0 9000 0 1000 0 0949 1 9 2 0000 1 8 3 5000 1 7 6 0000 0 7714 0 2286 0 1442 2 6 6 5000 0 6429 0 3571 0 1679 3 5 6 5000 3 4 11 0000 0 4821 0 5179 0 1877 4 3 11 0000 4 2 13 0000 0 2411 0 7589 0 1946 5 1 17 0000 0 1 0000 0 6 0 Censored Observation Censored Observation Summary Statistics for Time Variable T Point 95 Confidence Interval Quantile Estimate Lower Upper 75 13 0000 11 0000 17 0000 50 11 0000 6 0000 17 0000 25 6 5000 1 5000 13 0000 Mean 10 7571 Standard Error 1 9434 Summary of the Number of Censored and Uncensored Values GROUP Total Failed Censored Censored 1 16 8 8 50 0000 2 10 6 4 40 0000 Total 26 14 12 46 1538 以下输出 Kaplan Meier 生存曲线图 The LIFETEST Procedure Survival Function Estimates SDF 1 0 A A A S B B u r A A v i 0 8 v BB a A A l D B B A A i s 0 6 t r A A i b B B u t i 0 4 A A o n F u A A n B B c 0 2 t i A o n 0 0 B 0 5 10 15 20 25 30 35 40 45 T Strata B B B B B A A A A A A A 0 5 10 15 20 25 30 35 40 45 T Legend for Strata Symbols A GROUP 1 B GROUP 2 最后输出 Logrank 检验和 Wilcoxon 检验过程中的秩统计量 S 和方差协方差 V 以 及 2统计量和 p 值 其结果与例 11 3 的手工计算结果一致 Testing Homogeneity of Survival Curves over Strata Time Variable T Rank Statistics GROUP Log Rank Wilcoxon 1 3 1355 46 000 2 3 1355 46 000 Covariance Matrix for the Log Rank Statistics GROUP 1 2 1 1 95379 1 95379 2 1 95379 1 95379 Covariance Matrix for the Wilcoxon Statistics GROUP 1 2 1 577 500 577 500 2 577 500 577 500 Test of Equality over Strata Pr Test Chi Square DF Chi Square Log Rank 5 0321 1 0 0249 Wilcoxon 3 6641 1 0 0556 2Log LR 2 2239 1 0 1359 11 3 Cox 回归模型与 PHREG 过程 1972 年 英国统计学家 DR Cox 提出了半参数生存分析数学模型 Cox 回归模型 Cox regression medel 在以后二十多年中 众多的生物统计学家在理论和应用方面作了大量 的研究 目前 Cox 回归模型已成为生存分析中理论较为完善 应用最广泛的统计模型 我国随着计算机和统计软件的推广应用也已开始应用 Cox 模型进行慢性病的预后分析 实践证明 Cox 模型对许多生存资料都有用 而且行之有效 本节简要介绍 Cox 模型结构 意义 参数估计和假设检验方法 以及 SAS 软件中如何 应用 PHREG 过程作 Cox 模型分析 11 3 1 Cox 模型的结构及流行病学意义模型的结构及流行病学意义 在随访研究中 要考察和比较不同的治疗方法 不同病理类型 病人的某些特征对疾病 预后的影响 可运用非参数方法进行组间比较 这类方法虽然使用简单 但在多因素共 存条件下 单一因素的比较会受到其他因素干扰和混杂 组间难以达到均衡 而且不能 分析和考察因素间的关系 如交互作用 和进行定量评价 因此 需用多因素回归模型 的方法 生存分析中 最典型 最常用的就是 Cox 回归模型 1 Cox 模型的结构 对于一组带有多个预后因素 称为解释变量或协变量 的随访病人 假设病人不受到任 何协变量的作用 那末 病人在整个随访期间会显现一定的生存或死亡规律 这一死亡 规律用 h0 t 来描述和表示 而且所有病人的死亡规律应该是一致的 而病人在随访期间 实际上是受到多个预后因素的综合作用 预后因素对病人的生存或死亡规律的影响可以 看作是协变量对 h0 t 的修改 并且病人各自所带有的预后因素是不同的 因而 病人实 际上呈现出不同的生存或死亡规律 病人的实际死亡规律用 h t x 表示 x 表示协变量 由此可见 生存分析回归模型由基本部分 h0 t 和修改部分 x x 为预后因素 x 的函 数 两部分组成 这里协变量 x x1 x2 xp可以是不同的治疗方法或手术方法 不同病理 类型或病人的性别 年龄等特征 也可以是复合变量 如交互作用 当基本部分 h0 t 与修改部分为乘积关系时 即危险函数由基本危险率函数乘上一个 x 协变量的常数因子组成 h t x h0 t 11 23 x 称此模型为成比例危险模型 Proportional hazards model 带协变量的修改部分最 x 常见的形式为 exp x 则模型为 h t x h0 t exp x h0 t exp 1x1 2x2 pxp 11 24 其中 1 2 p 是 p 个未知待估计的回归系数 它是描述各个因素对生存 期影响大小的参数向量 x x1 x2 xp 是一个 p 维协变量向量 xp可以是影响预后的治 疗方法 病人的某些特征 也可以是这些因素的交互作用项 h0 t 是未知的 任意分布 的基准危险函数 对比例危险模型 若假定 h0 t 为特定的时间分布函数 则称为全参数生存分析模型 可 以用极大似然法作参数 估计和检验 当 h0 t 未知分布时 用部分似然函数 Partial likelihood function 对参数 作出估计 则称此模型为半参数模型 也就是通常意义下 的 Cox 回归模型 2 流行病学意义 应用部分似然函数理论对模型中的 作出估计时 其意义在于 当某个协变量 xj为 二分变量时 xj 1 的危险性相对于 xj 0 为 RR 11 25 x 0 xexp t h x 1 xexp t h ppj110 ppj110 e j 从上式可以发现 xj 1 相对 xj 0 在任意时间的危险性是一个常数 也即两组病人的 e j 相对危险性为一个常数 成比例危险模型由此而得名 当 xj为多分类有序变量或连续变量时 为变量每增加一个单位 个体增加或减少的 e j 相对危险性 变量 xj是危险因素还是保护因素要视 j的正负值和变量的取值而定 11 3 2 Cox 模型的参数估计与统计推断模型的参数估计与统计推断 1 Cox 模型的似然函数 当生存数据中没有重合死亡 失效 点 无 ties 时 所有样本个体 i 1 2 n 的生 存时间 i 1 2 d 和截尾时间可以看作是同一始点开始随访观察 按从小到大排序记 为 t1 t2 tn 可以看作这一群个体须经过不同的死亡点 Failure 而在各死亡点必须有一 个个体死亡 失效 经过死亡点的所有个体记作死亡点 ti的危险集 R ti 已经生存到 ti 时刻的个体 R ti 在 ti时刻必须有一个个体死亡 而这个个体就是第 i 个体的概率 称 为条件死亡概率或条件危险率 记为 P i R ti 给定 R ti 且已知 ti时刻有一个病人死亡的条件下 则病人是第 i 个体的概率是 P i R ti hi t x 11 26 t Rk k i x t h 根据似然函数构造理论 Cox 模型的似然函数应是各个死亡点对它的贡献而构成 即各 死亡点 ti的条件危险率的乘积 L d 1i t Rk k i l x t h x t h d 1i t Rk kk i0 l x t h x t h d 1i t Rk k i l xexp xexp 有趣的是在有截尾数据时 其似然函数仅为死亡病人在各死亡点上条件死亡概率的乘积 这一似然函数并不是通常意义下的似然函数 DR Cox 于 1975 年证明了这是特定意义 下全似然函数中的一部分 故称部分似然函数 Partial likelihood fanction 2 参数 及其方差的估计 Cox 也证明了可将一般极大似然理论应用到部分似然函数的参数 及其方差的估计 为计算方便 对 Cox 模型的部分似然函数两边取对数 得对数部分似然函数 lnL 11 27 d 1i t Rk ki i xexp lnx 应用极大似然估计理论 将对数部分似然函数 lnL 对 j求一阶偏导数并令其等于 0 即 0 1 Lln 0 P 个非线性方程组 2 Lln 0 p Lln 对这一似然方程组 可用迭代法 例如 Newton Raphson 法 求 j的估计值 迭代求解 过程从略 参数 j的估计值的方差和标准误可用二阶偏导数的 j var j SE j Lln 2 负值组成的信息矩阵及其逆矩阵 方差协方差矩阵 的主对角线元素求 pp I pp 1 pp VI 得 3 参数 及回归效果检验 在 cox 模型的参数估计和模型建立过程中 要对变量进行筛选并作出检验 DR Cox 的 1972 和 1975 年的文章以及以后十多年中许多学者对部分似然函数的大样本性质作了 证明 提出部分似然可代替普通似然对参数 进行统计推断 关于极大部分似然估计 的统计推断有三种渐近方法 1 似然比检验 Likelihood ratio test 设原模型参数为 1 2 p 并求得对数部分的似然函数值 需检验的 Lln 参数 1 1 2 k k Variable DF Estimate Error Chi Square Chi Square RENAL 1 3 591054 1 07898 11 07687 0 0009 Analysis of Maximum Likelihood Estimates Risk Variable Ratio RENAL 36 272 Summary of Stepwise Procedure Variable Number Score Wald Pr Step Entered Removed In Chi Square Chi Square Chi Square 1 RENAL 1 27 0399 0 0001 Data Set WORK COX Dependent Variable T Censoring Variable D Censoring Value s 0 Ties Handling BRESLOW Summary of the Number of Event and Censored Values Percent Total Event Censored Censored 25 21 4 16 00 Step 1 Variable RENAL is entered The model contains the following explanatory variables RENAL Testing Global Null Hypothesis BETA 0 Without With Criterion Covariates Covariates Model Chi Square 2 LOG L 108 948 89 817 19 132 with 1 DF p 0 0001 Score 27 040 with 1 DF p 0 0001 Wald 11 077 with 1 DF p 0 0009 NOTE No additional variables met the 0 05 level for entry into the model Analysis of Maximum Likelihood Estimates Parameter Standard Wald Pr Variable DF Estimate Error Chi Square Chi Square RENAL 1 3 591054 1 07898 11 07687 0 0009 Analysis of Maximum Likelihood Estimates Risk Variable Ratio RENAL 36 272 Su