R语言系列第五期：④R语言与生存分析
对于寿命数据的分析，在生物学和医药学中是非常重要的话题。除此之外，在工程应用中的可靠性分析中也非常重要。寿命数据往往是高度非正态数据，因此使用标准的线性模型可能会有很多问题。
而相对于逻辑回归的只有分类结局，只考虑终点事件是否出现的情况。
生存分析的结局还会考虑观察对象达到终点所经历的时间长短。生存分析是将终点时间的出现与否和达到终点所经历的时间结合起来分析的一类统计分析方法。
我们使用包survival，作者是Terry Therneau，这个包继承了一批生存分析的先进工具。我们这里只会为大家介绍其中的一小部分。
首先，载入survival包：
> library(survival)
#Tips：同样，如果没有安装这个包的，请先安装
> install.packages(“survival”)
Survival包中最常用的是“Surv”类的对象，它是时间和状态信息合并在一起的一种数据结构，这种对象由函数Surv()生成，该函数带有两个参数，其一是观测到的时间，其二是事件状态标志。后者被编码为逻辑变量，0/1或者1/2。当然还是前一个编码方式比较推荐。
#Tips：其实，Surv()函数还有3个参数，用来处理开始时间、结束时间以及时间区间内截断事件的数据。
我们使用K.T.Drzewiecki收集的melanom（黑色素瘤）数据集，数据可以通过以下方法获取：
> library(ISwR)
载入程辑包：‘ISwR’
The following object is masked from‘package:survival’:
lung
#Tips：这里有个无关紧要的提醒，当先挂载ISwR,然后挂载survival也会出现关于lung数据集的提醒。
> names(melanom)
[1] "no" "status" "days" "ulc" "thick" "sex"
#Tips：names()函数可以告诉我们某个数据集包含的变量名称。变量status是病人在研究期结束时的状态：1表示“死于恶性黑色素瘤”，2表示“1978年1月1日的时候还是存活的”，3表示“死于其他原因”。变量days是观测的生存日期，ulc指是否存在溃疡性肿瘤（1表示是，2表示否），thick是以1/100mm计量的厚度，sex指病人性别（1表示女性，2表示男性）。
然后我们把melanom放在检索路径上：
> attach(melanom)
我们希望创建一个Surv对象，其中变量status的值2和3作为删失。由以下代码实现：
> Surv(days,status==1)
[1] 10+ 30+ 35+ 99+ 185 204 210 232 232+
[10] 279 295 355+ 386 426 469 493+ 529 621
[19] 629 659 667 718 752 779 793 817 826+
[28] 833 858 869 872 967 977 982 1041 1055
......
[181] 3476+ 3523+ 3667+ 3695+ 3695+ 3776+ 3776+ 3830+ 3856+
[190] 3872+ 3909+ 3968+ 4001+ 4103+ 4119+ 4124+ 4207+ 4310+
[199] 4390+ 4479+ 4492+ 4668+ 4688+ 4926+ 5565+
#Tips：上述代码直接打印出了Surv对象，其中‘+’号标记出删失的观测。比如10+就表示这个病人并未在10天内死于黑色素瘤，然而无法继续进行跟踪试验（实际上，这个患者是死于其他原因），185表示这个病人在术后半年多时间就死于黑色素瘤了。Surv()函数的第二个参数是一个逻辑向量：status==1对于死于黑色素瘤的患者观测为TRUE，其他为FALSE。
A. Kaplan—Meier估计
Kaplan-Meier估计（乘积极限法）用以计算右侧截断数据的生存函数的估计，这个估计是一个阶梯函数，它的跳跃点是给定的时间点。
生存函数的Kaplan-Meier估计的计算可以通过调用函数survfit()来实现。该函数最简单的形式只带有一个参数，即为Surv对象。函数返回一个survfit对象。代码如下：
> surv.all<-survfit(Surv(days,status==1)~1)
> surv.all
Call: survfit(formula = Surv(days, status == 1) ~ 1)
n events median 0.95LCL 0.95UCL
205 57 NA NA NA
#Tips：可以看到单纯使用survfit()函数并没有提供多少信息，你获得的信息包括一些汇总的统计量，以及对中位生存中位数的一个估计。要看真正的估计，我们需要对survfit对象使用summary()函数。
> summary(surv.all)
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
...
3042 52 1 0.664 0.03994 0.590 0.747
3338 35 1 0.645 0.04307 0.566 0.735
通常，相比于数值，你可能对Kaplan-Meier估计的图像更感兴趣。你只需要输入一下代码，就可以做出图形：
> plot(surv.all)
#Tips：曲线上的记录表示截断时间，两侧的虚线围成的就是置信区间。
将多条生存曲线同时画在一个图上有时候更有用，这样有助于对其进行直接比较。我们要获取不同性别的生存曲线，可以输入如下代码：
> surv.bysex<-survfit(Surv(days,status==1)~sex)
> plot(surv.bysex)
#Tips：原来的1被改成sex变量了，一开始是想所有数据凑在一起，所以写1，后来我们要按照sex分组，所有“~”后要放上sex变量。另外，我们可以注意到置信带没有了，因为默认做分组的时候会关闭，只要我们改变参数就可以，另外，当两个置信区间带都画出来可能就容易混，所以最好更改颜色。
> plot(surv.bysex,conf.int=T,col=c("black","grey"))
#Tips：同样，你可以在不分组的时候设置conf.int=F来避免画置信区间，如果希望置信度为99%，可以设置conf.int=0.99.
B. 对数秩检验
对数秩检验可以检验两条或者多条生存曲线是否相同，是典型的非参数检验。对数秩检验的计算是通过survdiff()来计算得到的：
> 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
#Tips：对模型的指定和线性模型或者广义线性模型中类似，但是这个检验只能处理分组变量，所以如果你在方程右边指定了多个变量，则检验是对由这些变量所有取值组合形成的分组进行的。该检验对于因子型和数值型的编码方式不做区分，对于函数survfit()也是一样。
指定分层分析也可以，比如，可以按照是否是溃疡型（ulc）分层，然后按照性别分组的对数秩检验：
> 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
#Tips：strata()就是一个把变量做成分层变量的函数，另外，这样做就使得性别的效应不明显了，可能是因为，相比于女性，男性在疾病发展到比较严重的阶段才会寻求医疗救助，所以对病情的不同阶段分别进行假设检验的时候，性别之间的差别的显著性就下降了。
C. Cox比例风险模型
比例风险模型允许用类似lm或者glm的回归模型来分析数据，并且假设在对数风险这一刻度上，关系是线性的。模型可以通过用极大似然Cox函数拟合得到。作为例子，我们仅仅考虑包含一个回归变量sex的模型：
> 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 )
Rsquare= 0.03 (max possible= 0.937 )
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
#Tips：其中coef是估计得到的两个组之间的风险比的对数—ln(RR)，因此真正的风险比是exp(coef)，其实得到的就是RR（相对危险度）。最后三行是三种检验结果，在大样本的时候，这三种检验结果是一样，但是对于小样本也许会出现差别。
同样分层分析：
>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 )
Rsquare= 0.063 (max possible= 0.9 )
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
Tips：可以看到sex的显著性被大大削减了。
Cox模型假设一个潜在的基线模型对应一条生存曲线。在分层分析中，每一个层中都会有一条如此的曲线。可以通过在coxph的输出中用survfit()函数得到该曲线：
> plot(survfit(coxph(Surv(days,status==1)~sex+log(thick)+strata(ulc))))
