生存分析笔记

生存分析

使用对象：生存分析（survival analysis）使是将事件的结果和出现这一结果所经历的时间结合起来分析的一类统计分析方法，不仅考虑事件是否出现，而且也考虑事件出现的时间长短，因此这类方法也称为事件时间分析（time-to-event analysis）。

生存分析的资料通常有以下特点：①同时考虑生存时间和生存结局；②通常含有删失数据；③生存时间的分布通常不服从正态分析。

以案例说明各种术语

以教材中的案例一下：
例19-1 为了估计HIV（人类免疫缺陷病毒）阳性患者的生存时间，某研究者进行了临床随访研究。研究对象是于2002年1月1日至2004年12月31日期间在某市确诊为HIV阳性者，随访这些对象直至死于AIDS（艾滋病）或其并发症（status=1为死亡，0为删失），研究截至日期为2008年12月31日。并记录每个研究对象的性别（sex=1为男， 0为女）、年龄（age，岁）、是否用药（drug=1为用药，0为不用），2017-09-26 15:09:14结果见表19-1所示。

起始事件、终点事件、生存时间

事件可以分为起始事件（initial event）和终点事件（terminal event），从起始事件到终点事件之间所经历的时间跨度为生存时间（survival time），常记为随机㸄T，其中T>=0，其值记为t。

在例19-1中的起始事件就是“确认HIV明性”，感兴趣的终点事件为“死于AIDS或其并发症”。表19-1中的起始日期（entdate）至终点日期（enddate）之间的月数就是生存时间（time），如编号ID个体的生存时间t=10个月。

生存时间的分布

生存时间的分布通常不呈正态分布，而呈偏态分布，如呈指数分析、Weibul分布、Gompertz分布、对数Logistic分布等，下图就是19-1中的生存间分布图，其中组距为6个月，横轴是生存时间，纵轴是频率（%，其分子为频数，分母为总例数）：

删失

生存结局（status）分为“死亡”和“删失”两类，“死亡”是感兴趣的终点事件，其它的终点事件或生存结局都归于删失（censoring，也称为示截尾或终检）。产生删失的原因可能为：①研究截至日期时，感兴趣终点事件仍未出现；②因搬迁等原因失去联系，不知感兴趣终点事件何时发生或是否会发生；③困其它各种原因中途退出；④死亡其它“事件”，例如交通事故或其它产生。

删失分为左删失、区间删失、右删失3类。

左删失：只知道感兴趣终点事件会在目前知晓时间（如截至时间、失访时间、死于其它疾病时间）之前发生。

区间删失：只知道感兴趣终点事件会在某一区间内发生。

右删失：只知道感兴趣终点事件会在知晓昱是之后发生。例如表19-1的ID=2个体的生存时间为1个月，因为是删失个体，所以该个体至少可以存活1个月，记为t=1+(月)，采用“+”就表示个体为右删失。

生存时间图示

每一个体进入随访队列的时间可以相同也可以不相同，当起时间不相同的时候，一般应规定一个时间范围，如例19-1就规定2002年1月1日至2004年12月31日。通常研究也会规研究截至时间，例19-1就规定为2008年12月31日。为了采用图形反映生存时间，特从表19-1中取出了个5个个体。生存结局根据记录的“终点事件”获得。生存时间通地客经2列和第3列的起始日期与终点日期计算获得。

图19-2分别采用年历时间（图19-2a）和生存时间（图19-2b）反映了表19-2中5个个体的生存情况，用“X”表示“死亡”，用“o”表示“删失”。第4、29、78号个体“死于AIDS”，即死于感兴趣的终点事件，其终点标记为“X”。第38号、39号个体为删失数据，其终点标记为“o”。第38号个体删失的原因为“失访”；第89号个体删失的原因为已到“截至日期”，感兴趣终点事件仍然没有发生。

生存率

如果对图19-1的生存时间直方图拟合一条光滑曲线（图19-3），如下所示：
则这条曲线就是死亡概率密度函数（probability density funciotn,PDF）曲线，死亡概率密度函数记为f(t)，其定义为：

$f(t)=\lim_{\Delta \to 0} \frac{P(t\leq T< t + \Delta t)}{\Delta t}$

如果资料中没有删失数据，则生存率如下所示：

$S(t)=P(Tt)=\frac{t时刻仍存活的例数}{观察总例数}$

如果含有删失数据，则要分时段计算生存概率。假定观察对象在各个时段的生存事件独立，应用概率乘法定量将各时段的生存概率相乘得到生存率，如下所示：

$S(t)=p_{1}\cdot p_{2} \cdots p_{i} \cdots p_{k}=S(t_{k-1})\cdot p_{k}$

其中

$p_{i}(i=1,2, \cdots,k)$$是各时段的生存概率。对于以区间陈列的生存时间资料，可以采用以下的公式来计算死亡概率密度函数值： $$f(t)=[S(t_{i-1})-S(t_{i})]$

其中，$wi$为第i区间的时间宽度，$S(t_{i-1})$、$S(t_{i})$分别是时间$t_{i-1})$和$t_{i})$的生存率。

生存概率与死亡概率

生存概率（probablity of survival）表示某时段开始时存活的个体，到该时段估结束时仍存活的可能性。如年生存要纺表示年初尚存人口存活满一年的可能性。如下所示：

$p=\frac{活过该年人数}{某年年初人数}$

死亡概率（probablity of death）表示某时段开始时存活的个体，在该时段内死亡的可能性。如年死亡概率表示年初尚存人口在今后1年内死亡的可能性，如下所示：

$q=\frac{该年内死亡人数}{某年年初人口数}$

显然，死亡概率$q=1-p$。

注：死亡概率与一般死亡率（death rate或mortality rate）有所不同，如年死亡率：

$m=\frac{该年内死亡人数}{某年平均人口数}$

其分母不同，为”某年平均人口数”，但分子相同。

风险函数

风险率（hazard rate）或风险函数（hazard funciton）记为h(t)，定义如下所示：

$f(t)=\lim_{\Delta \to 0} \frac{P[(t\leq T< t + \Delta t|(T\geq t))}{\Delta t}=\frac{f(t)}{S(t)}$

公式解读：已生存到时间t的观察对象在t时刻的瞬时死亡率。为非负数，且可以$>1$。当Δt=1，$h(t)\approx P(t\leq T < t+1|T \geq t)$，此时的风险率就是时刻t存活的个体在此后一个单位时段内的死亡概率。

累积风险函数（cumulative hazard function），记为H(t)，与和正函数之间的关系为H(t)=-lnS(t)

生存分析的两种方法

Kaplan-Meier法

R中与生存分析有关的包是survival，会涉及到三个函数，如下所示：

Surv：用于创建生存数据对象
survfit：创建KM生存曲线或是Cox调整生存曲线
survdiff：用于不同组的统计检验，其中rho数量参数用以指定检验类型，0为log-rank法或Mantel-Haenszel法，1为Wilcoxon法，默认为0。

案例1：Kaplan-Meier法，生存曲线，

例19-2 为了比较不同手术方法治疗肾上腺肿瘤的疗效，某研究者随机将43例病人分成两组，甲组23例、乙组20例的生存时间（月）如下所示：

原始数据如下所示：

甲组：1，3，5（3），6（3），7，8，10（2），14+，17，19+，20+ ，22+，26+，31+，34， 34+，44，59

乙组：1（2），2，3（2），4（3）， 6（2），8，9（2），10，11，12，13，15，17，18
其中有“+”者是删失数据，表示病人仍生存或失访，括号内为重复死亡数。试计算甲组的生存率与标准误。

计算过程如下所示：

导入原始数据

library(survival)
library(coin)
# 导入例题19-02的数据文件
data1902 <- read.csv("https://raw.githubusercontent.com/20170505a/raw_data/master/data_szq_1902.csv",sep=",")
str(data1902) # 查看数据结构

如下所示：

可以发现，group都是整数型，需要将其转换为因子类型的数据：

1 2	data1902$group <-as.factor(data1902$group) #将分组转化为因子类型的数据 str(data1902)

转换结果如下所示：

创建生存分析的对象：

Surv(data1902$t,data1902$status==1) 
#此语句的意思是创建一个生存分析对象；
#status==1表示完整数据，0为删失数据；
# 其中"+"是删失的数据，表示病人仍生存或失访

生存分析对象如下所示：

生存分析各种参数进行计算：

1
2
3

data.fit <- survfit(Surv(t,status==1)~group,data=data1902) 
summary(data.fit)
# 查看分析结果，group=1表示甲组，2表示乙组

计算结果如下所示：

Surv()结果说明：

time表示的是生存时间，这里与教材中的不太一样，教材中的信息比较全，而R中只列举了有死亡事件出现的生存时间，例如教材中的14个月有数据，而R中没有，因为在14月这个生存时间内，没有死亡，但R中没有列出，没有教材中的详细
n代表处于风险中的人数，这个数字就是上一个生存时间内的人数减去上一个生存时间内的死亡人数
n.event代表的是事件数，默认的status==1表示的是1代表死亡；
survival表示的是生存率
std.err表示的是生存率的标准误；
lower 95% CI：下95%置信区间；
upper 95% CI：上95%置信区间。

Log-Rank检验：

这里使用的是Log-Rank检验法，

1 2	survdiff(Surv(t,status==1)~group,data=data1902) # 比较两条生存曲线的差异，用的是Log-rank检验

如下所示
：

可以发现，${\chi}^2$ 值为8.8，p值为0.00309，可见这两种组别的数据有差异。

绘制生存曲线：

plot(data.fit,col=c("red","blue"),lwd=1) #绘图
legend("topright",legend=c("甲组","乙组"),col=c("red","blue"),lty=1,cex=0.8,adj=-0.7,bty="o") 
#加图例

如下所示：

或者使用survminer包中的ggsurvplot生成图片：

1
2
3

library(survminer)
data.fit <- survfit(Surv(data1$time~data1$status)~data1$group)
ggsurvplot(data.fit,legend.title = "Sex",legend.labs = c("甲组", "乙组"))

如下所示：

中位数生存时间

中位数生存时间（median survival time）又称为生存时间的中位数，表示刚好有50%的个体其存活期大于该时间。

计算中位数生存时间的方法有两种：

①在纵轴生存率为0.5处画一条与横轴平等的线，并与生存率曲线相交，然后自交点画垂线与横轴相交，此效应对应的时间点即为中位生存时间，但结果较粗略，如下所示：

plot(data.fit,col=c("red","blue"),lwd=1) #绘图
legend("topright",legend=c("甲组","乙组"),col=c("red","blue"),lty=1,cex=0.8,adj=-0.7,bty="o") 
#加图例 
abline(h=0.5)

②线性内插法，原理见书，但R中有类似的方法，如下所示：

a <- which(data1902$group==1)
b <- which(data1902$group==2)
data1902[a]
data1902[b]
median_a <- Surv(data1902[a,]$t,data1902[a,]$status==1)
median_b <- Surv(data1902[b,]$t,data1902[b,]$status==1)
summary(median_a)
summary(median_b)

可以看到，甲组的中位数生存时间为10个月；乙组为7个月，跟书中略有出入。

案例分析补充

案例1

某研究者分别用免疫疗法、药物与免疫结合疗法治疗黑色素瘤患者，经随访得到各患的生存时间（月）见下表（孙振球. 医学统计学习题解答-第4版[M]. 人民卫生出版社, 2015.P136）。

试彩乘积极限法计算其生存率及其标准误。
对两组的生存率进行log-rank可检验。
绘制生存曲线。

计算过程如下所示：

library(survival)
library(coin)
data1<- read.csv("https://raw.githubusercontent.com/20170505a/raw_data/master/data_szq_survival_case1.csv",sep=",")
data1$group <-as.factor(data1$group) #将分组转化为因子类型的数据
Surv(data1$time,data1$status==1)
data.fit <- survfit(Surv(data1$t,data1$status==1)~data1$group) 
summary(data.fit)
survdiff(Surv(time,status==1)~group,data=data1) 
plot(data.fit,col=c("red","blue"),lwd=1) #绘图
legend("topright",legend=c("BCG治疗组","药物和BCG治疗组"),col=c("red","blue"),lty=1,cex=0.8,x.intersp=1,bty="o")

标准误与生存率：

两组进行比较：

p值为0.811，还不能认为两种治疗方法有差异。

生存曲线：