生存分析
2018-09-03
2018-09-03
1 基本概念
生存分析可以用来估计总体的生存率及其有关指标,如癌症患者治疗后的生存时间资料;还可以对不同处理组生存率进行比较,如比较哪种癌症治疗方法更有效果;同时还可探索生存时间资料的影响因素;预测不同的影响因素水平的个体生存时间和一定时间的生存率。近年来也越来越多人用在互联网数据挖掘中,例如去预测信息在社交网络的传播程度,或者去预测用户流失的概率。
1. 生存分析
是研究生存现象和响应时间数据及其统计规律的一门学科。将事件的结果(终点事件)和出现这一结果所经历的时间结合起来分析的一种统计分析方法。生存分析与其他多因素分析最大的区别,就是生存分析考虑了观测结果出现的时间长短。
2. 事件
指研究中规定的生存研究的终点,在研究开始之前就已经制定好。根据研究性质的不同,事件可以是患者的死亡、疾病的复发、仪器的故障,也可以是下岗工人的再就业等等。
3. 生存时间
指从某一起点到事件发生所经过的时间。生存时间不仅仅指医学中的存活,也可以是机器出故障前的正常运行时间等等。
终点事件与起始事件之间的时间间隔。
终点事件:研究者所关系的特定结局。
起始事件:反应研究对象生存过程的起始特征的事件。
终点事件与起始事件是相对而言的,都是由特定的研究目的所决定的,是整个研究过程的标尺,需要在设计时明确规定,并在研究期间严格遵守,不能随意改变。比如研究药品是否起作用和起作用的时间问题,可以将服药作为起始事件,终点事件设为痊愈。
4. 生存时间的类型
广义上指某个起点事件开始到某个终点事件发生所经历的时间,度量单位可以是年、月、日、小时等,常用符号t所示。这个时间也未必是通常意义上的时间,也可以是和时间相关的变量。比如距离等,具体要根据研究目的而定义。根据研究对象的结局,生存时间数据可分为两种类型:
完全数据:从观察起点到发生死亡事件所经历的时间。
不完全数据:生存时间观察过程的截止不是由于死亡事件,而是由其他原因引起的。不完全数据分为:删失数据,截断数据。
不完全主要原因:
- 失访:指失去联系;
- 退出:死于非研究因素或非处理因素而退出研究;
- 终止:设计时规定的时间已到而终止观察,但研究对象仍然存活。
删失的表现形式:
- 右删失:只知道实际寿命大于某数;
- 左删失:只知道实际寿命小于某数;
- 区间删失:只知道实际寿命在一个时间区间内。
截尾数据的产生:
往往是因为实验设计的要求使得数据天然具有上界或者下界。
- 左截尾:数据都大于某个值。
- 右截尾:数据都小于某个值。
如一个实验研究退休职工的生存情况,那么显然这些数据都是左截尾的,因为所有个体的年龄都大于退休年龄(如t≥60)
5. 条件生存概率
生存时间的分布一般不呈现正态分布。条件生存概率表示某时段开始存活的个体,到该时段结束时仍存活的可能性。
6. 生存率
生存概率又称为生存率或生存函数,表示观察对象经历t个单位时间段后仍存活的可能性。
生存率与条件生存概率不同。条件生存概率是单个时段的结果,而生存率是条件生存概率的累积,是多个时段的累积结果。例如,3年生存率是第1年存活,第2年也存活,第3年还存活的可能性。
7. 生存曲线
以观察时间为横轴,以生存率为纵轴,将各个时间点所对应的生存率连接在一起的曲线图。
生存曲线是一条下降的曲线,分析时应注意曲线的高度和斜率。平缓的生存曲线表示高生存率或较长生存期,陡峭的生存曲线表示低生存率或较短生存期。
8. 中位生存期
又称半数生存期,表示恰好有50%的个体尚存活的时间。中位生存期越长,表示疾病的预后越好。
9.生存函数S(t)
观察对象的生存时间T大于某时刻t的概率。
\[S(t)=\frac{时刻t的存活例数}{期初观察例数}\]
10. 死亡函数F(t)
观察对象的生存时间T不大于某时刻t的概率。
\[F(t)=1-S(t)\]
11. 死亡密度函数f(t)
观察对象在某时刻t的瞬时死亡率。
\[f(t)=\frac{\nabla{F(t)}}{\nabla{t}}\]
12. 风险函数h(t)
生存到时刻t的观察对象在时刻t的瞬时死亡率。 \[h(t)=\frac{f(t)}{S(t)}\]
2 生存分析的方法
1.描述法
根据样本观测值直接计算每一个时间点或时间区间上的生存函数、死亡函数、风险函数等,并采用列表或绘图的形式显示生存时间的分布规律。方法对数据分布无要求。
2.非参数法
生存时间的分布未知,估计生存函数;或检验危险因素对生存时间的影响。常用乘积极限法和寿命表法,可以用来估计生存函数,比较两组或多组生存分布函数。分析危险因素对生存时间的影响。
3.参数法
根据样本观测值来估计假设分布模型中的参数,获得生存时间的概率分布模型。生存时间经常服从的分布有:指数分布、Weibull分布、对数正态分布、对数Logistic分布、Gamma分布。不仅能处理非参数法解决的问题,同时还可建立生存时间与危险因素之间的关系模型。但是我们必须提前知道生存时间的分布。
4.半参数法
不需要对生存时间的分布做出假定,但是却可以通过一个模型来分析生存时间的分布规律,以及危险因素对生存时间的影响,最著名的就是COX回归。不仅能处理参数法解决的问题,还不需要事先知道生存时间的分布。
3 R实现
R实现生存分析涉及survival包和OIsurv包,同时数据集来自于KMsurv包。
# 在第一次导入包之前记得安装相关的包
library(KMsurv)
library(survival)
library(OIsurv)
# 查看数据集
data(aids)
head(aids)##   infect induct adult
## 1   0.00   5.00     1
## 2   0.25   6.75     1
## 3   0.75   5.00     1
## 4   0.75   5.00     1
## 5   0.75   7.25     1
## 6   1.00   4.25     13.1 生存类型数据
survival包中的许多函数将方法应用到Surv对象,这些对象是使用Surv()函数创建的生存类型对象。在这里,我们讨论了右删失Surv对象的构造和左截尾右删失Surv对象的构造。
具体形式:
Surv(time, event), Surv(time, time2, event, type)
1. 右截尾数据
在Surv()函数中仅需要两个参数:一个时间向量和一个表示哪些时间被观察和删失的向量。
# 查看数据集
data(tongue)
head(tongue)##   type time delta
## 1    1    1     1
## 2    1    3     1
## 3    1    3     1
## 4    1    4     1
## 5    1   10     1
## 6    1   13     1attach(tongue)
my.surv.object <- Surv(time[type==1], delta[type==1])
my.surv.object##  [1]   1    3    3    4   10   13   13   16   16   24   26   27   28   30 
## [15]  30   32   41   51   65   67   70   72   73   77   91   93   96  100 
## [29] 104  157  167   61+  74+  79+  80+  81+  87+  87+  88+  89+  93+  97+
## [43] 101+ 104+ 108+ 109+ 120+ 131+ 150+ 231+ 240+ 400+detach(tongue)标识加号是指被右截的观察结果。Surv()中的第一个参数应该输入观察时间和右截尾时间的向量。第二个参数使用指标向量来表示事件是否被观察。
2. 左截尾右删失数据
第一个参数输入左截时间,第二个参数输入生存事件向量和删失时间,第三个参数输入右删失指示向量。
# 查看数据集
data(psych)
head(psych)##   sex age time death
## 1   2  51    1     1
## 2   2  58    1     1
## 3   2  55    2     1
## 4   2  28   22     1
## 5   1  21   30     0
## 6   1  19   28     1attach(psych)
my.surv.object <- Surv(age, age+time, death)
my.surv.object##  [1] (51,52]  (58,59]  (55,57]  (28,50]  (21,51+] (19,47]  (25,57] 
##  [8] (48,59]  (47,61]  (25,61+] (31,62+] (24,57+] (25,58+] (30,67+]
## [15] (33,68+] (36,61]  (30,61+] (41,63]  (43,69]  (45,69]  (35,70+]
## [22] (29,63+] (35,65+] (32,67]  (36,76]  (32,71+]detach(psych)3.2 Kaplan-Meier估计
生存分析的主要目的是估计生存函数,常用的方法有Kaplan-Meier法。Kaplan-Meier估计是生存函数S(t)的非参数极大似然估计(MLE)。该估计是一个在观测事件时间\(t_i\)上跳跃的阶跃函数。
具体形式:
survfit(formula, conf.int = 0.95, conf.type = “log”)
# 查看数据集
data(tongue)
attach(tongue)
my.surv <- Surv(time[type==1], delta[type==1])
survfit(my.surv~1)## Call: survfit(formula = my.surv ~ 1)
## 
##       n  events  median 0.95LCL 0.95UCL 
##      52      31      93      67      NASurvfit()也有一些可选参数。例如,置信水平使用第二个参数conf.int更改(例如,90%的置信范围,conf.int=0.90)。conf.type参数描述了置信区间的类型。默认的是“log”,它等于对转换函数g(t) =log(t)。“log-log”选项使用g(t)=log(−log(t))。一个线性置信区间是使用参数conf.type=“plain”创建的。与R中的许多函数一样,survfit()函数返回可以通过正确的命令访问隐藏信息。下面我们考虑这个隐藏信息的几个元素,它们存储在一个列表中。若需要查看整个信息列表的结构,可对my.fit和summary(my.fit)应用str函数。
# 查看数据集
my.fit <- survfit(my.surv~1)
summary(my.fit)$surv# 生存函数值##  [1] 0.9807692 0.9423077 0.9230769 0.9038462 0.8653846 0.8269231 0.8076923
##  [8] 0.7884615 0.7692308 0.7500000 0.7115385 0.6923077 0.6730769 0.6538462
## [15] 0.6340326 0.6142191 0.5944056 0.5745921 0.5547786 0.5342312 0.5061138
## [22] 0.4779963 0.4481216 0.4161129 0.3814368 0.3051494 0.2288621summary(my.fit)$time##  [1]   1   3   4  10  13  16  24  26  27  28  30  32  41  51  65  67  70
## [18]  72  73  77  91  93  96 100 104 157 167summary(my.fit)$n.risk##  [1] 52 51 49 48 47 45 43 42 41 40 39 37 36 35 33 32 31 30 29 27 19 18 16
## [24] 14 12  5  4summary(my.fit)$n.event##  [1] 1 2 1 1 2 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1summary(my.fit)$sstd.err## NULLsummary(my.fit)$lower##  [1] 0.9441432 0.8810191 0.8534195 0.8271685 0.7774159 0.7302341 0.7073724
##  [8] 0.6849189 0.6628317 0.6410776 0.5984672 0.5775712 0.5569274 0.5365233
## [15] 0.5156073 0.4949392 0.4745101 0.4543127 0.4343412 0.4137057 0.3837502
## [22] 0.3546085 0.3240114 0.2916674 0.2571797 0.1692469 0.1010963str(my.fit)## List of 13
##  $ n        : int 52
##  $ time     : num [1:45] 1 3 4 10 13 16 24 26 27 28 ...
##  $ n.risk   : num [1:45] 52 51 49 48 47 45 43 42 41 40 ...
##  $ n.event  : num [1:45] 1 2 1 1 2 2 1 1 1 1 ...
##  $ n.censor : num [1:45] 0 0 0 0 0 0 0 0 0 0 ...
##  $ surv     : num [1:45] 0.981 0.942 0.923 0.904 0.865 ...
##  $ type     : chr "right"
##  $ std.err  : num [1:45] 0.0194 0.0343 0.04 0.0452 0.0547 ...
##  $ lower    : num [1:45] 0.944 0.881 0.853 0.827 0.777 ...
##  $ upper    : num [1:45] 1 1 0.998 0.988 0.963 ...
##  $ conf.type: chr "log"
##  $ conf.int : num 0.95
##  $ call     : language survfit(formula = my.surv ~ 1)
##  - attr(*, "class")= chr "survfit"str(summary(my.fit))## List of 15
##  $ n            : int 52
##  $ time         : num [1:27] 1 3 4 10 13 16 24 26 27 28 ...
##  $ n.risk       : num [1:27] 52 51 49 48 47 45 43 42 41 40 ...
##  $ n.event      : num [1:27] 1 2 1 1 2 2 1 1 1 1 ...
##  $ n.censor     : num [1:27] 0 0 0 0 0 0 0 0 0 0 ...
##  $ surv         : num [1:27] 0.981 0.942 0.923 0.904 0.865 ...
##  $ type         : chr "right"
##  $ std.err      : num [1:27] 0.019 0.0323 0.037 0.0409 0.0473 ...
##  $ lower        : num [1:27] 0.944 0.881 0.853 0.827 0.777 ...
##  $ upper        : num [1:27] 1 1 0.998 0.988 0.963 ...
##  $ conf.type    : chr "log"
##  $ conf.int     : num 0.95
##  $ call         : language survfit(formula = my.surv ~ 1)
##  $ table        : Named num [1:9] 52 52 52 31 147 ...
##   ..- attr(*, "names")= chr [1:9] "records" "n.max" "n.start" "events" ...
##  $ rmean.endtime: num 400
##  - attr(*, "class")= chr "summary.survfit"Kaplan-Meier估计可以使用plot(my.fit)绘制。
plot(my.fit,main="Kaplan-Meier estimate with 95% confidence bounds", xlab = "time",ylab="survival function")上面分析不同的组被包含在一个单一的Surv对象中。例如,tongue数据集中的type变量描述了患者的DNA谱。我们可以通过对type变量上的Surv对象进行回归得到这些组的Kaplan-Meier估计:
my.fit1 <- survfit( Surv(time, delta) ~ type )
summary(my.fit1)## Call: survfit(formula = Surv(time, delta) ~ type)
## 
##                 type=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1     52       1    0.981  0.0190        0.944        1.000
##     3     51       2    0.942  0.0323        0.881        1.000
##     4     49       1    0.923  0.0370        0.853        0.998
##    10     48       1    0.904  0.0409        0.827        0.988
##    13     47       2    0.865  0.0473        0.777        0.963
##    16     45       2    0.827  0.0525        0.730        0.936
##    24     43       1    0.808  0.0547        0.707        0.922
##    26     42       1    0.788  0.0566        0.685        0.908
##    27     41       1    0.769  0.0584        0.663        0.893
##    28     40       1    0.750  0.0600        0.641        0.877
##    30     39       2    0.712  0.0628        0.598        0.846
##    32     37       1    0.692  0.0640        0.578        0.830
##    41     36       1    0.673  0.0651        0.557        0.813
##    51     35       1    0.654  0.0660        0.537        0.797
##    65     33       1    0.634  0.0669        0.516        0.780
##    67     32       1    0.614  0.0677        0.495        0.762
##    70     31       1    0.594  0.0683        0.475        0.745
##    72     30       1    0.575  0.0689        0.454        0.727
##    73     29       1    0.555  0.0693        0.434        0.709
##    77     27       1    0.534  0.0697        0.414        0.690
##    91     19       1    0.506  0.0715        0.384        0.667
##    93     18       1    0.478  0.0728        0.355        0.644
##    96     16       1    0.448  0.0741        0.324        0.620
##   100     14       1    0.416  0.0754        0.292        0.594
##   104     12       1    0.381  0.0767        0.257        0.566
##   157      5       1    0.305  0.0918        0.169        0.550
##   167      4       1    0.229  0.0954        0.101        0.518
## 
##                 type=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1     28       1   0.9643  0.0351       0.8979        1.000
##     3     27       1   0.9286  0.0487       0.8379        1.000
##     4     26       1   0.8929  0.0585       0.7853        1.000
##     5     25       2   0.8214  0.0724       0.6911        0.976
##     8     23       1   0.7857  0.0775       0.6475        0.953
##    12     21       1   0.7483  0.0824       0.6031        0.929
##    13     20       1   0.7109  0.0863       0.5603        0.902
##    18     19       1   0.6735  0.0895       0.5190        0.874
##    23     18       1   0.6361  0.0921       0.4790        0.845
##    26     17       1   0.5986  0.0939       0.4402        0.814
##    27     16       1   0.5612  0.0952       0.4025        0.783
##    30     15       1   0.5238  0.0959       0.3658        0.750
##    42     14       1   0.4864  0.0961       0.3302        0.716
##    56     13       1   0.4490  0.0957       0.2956        0.682
##    62     12       1   0.4116  0.0948       0.2621        0.646
##    69     10       1   0.3704  0.0938       0.2255        0.608
##   104      8       2   0.2778  0.0904       0.1468        0.526
##   112      5       1   0.2222  0.0877       0.1025        0.482
##   129      4       1   0.1667  0.0815       0.0639        0.435
##   181      2       1   0.0833  0.0717       0.0155        0.449detach(tongue)3.3 累计风险函数
累积风险函数在连续数据上和生存函数的关系如下:
\[S(t)=exp{(-H(t))}\] 风险函数的最大似然估计可以通过kaplan-Meier估计的逆变换得到:
\[\hat{H(t)}=-log\hat{S(t)}\]
尽管survival包中没有函数自动计算任何一个表单,但是通过summary(survfit())返回的对象可以用来计算估计值:
data(tongue)
attach(tongue)
my.surv <- Surv(time[type==1], delta[type==1])
my.fit <- summary(survfit(my.surv~1))
hat.H <- -log(my.fit$surv)
hat.H##  [1] 0.01941809 0.05942342 0.08004271 0.10109612 0.14458123 0.19004360
##  [7] 0.21357410 0.23767165 0.26236426 0.28768207 0.34032581 0.36772478
## [13] 0.39589566 0.42488319 0.45565485 0.48740355 0.52019337 0.55409493
## [19] 0.58918625 0.62692657 0.68099379 0.73815221 0.80269073 0.87679870
## [25] 0.96381008 1.18695363 1.47463570上面得到的结果即是累积风险函数值。下图是累积风险函数值在时间上的走势。
plot(my.fit$time, hat.H, xlab="time", ylab="cumulative hazard value",main="cumulative hazards", type="s")3.4 多组生存时间比较
survdiff(formula, rho=0)
给定两个或两个以上的样本,生存时间是否存在显著差异? 我们可以结合风险函数假设为:
\[H_0:对所有的t有 h_1(t)=h_2(t)=...=h_n(t)\\H_1:h_i(t_0) \neq h_j(t_0)\] 函数survdiff()用于这个假设检验。第一个参数是生存对象对一个分类的协变量,通常是一个变量指定哪个组对应哪个生存时间。survdiff()返回的对象提供了有用的信息。
data(btrial)
attach(btrial)
survdiff(Surv(time, death)~im)## Call:
## survdiff(formula = Surv(time, death) ~ im)
## 
##       N Observed Expected (O-E)^2/E (O-E)^2/V
## im=1 36       16    20.19     0.869      5.49
## im=2  9        8     3.81     4.599      5.49
## 
##  Chisq= 5.5  on 1 degrees of freedom, p= 0.02上面假设检验结果P值等于0.02小于0.05,所以拒绝原假设,两组的样本的生存时间存在显著差异。
第二个survdiff()参数rho,默认等于0,参数rho用来设置权重,如果给生存曲线的第一部分赋予更大的权重,rho大于0。如果更大权重在生存曲线的后面部分,用小于0的。
3.5 Cox比例风险模型
如果考虑其他影响生存时间分布的因素,可以使用Cox回归模型(也叫比例风险模型)它是一种半参数回归模型,利用数学模型拟合生存分布与影响因子之间的关系,评价影响因子对生存函数分布的影响程度。该模型以生存结局和生存时间为应变量,可同时分析众多因素对生存期的影响,能分析带有截尾生存时间的资料,且不要求估计资料的生存分布类型。这里的前提是影响因素的作用不随时间改变,如果不满足这个条件,则应使用含有时间依存协变量的Cox回归模型。
现在我们考虑时间依存的协变量,如干预措施或可能的环境因素在研究过程中产生变化。为了在R中使用时间相关的协变量,我们大量地应用左截断。例如,如果有一个受干预的病人\(i\),然后我们将病人\(i\)分为两个独立的观察:干预前和干预后。假设患者干预发生在\(t_i = 45\)时,观察病人\(i\)在\(t_j = 58\)时发生的事件。我们将这个病人的病历分成两片段:0到45和45到58。干预的协变量可以在每个时间间隔被赋予不同的值。
我们考虑以下示例(模拟数据)。病人记录来自一个诊所。在患者清醒后的150天内获得。数据的事件是酗酒,有四个变量。
* event变量描述被观察的或被删失的复发时间;
* delta变量描述该事件是被观察(真实)还是被删失(假);
* gender变量是一个独立于时间的协变量;
* int是一个依赖时间的协变量表示患者是否进行了干预,有些患者在计划的干预时间前复发,所以里面有NA值。所有干预均发生在患者清醒后。也就是说,只要没有复发,干预协变量随着时间的推移而变化。数据relapse保存在OIsurv包中。
data(relapse)
head(relapse)##   event delta gender inter
## 1   150 FALSE      0    84
## 2    53  TRUE      1    50
## 3    12  TRUE      1    NA
## 4   150 FALSE      0    89
## 5   150 FALSE      1    77
## 6   135  TRUE      1     7这些数据可以通过以下两个步骤进行建模:
1. 构造生存记录,包括左和右删失的观察。每一个接受干预的病人的生存记录分成两个:干预前和干预后。
2. 新的生存记录可以通过coxph()函数运行。
我们首先初始化新的向量来表示生存记录的变量。变量\(t_1\)和 \(t_2\)分别表示开始和结束时间,\(d\)分别表示观察到复发(TRUE)还是右删失(FALSE),\(g\)代表性别,\(i\)代表病人是否正在经历干预治疗。
N <- dim(relapse)[1]
t1 <- rep(0, N+sum(!is.na(relapse$int)))# 初始化开始时间
t2 <- rep(-1, length(t1))# 构建结束时间向量
d <- rep(-1, length(t1))# 时间是否删失
g <- rep(-1, length(t1))# 性别
i <- rep(FALSE, length(t1))# 初始化干预为FALSE接下来,对每个病人进行检查,看他们是否接受过干预。如果他们没有,那么他们的记录被简单地复制到所使用的新变量中。如果有干预,他们的观察分为两部分:干预前和干预后。
j <- 1
for(ii in 1:dim(relapse)[1]){ 
  if(is.na(relapse$int[ii])){
   t2[j] <- relapse$event[ii]
   d[j] <- relapse$delta[ii]
   g[j] <- relapse$gender[ii]
   j <- j+1
  } else {          # 若发生干预事件,分裂数据
    g[j+0:1] <- relapse$gender[ii] # 每个时间的性别是一样的
    d[j] <- 0 # 干预前没有复发
    d[j+1] <- relapse$delta[ii] # 干预后是否复发
    i[j+1] <- TRUE # 干预事件发生
    t2[j] <- relapse$int[ii]-1 # 干预前的结束时间
    t1[j+1] <- relapse$int[ii]-1 # 干预后的开始时间
    t2[j+1] <- relapse$event[ii] # 干预后的结束时间
    j <- j+2 # 两条记录增加
  }
}虽然许多患者有随时间变化的协变量(即干预),但每一个新的间隔存在不改变的协变量,生存对象就可以建立Cox PH模型。
mySurv <- Surv(t1, t2, d)# 构建生存对象
myCPH <- coxph(mySurv ~ g + i)# cox回归
test.ph <- cox.zph(myCPH)
test.ph##            rho chisq     p
## g      -0.0325 0.216 0.642
## iTRUE  -0.0225 0.109 0.741
## GLOBAL      NA 0.314 0.855使用R语言survival包中的cox.zph函数进行模型诊断,从上面的结果可以看出,两个变量的P值都大于0.05,说明每个变量均满足PH检验,而模型的整体检验P值0.855,模型整体满足PH检验。
这个例子是一个简单的例子;只有一个与时间相关的协变量,它每次最多改变一次。在某些情况下,可能有许多随时间变化的协变量。在大多数情况下,可以使用相同的方法。每当一个协变量从一个时间单位变化到下一个时间单位,将间隔分割为两个,根据需要使用左截断和右删失。
4 本章汇总
| 名称 | 类别 | 功能 | 
|---|---|---|
| KMsurv | 包 | 生存分析的数据集 | 
| survival | 包 | 生存分析包 | 
| OIsurv | 包 | 生存分析包 | 
| data() | 函数 | 载入R内置数据集 | 
| head() | 函数 | 查看数据 | 
| attach() | 函数 | 减少重复键入表名 | 
| detach() | 函数 | 解除attach函数 | 
| Surv() | 函数 | 生存对象产生 | 
| survfit() | 函数 | 生存函数拟合 | 
| survdiff() | 函数 | 生存时间比较 | 
| dim(relapse) | 函数 | 得到对象的大小 | 
| rep() | 函数 | 将对象重复多少次产生新对象 | 
| is.na() | 函数 | 判断是否为空值 | 
| coxph() | 函数 | cox回归拟合函数 | 
| cox.zph() | 函数 | cox回归模型检验 | 
5 参考文献
[1] 彭非. 生存分析[M]. 中国人民大学出版社, 2004.
[2] http://www.openintro.org/stat/down/Survival-Analysis-in-R.pdf
