温馨提示：
文中链接在微信中无法生效。请点击底部
「阅读原文」
。或直接长按/扫描如下二维码，直达原文：
作者：
雷诺 (新加坡国立大学)
邮箱：
[email protected]
Source:
Mark Bounthavong. (2023). “Survival Analysis in R.” February 13, 2023. -Link-.
Keywords
: 生存函数, 生存曲线, 对数秩检验, cox, duration analysis, survival analysis, time-to-event, Kaplan-Meier, survdiff, coxph
1. 引言
在实证分析中，我们常用二元变量来描述两种结果，如“是或否”，“存活或死亡”，“治愈或未治愈”。我们还可以通过相对风险比或比值比来比较两组在特定事件发生方面的概率或风险差异。尽管这些方法能提供关于效应的大小、方向以及确定性（例如，95% 置信区间）的信息，但它们并不提供关于个体存活或处于风险状态的时间长度这一关键信息。为了解决这一问题，我们需要运用生存分析。
生存分析（survival analysis），亦称时至事件分析（time-to-event analysis）或持续时间分析（duration analysis），是一种专门研究事件发生之前时间长度的统计方法。与仅确定事件发生与否相比，它更关注时间维度的分析。生存数据通常不符合正态分布，并且经常含有不完整信息，比如存在被删失的受试者。这些独特的数据特性要求我们采用特定的方法和统计技巧来进行综合分析。
在社会经济学研究领域，生存分析被广泛应用于探究一系列复杂问题，如失业与就业的持续时间、通货膨胀的影响周期、银行贷款的供求关系，以及产品、供应商和顾客预期寿命的分析等。它旨在解答一系列问题：例如，在总体中，有多少比例的个体能够存活超过一定时间？在存活的样本中，他们会以何种速率经历死亡或其他失败事件？是否可以同时考虑多个导致死亡或失败的原因？特定情境或特征如何影响生存概率？通过这些问题的探索，生存分析帮助我们深入理解并分析经济行为和社会趋势。
2. 相关概念和模型
2.1 生存函数和危险函数 (survivor and Hazard functions)
生存分析的两个主要特征是生存函数 (survivor function) 和危险函数 (hazard function)。生存函数用于描述在特定时间点之后，主体继续存活的概率，表示为：
其中，
是生存概率 (
)，
是事件
的时间大于某一时间
的概率。
危险函数用于描述给定主体仍然存活的条件下，单位时间内事件发生的瞬时潜在可能性（率），表示为：
具体而言，在某一瞬间点
上的危险率是指，在主体存活的条件下 (
)，事件在极短的未来时间内发生的概率，即
。
生存函数和危险函数的关系可表示为：
这表明，我们可以从生存函数
推导出危险函数
，反之亦然。
2.2 Kaplan-Meier 生存曲线 (Kaplan-Meirer Curve)
Kaplan-Meier 生存曲线（KM 曲线）是一种非参数方法，可以直观地比较不同组别的生存情况。这种曲线的特点是其“阶梯状”外观，每一个“台阶”表示受试者在某个时间点发生重要事件（如死亡）或数据被“删失” (censoring)。所谓“删失”，指的是某些情况下无法获得完整数据，比如某人在研究中途退出或在特定时间点还未经历关键事件（如75岁时仍然存活）。
要绘制 KM 曲线，必须至少包含以下变量：
在本例中，我们有两个不同的组别（group），每组包含10名受试者（subject）。如果受试者在研究期间死亡，则
，若受试者在研究期间退出或失联，则事件状态标记为删失，即
。“生存时间” (
) 指的是受试者在风险中度过或对生存时间作出贡献的具体时间长度。本例中的数据可以从 GitHub site 中获得。
kmplot "kmcurve.csv")"Table 1. Variables for constructing KM curves" )Table: Table 1. Variables for constructing KM curves1 | 10 | 1 | 1 |2 | 20 | 1 | 1 |3 | 30 | 1 | 1 |4 | 35 | 0 | 1 |5 | 40 | 1 | 1 |6 | 50 | 1 | 1 |7 | 60 | 0 | 1 |8 | 70 | 1 | 1 |9 | 80 | 1 | 1 |10 | 90 | 1 | 1 |11 | 5 | 1 | 0 |12 | 8 | 1 | 0 |13 | 14 | 1 | 0 |14 | 20 | 1 | 0 |15 | 22 | 0 | 0 |16 | 23 | 0 | 0 |17 | 24 | 0 | 0 |18 | 25 | 1 | 0 |19 | 45 | 1 | 0 |20 | 60 | 1 | 0 |
从图 1 中，我们可以发现，随着时间的推移，每个时间点的生存可能性（survival probability）逐渐下降。例如，
时，每组有 10 名受试者。
时，处理组 0 的受试者减少至 3 人，而处理组 1 仍有 5 人。
时，处理组 0 的受试者降至 0，而处理组 1 则剩下 2 人。图中的刻度标记指示了各个时间点的受试者删失数量。在此示例中，处理组 0 有 3 个刻度标记，而处理组 1 有 2 个。
library (survival)library (survminer)TRUE ,"treatment 0" , "treatment 1" ),"Kaplan-Meier curves" )
图 1 KM 曲线
2.3 对数秩检验 (Log rank test)
对数秩检验（log rank test）是一种统计假设检验，用于比较两组样本在生存时间上的分布差异。作为非参数检验，它特别适用于数据右偏和存在删失的情况。虽然对数秩检验能够检测生存曲线之间的差异，但它不提供关于效应大小、方向或不确定性的具体信息。对数秩检验的基本统计假设包括：
: 两个群体在任何时间点上事件发生概率上没有差异。
: 两个群体在任何时间点上事件发生概率上存在差异。
在统计显著性水平设定为
的前提下，如果 p 值小于 0.05，我们将拒绝零假设，即认为两个群体间的生存情况存在显著差异。相反，如果 p 值大于或等于 0.05，我们将无法拒绝零假设，这意味着没有足够的证据支持两个群体间的生存情况存在显著性差异。接下来，我们用
survival
包中的
survdiff
函数进行对数秩检验：
survdiff(Surv(survivaltime, event) ~ group, data = kmplot)## Call: ## survdiff(formula = Surv(survivaltime, event) ~ group, data = kmplot) ## ## N Observed Expected (O-E)^2/E (O-E)^2/V ## group=0 10 7 4.14 1.978 3.27 ## group=1 10 8 10.86 0.754 3.27 ## Chisq= 3.3 on 1 degrees of freedom, p= 0.07
正如卡方检验一样，对数秩检验通过比较预期和观察到的死亡受试者数量来计算 p 值。本例中的 p 值是 0.07，意味着我们没有足够的证据拒绝零假设，即我们不能断定处理组 0 和处理组 1 之间的生存率存在统计学上的显著差异。如前文所述，由于对数秩检验是一个显著性检验，它并不提供关于处理效果大小、方向和不确定性的信息。为此，为了更深入地分析处理效果，我们需要构建一个 Cox 比例风险模型。
2.4 Cox 比例风险模型 (Cox proportional hazards model)
Cox 比例风险回归模型 (Cox proportional hazards model) 是一种生存回归，它假设两组之间的风险是成比例的（比例风险回归）。这意味着预测变量对结果（生存）的影响在时间上是恒定的，并且是线性累加的。例如，如果某因素在
时将风险增加一倍，那么在
、
等也会增加一倍。确定是否违背比例风险假设的一种方法是检查两条 KM 曲线是否在多个时间点相交。
Cox 比例风险模型没有常数项，而是以一个随时间变化的基线风险
作为基础，可以表示为：
风险比 (
) 是两组风险的指数，可以表示为：
接下来，我们将利用
survival
包中的
coxph
函数来估计 Cox 比例风险模型。为此，我们重写风险比 (
) 公式，以反映主要的预测变量，即不同的处理组（处理组 0 vs 处理组 1）：
这里，我们感兴趣的系数是
。
coxmodel coxmodel ## 生成对数风险比 ## Call: ## coxph(formula = Surv(survivaltime, event) ~ group, data = kmplot) ## ## coef exp(coef) se(coef) z p ## group -1.0554 0.3481 0.6084 -1.735 0.0828 ## ## Likelihood ratio test=3.06 on 1 df, p=0.08047 ## n= 20, number of events= 15
我们需要对对数风险比取指数来得到风险比，以及相应的 95% 置信区间：
exp(coef(coxmodel)) ## 生成指数系数 (或风险比) ## group ## 0.34805 exp(confint(coxmodel)) ## 生成 95% 置信区间 ## 2.5 % 97.5 % ## group 0.1056181 1.146952
这里我们感兴趣的是组系数，即对数风险比。当我们对其取指数时，我们得到的风险比为 0.35。95% 置信区间为 0.11 到 1.15。由此可知，与处理组 0 相比，处理组 1 中受试者发生事件的风险降低了 65%，其置信区间同样为 0.11 到 1.15。由于 95% 置信区间包含了零值 (
) 或
（无效果），我们无法判定该结果在统计学上具有显著性。
3. R 案例：老兵肺癌研究
3.1 数据
本节使用 Prentice (1973) 文中的数据来演示上文介绍的主要方法和模型在 R 中的实操。
该研究涵盖了 137 名晚期、无法进行手术的肺癌患者。这些患者被随机分配到两组：一组接受化疗加上新药治疗（即实验组），另一组仅接受化疗（即对照组）。需要注意的是，本例中 treatment 1 是指对照组， treatment 2 是指实验组。
首先，观察老兵数据的变量：
## 加载老兵数据，并展示其变量列表 "Table 2. Variables in the Veterans Affairs Lung Cancer Trial" 2. Variables in the Veterans Affairs Lung Cancer Trial## |variable |description | ## |:--------|:-------------------------------------------------------------------------| ## |trt |treatment assignment (1 = control, 2 = experimental) | ## |celltype |lung cancer cell type (1 = squamous, 2 = smallcell, 3 = adeno, 4 = large) | ## |time |survival time in days | ## |status |censoring status (1 = death occurred, 0 = censored) | ## |karno |Karnofsky performance score (100 = good) | ## |diagtime |months from diagnosis to randomization | ## |age |in years | ## |prior |prior therapy (0 = no, 10 = yes) |
观察前六行数据：
### 从 "survival" 包中加载老兵数据 "Table 3. First six rows of the Veterans Affairs Lung Cancer Trial" ## Table: Table 3. First six rows of the Veterans Affairs Lung Cancer Trial ## ## | trt|celltype | time| status| karno| diagtime| age| prior| ## |---:|:--------|----:|------:|-----:|--------:|---:|-----:| ## | 1|squamous | 72| 1| 60| 7| 69| 0| ## | 1|squamous | 411| 1| 70| 5| 64| 10| ## | 1|squamous | 228| 1| 60| 3| 38| 0| ## | 1|squamous | 126| 1| 60| 9| 63| 10| ## | 1|squamous | 118| 1| 70| 11| 65| 10| ## | 1|squamous | 10| 1| 20| 5| 49| 0|
从老兵数据集的汇总统计信息可知：该样本中包含 137 名受试者，他们的平均生存时间是 121.63 天，平均年龄是 58.3 岁。至于其它的离散变量，由于它们不适合用平均值来描述，因此我们需要采用
gmodels
包提供的
CrossTable
命令来进行统计性描述：
## 描述所有变量 library (psych)2 )) ## | | vars| n| mean| sd| median| trimmed| mad| min| max| range| skew| kurtosis| se| ## |:---------|----:|---:|------:|------:|------:|-------:|-----:|---:|---:|-----:|-----:|--------:|-----:| ## |trt | 1| 137| 1.50| 0.50| 1| 1.50| 0.00| 1| 2| 1| 0.01| -2.01| 0.04| ## |celltype* | 2| 137| 2.34| 1.07| 2| 2.30| 1.48| 1| 4| 3| 0.28| -1.17| 0.09| ## |time | 3| 137| 121.63| 157.82| 80| 90.33| 87.47| 1| 999| 998| 3.06| 12.33| 13.48| ## |status | 4| 137| 0.93| 0.25| 1| 1.00| 0.00| 0| 1| 1| -3.47| 10.10| 0.02| ## |karno | 5| 137| 58.57| 20.04| 60| 59.37| 29.65| 10| 99| 89| -0.34| -0.83| 1.71| ## |diagtime | 6| 137| 8.77| 10.61| 5| 6.73| 4.45| 1| 87| 86| 4.06| 23.00| 0.91| ## |age | 7| 137| 58.31| 10.54| 62| 59.13| 8.90| 34| 81| 47| -0.63| -0.49| 0.90| ## |prior | 8| 137| 2.92| 4.56| 0| 2.43| 0.00| 0| 10| 10| 0.91| -1.19| 0.39| ## 描述离散变量 library (gmodels)## Cell Contents ## |-------------------------| ## | N | ## | N / Table Total | ## |-------------------------| ## Total Observations in Table: 137 ## ## | squamous | smallcell | adeno | large | ## |-----------|-----------|-----------|-----------| ## | 35 | 48 | 27 | 27 | ## | 0.255 | 0.350 | 0.197 | 0.197 | ## |-----------|-----------|-----------|-----------|
根据
CrossTable
函数的结果可知，样本中有 35 位（25.5%）受试者患有鳞状细胞肺癌（squamous），48 位（35.0%）患有小细胞肺癌（smallcell），27 位（19.7%）患有腺癌（adeno），以及 27 位（19.7%）患有大细胞肺癌（large）。
我们可以使用
gtsummary
包和
tbl_summary
函数来创建一个人口统计表：
## 创建人口统计表 library (gtsummary)"no" ) %>%2
)) %>% "**Variable**" ) %>% "**Table 4. Patient characteristics from the Veterans Lung Cancer Trial**" ) %>%"stat_1" , "stat_2" ) ~ "**Treatment Received**" )
3.2 Kaplan-Meier 生存曲线 & 对数秩检验
接下来，我们为实验组和对照组绘制 KM 曲线，并进行对数秩检验。我们首先使用基本的绘图函数，此时不包括任何颜色或 95% 置信区间，以及删失数据的标记：
kmcurve plot(kmcurve)
图 2 案例：基本 KM 曲线
我们还可以安装
survminer
包来使用
ggsurvplot
函数为 KM 曲线增添额外信息，如 95% 置信区间、对数秩检验结果和风险表。对于对数秩检验，既可以通过在
ggsurvplot
函数中设置
pval = TRUE
选项来直接实现，也可以独立运用
survdiff
函数来进行：
ggsurvplot(kmcurve, TRUE , ## 置信区间 TRUE , ## 对数秩检验 TRUE , ## 风险表 "treatment 1" , "treatment 2" ), ## 图例 "Kaplan-Meier curves" ) ## 表头
图 3 案例：KM 曲线
survdiff(Surv(time, status) ~ trt, data = veteran.data)## Call: ## survdiff(formula = Surv(time, status) ~ trt, data = veteran.data) ## ## N Observed Expected (O-E)^2/E (O-E)^2/V ## trt=1 69 64 64.5 0.00388 0.00823 ## trt=2 68 64 63.5 0.00394 0.00823 ## ## Chisq= 0 on 1 degrees of freedom, p= 0.9
ggsurvplot
和
survdiff
函数在对数秩检验上给出了相同的 p 值，为 0.93 或四舍五入到 0.9。因此，我们不能拒绝零假设，即需要接受对照组（处理组 1）和实验组（处理组 2）之间的生存情况没有统计学上的显著差异的备择假设。
3.3 Cox 比例风险模型
我们可以用 Cox 比例风险函数模型估计对照组（处理组 1）和实验组（处理组 2）之间的风险比：
coxmodel2 exp(coef(coxmodel2))## trt ## 1.017901 ## 2.5 % 97.5 % ## trt 0.7143755 1.450389 ### 创建回归结果 "Cox proportional hazards regression model" ) %>% 'left' )
由此可知，风险比为 1.02 (
)，95% 置信区间为 0.71 (
0.71) 到 1.45 (
)。由于 95% 置信区间包含风险比
