临床试验-肿瘤：从生存数据谈ADTTE的构建，介绍了生存数据的定义与删失，以及ADTTE的构建。这篇文章介绍，内容包括生存函数与生存函数的Kaplan-Meier估计，第5节有相关SAS程序分享。
生存数据的分析包括，生存过程的描述、生存过程的比较以及影响因素的分析，这篇文章属于生存过程的描述。相应的统计分析方法分为参数方法和非参数方法。参数方法需要对生存分布、微积分有一定的了解，学习门槛会比较高。在临床研究和流行病学的刊物上，大多数研究工作都是使用非参数方法。
文章主要参考以下2本书籍：
- Statistical Methods for Survival Data Analysis (ELISA T. LEE)，
- Clinical Statistics: Introducing Clinical Trials, Survival Analysis, and Longitudinal Data Analysis (Olga Korosteleva)。
欢迎关注，SAS茶谈，后台回复：SA，获取相关资料。
1. 生存数据的简单介绍
生存数据可以广义地定义为到给定事件发生的时间 (the time to the occurrence of a given event)，也叫做生存时间 (survival time)，这类数据也被称为Time-to-event。
在研究中，所关注的目标事件可能不会发生，这种情况下，对应的数据观测称之为删失观测或截尾观测。生存数据由两个变量构成：一个是表示时间长短的数值型变量，另一个是表示事件状态的分类变量（事件发生、事件删失）。
通常，研究会为删失观测确定一个删失时间来表示时间长短。如果简单地把删失时间纳入分析，那结果就会偏保守（实际生存时间长于观察到的时间）。针对这种情况，统计学家发明了特定的统计方法来进行分析，即生存分析。
2. 生存时间的函数
用 T 表示生存时间，T 的分布可以用三个互相等价的函数来刻画：
- 生存函数 (Survival Function)
- 概率密度函数 (Probability Density Function)
- 风险函数 (Hazard Function)
2.1 生存函数 (Survival Function)
生存函数用 S(t) 表示，定义为个体存活时间长于 t 的概率，其为 t 的函数。
生存函数 S(t) 是 t 的不增函数，在时间为0时，生存的概率为1；在时间趋向于无穷时，生存概率为0。
举一个简单例子，假设一项肿瘤研究入组100人，第1年有10人死亡，第2年有20人死亡。在没有删失的情况下，S(1) = (100 - 10) / 100 = 0.9，即生存时间长于1年的概率为0.9；S(2) = (100 -10 - 20) = 0.7，即生存时间长于2年的概率为0.7。
2.2 生存曲线 (survival curve)
生存函数 S(t) 又叫做累积生存率 (cumulative survival rate)，其图形称为生存曲线 (survival curve)，用于描述生存过程。
陡峭的生存曲线，表示低的生存率或较短的生存时间；平缓的生存曲线，表示高的生存率或较长的生存时间。
生存函数或生存曲线可以用来确定生存时间的中位数以及其他分位数 (50%, 25%, 75%)，也可以用来比较两个或多个生存分布。
2.3 中位生存时间 (median survival time)
中位生存时间，又称半数生存时间，表示50%个体存活且有50%个体死亡的时间，即当 S(t) = 0.5 时所对应的 t 值。在上图 (a) (b) 中，中位生存时间约等于 5 和 36。
通常，均值是用来描述一个分布的中心趋势。生存时间的分布常为偏态分布，少数特别长或特别短的生存时间会对均值产生过大影响，所以用中位数来描述生存时间的中心趋势会更合适一些。
2.4 生存函数的估计
总体的生存函数是未知的，我们需要基于样本数据对总体生存函数进行估计。
研究中，如果不存在删失数据，生存函数可使用存活时间长于 t 的患者数所占比例来估计：
举例来说，生存数据4, 6, 6, 10, 15, 20。当 t = 4 时，生存函数估计为 5 / 6 = 0.833， 存活时间长于 4 的概率为0.833；当 t = 6 时，生存函数估计为 3 / 6 = 0.5，存活时间长于 6 的概率为0.5。
当有删失数据时，就无法用上式进行估计，因为存活时间长于 t 的患者数有时是不能确定的。
举例来说，生存数据4, 6, 6+, 10+, 15, 20，"+"表示删失。当 t = 5时，我们可以获取生存函数的估计 5 / 6 = 0.833；但是我们不能得到 t = 11 时的估计，因为删失值10+的存在使得存活时间长于11的患者数是不确定的。
所以，当删失值存在时，需要用其他方法来进行生存函数的估计。估计生存函数最常用的非参数方法，是Kaplan和Meier于1958年提出的乘积-极限法 (product-limit, PL)。
3. Kaplan-Meier估计
一句话概括，生存函数的Kaplan-Meier估计为，各时段生存率的乘积。
何为生存率？该时段处于风险下 (at risk) 的受试者中，存活受试者所占的比例。何为处于风险下 (at risk) ？上一时段存活且没有删失的进入当前时段的受试者。如果某一时段含有删失值，删失值会计入该时段的存活数，但是不会计入到下一时段处于风险下的人数以及存活数中。
假设有 k 个不同的生存时间观测，按增序排列，t1 < t2 < … < tk。在 ti 时段 ( ti-1< t ≤ ti )，处于风险下的受试者数记为 ni (不包括上一时段的删失受试者)，死亡人数记为 di，则活过ti时段的人数为 ni - di，ti时段的生存率为 (ni - di) / ni。生存函数的KM估计，如下表示：
从条件概率角度看，KM估计可以这样理解，生存时间长于 k 的概率 = 第1年受试者的生存率第1年存活受试者在第2年的生存率
······
第k-1年存活受试者在第k年的生存率。
那么两相邻时刻的生存函数有如下关系：
定义与公式比较抽象，下面以数据进行演示，生存数据4, 6, 6+, 10+, 15, 20，"+"表示删失。
当 t = 0 时，处于风险下的受试者数为6，死亡人数为0，删失人数为0，生存率为(6 - 0) / 6 = 1，生存函数估计为1，即存活时间长于0的概率为1；
当 t = 4 时（0 - 4时段），当前时段处于风险下的受试者数为6，死亡人数为1，删失人数为0，生存率为(6 - 1) / 6 = 0.833，生存函数估计为1
0.833 = 0.833，即存活时间长于4 的概率为0.833；
当 t = 6 时（4 - 6时段），前一时段有1人死亡，0人删失，当前时段处于风险下的受试者数为5，死亡人数为1，删失人数为1，生存率为(5 - 1) / 5 = 0.8，生存函数估计为0.833
当 t = 10 时（6 - 10时段），前一时段有1人死亡，1人删失，当前时段处于风险下的受试者数为3，死亡人数为0，删失人数为1，生存率为(3 - 0) / 3 = 1，生存函数估计为0.833
当 t = 15 时（10 - 15时段），前一时段有0人死亡，1人删失，当前时段处于风险下的受试者数为2，死亡人数为1，删失人数为0，生存率为(2 - 1) / 2 = 0.5，生存函数估计为0.833
当 t = 20时（15 - 20时段），前一时段有1人死亡，0人删失，当前时段处于风险下的受试者数为1，死亡人数为1，删失人数为0，生存率为(1 - 1) / 1 = 0，生存函数估计为0.833
以上生存数据对应的生存曲线图如下：
4. 生存数据实例
一项为期 2 年的临床试验，测试新心脏瓣膜的功效。10 名患者的生存时间如下（以月为单位），“+”表示结果已删失。
当 t = 0 时，处于风险下的受试者数为10，死亡人数为0，删失人数为0，生存率为(10- 0) / 10 = 1，KM估计为1，即存活时间长于0的概率为1；
当 t = 5 时（0 - 5时段），当前时段处于风险下的受试者数为10，死亡人数为1，删失人数为0，生存率为(10 - 1) / 10 = 0.9，KM估计为1
当 t = 8 时（5 - 8时段），当前时段处于风险下的受试者数为9，死亡人数为1，删失人数为1，生存率为(9 - 1) / 9 = 0.89，KM估计为1
当 t = 10 时（8 - 10时段），当前时段处于风险下的受试者数为7，死亡人数为2，删失人数为0，生存率为(7 - 2) / 7 = 0.71，KM估计为1
其他各时间点KM估计参考以下表格。
5. Kaplan-Meier估计的SAS程序
通过SAS中Proc Lifetest进行生存分析，KM估计以及对应的生存曲线都可以进行简单输出。不过，如果想要输出更复杂的图形，就需要借助SGPLOT或GTL语句进行实现，这一部分会另写文章进行介绍。
以上一节研究数据为例，下面分享一些简单的生存分析SAS语句。
5.1 数据输入
按照ADTTE IG 建议，以数值1表示数据删失，数值0表示数据未删失。数据输入程序如下：
**Data input;
data valves;
input duration status @@;
datalines;
24 1 16 1 8 0 19 0 10 0 8 1 5 0 17 0 20 0 10 0
;
run;
5.2 KM估计
过程步Proc Lifetest进行KM估计，需要指定删失状态的变量值。变量值可以是单个数值，也可以是值的序列。
**KM estimator;
proc lifetest data = valves method = km;
time duration * status(1);
run;
以下为SAS默认输出结果，Survival列为KM的估计值，SAS不会输出删失记录以及重复记录的KM估计。
5.3 KM Plot
若想要简单输出KM估计的生存曲线，即KM Plot， 可使用plots=选项，并打开ods graphics选项。
**KM plot;
ods graphics on;
proc lifetest data = valves method = km plots=( survival );
time duration * status(1);
run;
ods graphics off;
输出图形如下：
5.4 处于风险下的人数 (At risk)
若想要输出处于风险下的人数，需要使用atrisk选项。
**At risk;
ods graphics on;
proc lifetest data = valves method = km plots=( survival(atrisk) );
time duration * status(1);
run;
ods graphics off;
输出图形如下：
5.5 各数据时点的处于风险下的人数
SAS默认输出的时间轴刻度与数据中的时点并不一致，可以使用atrisktickonly选项以及设置ticks序列范围，使之保持一致。
**Ticks in data;
ods graphics on;
proc lifetest data = valves method = km plots=( survival( atrisk(atrisktickonly) =(0 5 6 10 16 17 19 20 24) ) );
time duration * status(1);
run;
ods graphics off;
输出图形如下：
其他处理，参考Proc Lifetest的SAS官方文档。
6. 总结
这篇文章介绍了，生存分析中的生存函数，以及对应的非参数估计——Kaplan-Meier估计。KM估计为生存数据各时段的生存率的乘积，读者可以通过文章中具体时点的估计值解说，加深对KM估计的理解。
感谢阅读， 欢迎关注：SAS茶谈！
若有疑问，欢迎评论交流！
梳理不易，转载请注明出处 (by Jihai / SAS茶谈)
