12-Survival analysis
2020-07-01 本文已影响0人
共由小兑
《Introductory Statistics With R》阅读笔记 ☀️☀️☀️☀️
寿命分析是生物学和医学领域的重要课题，尤其是在工程应用中的可靠性分析中。这样的数据通常高度非正态分布，因此使用标准线性模型是有问题的。生命周期数据经常被审查（censored），因为不知道确切的生命周期，只知道它比给定的值长。例如，在一项癌症试验中，有些人无法进行随访，或者干脆活过了研究期。在统计分析中忽略审查是错误的，有时会造成极端的后果。
12.1 Essential concepts
- X :true lifetime
-
T :censoring time
T可以是一个随机变量，也可以是一个固定时间，这取决于上下文，但如果它是随机的，那么它通常应该是非信息的，我们在这里描述的方法是适用的。
有时，“其他原因造成的死亡”被认为是给定疾病死亡率的censoring event，在这些情况下，确保这些其他原因与疾病状态无关尤为重要。 -
S(t):survival function
测量在给定时间存活的概率。 实际上，它只是1减去X的累积分布函数，1- F(t)。 -
h(t):hazard function or force of mortality
假设被试在t时刻还活着，测量在短时间t内死亡的(无穷小的)风险。 -
如果寿命分布具有密度f，则h(t) = f(t)/S(t)
这通常被认为是比(例如)生存分布的平均中值更基本的数量，并被用作建模的基础。
12.2 Survival objects
survival
这个包实现了大量的先进技术。对于目前的目的，只使用它的一个小子集。
> library(survival) #载入survival
- 生存与Surv类的对象一起工作，该类是结合时间和检查信息的数据结构。
- 类对象使用Surv函数构造，带两个参数：observation time和event indicator。后者可以编码为逻辑变量，0/1变量或1/2变量。实际上，Surv还可以与三个参数一起使用，用于处理具有开始时间和结束时间的数据(“交错输入”)和间隔审查数据(知道事件发生在两个日期之间，例如在疾病的重复测试中)。
> library(ISwR)
> data(melanom) #载入生存数据
> attach(melanom) #把数据集放在检索路径上
> names(melanom)
[1] "no" "status" "days" "ulc" "thick" "sex"
> head(melanom)
no status days ulc thick sex
1 789 3 10 1 676 2
2 13 3 30 2 65 2
3 97 2 35 2 134 2
- status变量是研究结束时患者状态的一个指标:1表示“死于恶性黑色素瘤”，2表示“在1978年1月1日还活着”，3表示“死于其他原因”；
- days为观察时间，以天为单位；
- ulc为(1/2)有/无肿瘤有无溃疡；
- thick为厚度,单位1/100 mm；
- sex为患者性别，1表示女性，2表示男性。
#创建一个Surv对象，将状态变量的值2和3视为检查
> Surv(days, status==1)
[1] 10+ 30+ 35+ 99+ 185 204 210 232 232+ 279 295 355+
[13] 386 426 469 493+ 529 621 629 659 667 718 752 779
- 带有“ +”标记检查观察结果。
- 注意：带有“+”的是状态1的，剩下的是状态2和3的。
12.3 Kaplan–Meier estimates
Kaplan-Meier estimator 允许在右截尾情况下计算估计的生存函数。它也被称为product-limit estimator（乘积极限估计），因为描述这个过程的一种方法是，它将没有截尾观察或没有死亡的区间内的条件生存曲线相乘。
使用称为 survfit 的函数来计算生存函数的 Kaplan-Meier 估计量。它只需要一个参数，即Surv对象。 它返回一个survfit对象。 如上所述，我们认为“其他原因导致的死亡”是一种检查，并执行以下操作：
> survfit(Surv(days,status == 1)~1)
Call: survfit(formula = Surv(days, status == 1) ~ 1)
n events median 0.95LCL 0.95UCL
205 57 NA NA NA
- 得到一些简要的统计数据和中值生存的估计，在这种情况下，后者甚至都不有趣，因为估计是无限的。
要查看实际的 Kaplan-Meier 估计，请在 survfit 对象上使用 summary。将 survfit 对象保存到一个名为 surv.all 变量中。这是因为它包含了所有病人的原始生存功能，而不考虑病人的特征。
> surv.all <- survfit(Surv(days,status==1)~1)
> summary(surv.all,censored=F)
Call: survfit(formula = Surv(days, status == 1) ~ 1)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
185 201 1 0.995 0.00496 0.985 1.000
204 200 1 0.990 0.00700 0.976 1.000
210 199 1 0.985 0.00855 0.968 1.000
- 它包含生存函数在事件时刻的值。censored=T 可以显示审查时间。
- KM是一个阶跃函数，它的跳跃点是在 time 上给出的，跳跃后的值是在 survival 上给出的。另外，给出了曲线的标准误差估计和真曲线的点态置信区间。
> plot(surv.all)
Figure 12.1. Kaplan–Meier plot for melanoma data (all observations).
在同一图上绘制两个或多个生存函数通常很有用，以便可以直接比较它们（图12.2）。 要获得按性别划分的生存功能，请执行以下操作：
> plot(surv.bysex, conf.int=T, col=c("black","red"))
Figure 12.2. Kaplan–Meier plots for melanoma data, grouped by gender.
- conf.int=T 打开置信区间，如果想要达到99% 的置信水平，设置 conf.int=0.99。
- col 给每组赋予不同的颜色。
12.4 The log-rank test
对数秩检验用于检验两条或多条生存曲线是否相同。它是基于观察每个死亡时间的人口，并根据每一组中处于危险中的个人的数量，计算预期死亡人数。然后将其与所有死亡时间相加，并通过与χ2检验相似（但不相同）的程序与观察到的死亡人数进行比较。
利用函数 survdiff 计算对数秩检验。这实际上实现了一系列由参数 rho 指定的测试，允许对零假设进行各种非比例的风险替代，但是 rho = 0 的默认值给出了 log-rank 测试。
> survdiff(Surv(days,status==1)~sex)
Call:
survdiff(formula = Surv(days, status == 1) ~ sex)
N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 126 28 37.1 2.25 6.47
sex=2 79 29 19.9 4.21 6.47
Chisq= 6.5 on 1 degrees of freedom, p= 0.01
- 该规范使用了线性和广义线性模型的模型公式。 然而，这个测试只能处理分组的数据，所以如果你在右边指定多个变量，它将处理所有预测变量组合产生的数据分组。 它也不区分因子和数字代码。 survfit 也是如此。
也可以指定分层分析，即在数据集的分层中分别进行观测值和预测值的计算。 例如，可以计算按溃疡分层的性别效应的对数等级检验如下:
> survdiff(Surv(days,status==1)~sex+strata(ulc))
Call:
survdiff(formula = Surv(days, status == 1) ~ sex + strata(ulc))
N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 126 28 34.7 1.28 3.31
sex=2 79 29 22.3 1.99 3.31
Chisq= 3.3 on 1 degrees of freedom, p= 0.07
- 请注意，这会使性的影响显得不那么显著。 一种可能的解释可能是，男性在疾病处于比女性更严重的状态时寻求治疗，因此当根据疾病进展的衡量标准进行调整时，性别差异会缩小。
12.5 The Cox proportional hazards model
作为第一个例子，考虑一个单回归变量的模型:
> summary(coxph(Surv(days,status==1)~sex))
Call:
coxph(formula = Surv(days, status == 1) ~ sex)
n= 205, number of events= 57
coef exp(coef) se(coef) z Pr(>|z|)
sex 0.6622 1.9390 0.2651 2.498 0.0125 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
sex 1.939 0.5157 1.153 3.26
Concordance= 0.59 (se = 0.034 )
Likelihood ratio test= 6.15 on 1 df, p=0.01
Wald test = 6.24 on 1 df, p=0.01
Score (logrank) test = 6.47 on 1 df, p=0.01
- coef 指两组之间危险比的估计对数,也将其作为实际危险比 exp(coef) 给出。
- 下面的一行还给出了危险比的倒置比率(交换组)和置信区间。 最后，对模型中的显著效应进行了三次全面的检验。
- 这些在大样本中都是等价的，但在小样本中可能有所不同。 注意，Wald 检验与基于估计系数除以标准误差的 z 检验相同，而得分检验等同于对数秩检验(只要模型只涉及一个简单的分组)。
更详细的例子，涉及连续协变量和分层变量：
> summary(coxph(Surv(days,status==1)~sex+log(thick)+strata(ulc)))
Call:
coxph(formula = Surv(days, status == 1) ~ sex + log(thick) +
strata(ulc))
n= 205, number of events= 57
coef exp(coef) se(coef) z Pr(>|z|)
sex 0.3600 1.4333 0.2702 1.332 0.1828
log(thick) 0.5599 1.7505 0.1784 3.139 0.0017 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
sex 1.433 0.6977 0.844 2.434
log(thick) 1.750 0.5713 1.234 2.483
Concordance= 0.673 (se = 0.039 )
Likelihood ratio test= 13.3 on 2 df, p=0.001
Wald test = 12.88 on 2 df, p=0.002
Score (logrank) test = 12.98 on 2 df, p=0.002
- 可以看出，性别变量的重要性已经进一步降低。
Cox模型假设了一个潜在的基线危险函数，以及一个对应的生存曲线。在分层分析中，每一层都有一条相同的曲线。它们可以通过使用 coxph 输出的 survfit 来提取，当然也可以使用针对 survfit 对象的 plot 方法来绘制(图12.3):
> plot(survfit(coxph(Surv(days,status==1)~log(thick)+sex+strata(ulc))))
Figure 12.3. Baseline survival curves (ulcerated and nonulcerated tumors) in
stratified Cox regression.
- 注意，survfit 的默认设置是为协变量在其平均值处的假个体生成曲线。 本例肿瘤厚度为1.86 mm，性别为1.39(!) 。
- 可以使用 survfit 的 newdata 参数来指定要为其计算生存曲线的数据框。
