12-Survival analysis

《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 参数来指定要为其计算生存曲线的数据框。