生存分析

http://www.bio-info-trainee.com/1313.html

http://www.bio-info-trainee.com/1371.html

http://www.dxy.cn/bbs/thread/3327209#3327209

http://blog.csdn.net/shmilyringpull/article/details/17529637

http://www.biomart.cn/experiment/430/586/588/240451.htm?trace=0420labstp

生存分析(survivalanalysis)适合于处理时间－事件数据。例如中风病人从首次发病到两次复发，其中就涉及到时间和事件。此例中时间就是复发的时间间隔，事件就是是否复发。

一般我们谈生存分析，就是说的KM（Kaplan Meier）方法估计生存函数，并且画出生存曲线，然后还可以根据分组检验一下它们的生存曲线是否有显著的差异，甚至，我们还可以用cox模型来分析这个因子是如何影响生存函数的。

Surv：用于创建生存数据对象
survfit：创建KM生存曲线或是Cox调整生存曲线
survdiff：用于不同组的统计检验
coxph：构建COX回归模型
cox.zph：检验PH假设是否成立
survreg：构建参数模型

一、生存分析基本概念

1、事件（Event）

指研究中规定的生存研究的终点，在研究开始之前就已经制定好。根据研究性质的不同，事件可以是患者的死亡、疾病的复发、仪器的故障，也可以是下岗工人的再就业等等。
2、生存时间(Survival time)
指从某一起点到事件发生所经过的时间。生存是一个广义的概念，不仅仅指医学中的存活，也可以是机器出故障前的正常运行时间，或者下岗工人再就业前的待业时间等等。有的时候甚至不是通用意义上的时间，比如汽车在出故障前的行驶里程，也可以作为生存时间来考虑。
3、删失（Sensoring）
指由于所关心的事件没有被观测到或者无法观测到，以至于生存时间无法记录的情况。常由两种情况导致：（1）失访；（2）在研究终止时，所关心的事件还未发生。
4、生存函数（Survivaldistribution function）又叫累积生存率，表达式为S（t）=P(T>t),其中T为生存时间，该函数的意义是生存时间大于时间点t的概率。t=0时 S(t)=1，随着t的增加S(t)递减（严格的说是不增），1-S(t)为累积分布函数，表示生存时间T不超过t的概率。
二、生存分析的方法
1、生存分析的主要目的是估计生存函数，常用的方法有Kaplan-Meier法和寿命表法。对于分组数据，在不考虑其他混杂因素的情况下，可以用这两种方法对生存函数进行组间比较。
2、如果考虑其他影响生存时间分布的因素，可以使用Cox回归模型（也叫比例风险模型），利用数学模型拟合生存分布与影响因子之间的关系，评价影响因子对生存函数分布的影响程度。这里的前提是影响因素的作用不随时间改变，如果不满足这个条件，则应使用含有时间依存协变量的Cox回归模型。
下面用一个例子来说明SPSS中Cox回归模型的操作方法。

1.      KM生存分析

只需要记住三个函数即可：

Surv：用于创建生存数据对象。

survfit：创建KM生存曲线或是Cox调整生存曲线。

survdiff：用于不同组的统计检验。

我们用下面的数据做例子:

无病生存期（Disease-free  survival,DFS）的定义是指从随机化开始至疾病复发或由于疾病进展导致患者死亡的时间。该指标也常作为抗肿瘤药物III期临床试验的主要终点。某些情况下，DFS与OS相比，作为终点比较难以记录，因为它要求认真随访，及时发现疾病复发，而且肿瘤患者的死亡原因也很难确定（16）。肿瘤患者常有合并症（如，心血管病），这些合并症可能会干扰对DFS的判断。并且，肿瘤患者常死于医院外，不能常规进行尸检。

总生存期（Overall survival,OS）的定义是指从随机化开始至因任何原因引起死亡的时间。该指标常常被认为是肿瘤临床试验中最佳的疗效终点。如果在生存期上有小幅度的提高，可以认为是有意义的临床受益证据。作为一个终点，生存期应每天进行评价，可通过在住院就诊时，通过与患者直接接触或者通过电话与患者交谈，这些相对比较容易记录。确认死亡的日期通常几乎没有困难，并且死亡的时间有其独立的因果关系。当记录至死亡之前的失访患者，通常截止到最后一次有记录的、与患者接触的时间。

library(survival)

my.surv <-Surv(OS_MONTHS,OS_STATUS==’LIVING’)#这个生存对象是看看病人的总生存期与死亡状态的关系, 这个Surv函数第一个参数必须是数值型的时间，第二个参数是逻辑向量，1,0表示死亡与否

kmfit1<- survfit(my.surv~1) #直接对生存对象拟合生存函数

summary(kmfit1)

plot(kmfit1) #画出生存曲线

##根据分组type估计KM生存曲线

my.surv<- Surv(OS_MONTHS,OS_STATUS==’LIVING’)

kmfit3<- survfit(my.surv~type)

summary(kmfit3)

plot(kmfit3,col= rainbow(4))

##显著性分析

survdiff(my.surv~type,data=dat)

2.     cox生存分析的风险因子（比例风险模型）:用cox方法看因子的权重！

数据：表达分组和甲基化分组

把这两种分组方式当做因子来探究它们对生存率的影响！

my.surv=Surv(sur_dat$OS_MONTHS.y,sur_dat$OS_STATUS.y==’LIVING’)

#plot(survfit(my.surv~1))

survfit(my.surv~1)

kmfit=survfit(my.surv~1)

plot(kmfit)

summary(kmfit)

用survdiff检验分组的显著性，结果如下：

如果用cox模型回归分析如下：

用Sorafenib比不用将使患者的死亡概率减少一半。可是要是结合生存曲线，我们可以说，用Sorafenib能延长一倍的生存时间。

关于HR进一步探讨

http://www.dxy.cn/bbs/topic/27383879

看见前面很多前辈们发过有关HR的帖子，深受启发！现将本人的学习经验心得上传，还望多多指教！
一、为何要选择风险比
生存资料的Meta分析，顾名思义，其基于的数据类型为生存资料。作为时间相关事件的数据(time-to-event data)主要类型之一, 生存资料主要是一组既能记录某一事件(如肿瘤死亡)的发生，同时也能反映出现这一结果所经历的时间的数据。该类数据的分析不仅考虑事件是否会出现，而且也考虑事件出现的时间长短。其中，该类数据最具代表性的是生存时间分析，如肿瘤患者化疗后的生存率，冠心病患者两次发作之间的间隔，艾滋病患者从确诊阳性到最后死于AIDS及其并发症等。通常，该类数据的来源是采用纵向观察随访获取，且其分布不服从正态分布，同时还存在着数据的删失(即随访时间内未观察到事件的发生)。
如有关肿瘤病人预后的Ⅲ期临床试验，主要的研究指标大多为总生存时间(OverallSurvival, OS)或疾病无进展时间(Progress Free Survival, PFS)等生存相关的数据。这就要求在进行数据合并时不仅仅是要考虑事件的发生，且还要考虑发生该事件所需要的时间。然而，当前有关生存资料Meta分析在多数情况下仅通过对纳入研究中的事件发生数进行某一个时间点上的点估计，通过计算不同组间事件发生率的比值进行描述，即运用相对危险度 (Relative Risk, RR)或者比值比(OddsRatio, OR)进行分析。这样单纯地运用某一个固定时间点两组的生存人数的事件发生率的比值来测量，并没有完全考虑到所有的因素，根本无法描述该类数据的全貌，且这种不全面的数据有可能会得出不恰当的结论。这就要求我们需要一种既要能记录发生事件的结果，同时也能体现该事件发生所经历的时间的测量指标---风险比 (HazardRatio, HR)。
在Meta分析中，HR的意义是指某一种干预措施(试验组)的应用所产生的风险率与不用该干预措施或者空白对照以及安慰剂等对照时所产生的风险率的比值。它本身是一个和时间无关的量。但是，相比于对中位生存时间进行T-检验或者线性回归，运用HR生存分析考虑了删失数据，且相对于单纯使用RR、OR或者是 Logistic回归而言，运用HR生存分析考虑了时间因素。进而能尽可能真实地反映该类数据所代表的结果。

二、风险比的原理及意义
对于原始研究的随访过程中，难免会出现数据的删失。当原始数据中不存在删失数据，其生存率(Survival Probability)的计算公式为 t时刻仍存活的病人数与观察总例数的比值，此类数据虽表示生存的概率，在大小上等同于该事件的发生率。当一个研究存在删失数据，则需分段计算各时间段的生存概率，运用概率乘积原理将各个时段的生存概率进行累积相乘而得。因为对于各个时段的事件发生是独立的，且每一个时段的事件发生的概率不全相等。对不同时段的不尽相同的生存概率进行累加，一方面考虑了事件的结局，同时也反映了发生该类事件的时间。

下面将用一个简单的例子进行说明：表1表示某位神经外科大夫收集的10例恶性肿瘤患者术后接受某干预措施的随访生存时间的结果，以此来估计该组患者的生存率。

由表2可以知道，每一个时间段的生存概率不尽相同，且总的生存率是各个时间段的生存概率的乘积。该种不同时间段生存概率的累积的方法，既考虑了事件的发生结果，同时也考虑了发生该事件所需的时间。故而，HR优于RR/OR。
HR的流行病学意义是指变量X_j暴 露水平时的风险率(HazardRate)与非暴露水平时的风险率之比，即分别具有协变量的两个体，其风险函数(HazardFunction)的比值。且该风险函数与相对的生存函数(SurvivalFunction)的对数之间存在负线性相关。

How to calculate HR?

摘抄自http://stats.stackexchange.com/questions/124489/how-to-calculate-the-hr-and-95ci-using-the-log-rank-test-in-r

library(survival)data.survdiff <- survdiff(Surv(time, status) ~ group)p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))

Output results:

> data.survdiffCall:survdiff(formula = Surv(data[, "os_whw"], data[, "status_whw"] ==     1) ~ data[, "pcascore"] >= median(data[, "pcascore"]))                                                       N Observed Expected (O-E)^2/E (O-E)^2/Vdata[, "pcascore"] >= median(data[, "pcascore"])=FALSE 4        3     4.33     0.411     0.974data[, "pcascore"] >= median(data[, "pcascore"])=TRUE  5        5     3.67     0.486     0.974 Chisq= 1  on 1 degrees of freedom, p= 0.324 > p.val[1] 0.3235935> HR[1] 1.970484> up95[1] 7.917248> low95[1] 0.4904239