8.常用统计分析方法——生存分析

生存分析基本概念

logistic回归中因变量是终点事件发生与否，而生存分析则关注的是终点事件所经历的时间。

生存资料的特点：1.随访资料，包括两个方面时间和结局；2.时间—事件变量；3.有不完全数据；
生存分析：就是用来研究“生存”过程的统计方法；研究生存时间；时间与事件的关系

生存资料的三要素 ：

研究因素
生存时间
生存结局

生存资料的分析 ：

（1）估计：根据样本生存资料估计总体生存率和其他有关指标；

（2）比较：对不同组生存曲线进行比较；

（3）影响因素分析：目的是探索和了解影响预后的因素，或平衡某些因素的影响后，研究某个或者某些因素对生存的影响。

生存时间（survival time） ：根据研究目的确定的观察起点到某一给定终点事件出现的时间。类型：

（1）完全数据：随访研究中，在规定的时间内，若观察到某些对象的死亡结局所经历的时间，称为生存时间的完全数据。

（2）删失数据：随访研究中，在规定的时间内，由于某些原因未能观察到某些对象的死亡结局，提供的是从手术到停止观察的时间长度，并非确切的生存时间，称为生存时间的删失数据（censored data）

删失常见的原因：

研究结束时终点事件未发生
由于未继续就诊、拒绝随访或失去联系等失访，未能观察到其死亡结局
因死于其他原因终止观察

起始事件 ：标志研究对象生存过程开始的特征事件，又称为起点事件，与终点事件相对应。
终点事件 ：指研究者所关心的特定事件，又称为失效事件，标志着“生存”的终止。失效事件不一定是不好额事件，只表示特点事件的发生。
截尾：又称为终检或删失，指在发生终点事件之前，被观察对象的观测过程终止了。没有观察到所关心终点事件的发生。
特点：不完全信息（不知道确切的生存时间），但可知真实的生存时间不会短于现在观察的时间。

生存概率和死亡概率：

生存概率：某单位时段开始存活的个体，到该时段结束时仍存活的可能性。

死亡概率：表示某单位时段开始存活的个体，在该时段内死亡的可能性。

若观察时段相同，p+q=1

生存率与风险率：

生存率（survival rate）也即是生存函数指观察对象活过t个单位时段的可能性，记为S(t)。如果没有删失数据，直接法计算生存率： $S(t)=P(T> t)\approx \frac{t时刻仍存活的例数}{观察总例数}$ 若有删失数据，须分段计算生存概率： $S(t)=P(T> t)=P_{1}\cdot P_{2}\cdot\cdot\cdot\ P_{t}$

风险率（hazard rate）或风险函数表示t时刻存活的个体在t时刻的瞬时死亡风险。

生存率估计

非参数法估计的生存率主要有寿命表法（life table method）和乘积极限法（product limit method），前者适用于 生存时间按区间分组 的大样本资料，后者适用于 仅含个体生存时间 的小样本或者大样本资料。

分析思想：先求出各个时段的生存概率，然后根据概率乘法定理估计生存率。

1. 乘积极限法

又称为K-M法，乘积的含义：生存率等于生存概率的乘积；极限的含义：寿命表法中的时间区间长度趋近于0，也即是有明确的生存时间或截尾时间。

条件：每个对象均有确切的生存时间/截尾时间；tk时刻的生存率等于tk时刻之前各点生存概率的乘积；不需要对资料分布类型进行假设。

基本计算步骤：

将生存时间（ $t_{i}$ ）由小到大排序，完全数据和删失数据相同时，删失数据在完全数据的后面；
列出时间区间[ $t_{i}, t_{i+1}$ ）上的死亡数 $d_{i}$ 和删失数 $c_{i}$ ;
计算期初例数，恰在某一时刻 $t_{i}$ 之前的生存人数 $n_{i}$ ;
计算死亡概率 $q_{i}$ 和生存概率 $p_{i}$ ;
计算生存率 $\hat{S}(t_{i})$ ，注意删失时间处不做变换，删失时间的影响只体现在期初例数的计算上。

（1）创建生存对象 ：将研究对象的结局变量由“time”（时间）和“event”（事件）组成。

生存对象可以由survival包中的 函数Surv（） 生成，是一种将事件和时间合并在一起的数据结局。

案例：比较卵巢癌患者在两种治疗方法下的生存率

library(survival)
library(survminer)
data(ovarian)
str(ovarian)
attr(ovarian,"var.labels")[1:6]<- c("随访时间",
                                   "结局状态：0存活，1死亡",
                                   "患者年龄",
                                   "疾病残留情况：1没有残留，2有残留",
                                   "治疗方法：1环磷酰胺，2环磷酰胺加阿霉素",
                                   "患者的ECOG评分，1较好，2较差")
library(epiDisplay)
des(ovarian)
attributes(ovarian)
summary(ovarian)

对分类变量进行因子化

ovarian$resid.ds <- factor(ovarian$resid.ds,
                           levels = c(1,2),
                           labels = c("no","yes"))
ovarian$rx <- factor(ovarian$rx,
                     levels = c(1,2),
                     labels = c("A","B"))
ovarian$ecog.ps <- factor(ovarian$ecog.ps,
                          levels = c(1,2),
                          labels = c("good","bad"))
hist(ovarian$age)#年龄不对称将，年龄进行因子化
ovarian$agegr <- cut(ovarian$age,
                     breaks = c(0,50,75),
                     labels = c("<=50",">50"))#按50为分界值将，年龄分成两组
table(ovarian$agegr)
创建生存对象： 
surv.obj <- Surv(time =ovarian$futime,event = ovarian$fustat);surv.obj 
在生存对象中，“+”号表示删失数据，例如，第四个值“421+”，表示这个患者并未在421天死于卵巢癌，只是没有继续随访 
(2) 生存率估计 
 生存率估计的基本格式： 
 survfit(Surv(生存时间，结局)~分组变量，data)，survfit根据KM法创建生存曲线，如果没有分组变量，就用"1"代替。 
 其中Surv（生存时间，结局）就是上面提到的生存对象。 
survfit(surv.obj~1) 
依次展示的是样本数，事件发生数，中位生存时间，中位生存时间的置信区间。因为数据样本较小，生存曲线没有穿过生存率在50%对应的水平先，无法计算生存时间置信区间的上限。 
通过summary（）函数对生存分析的数据进行展示，显示具体的各时间段的生存概率和生存率。 
fit <- survfit(surv.obj~1)
summary(fit)#不显示删失数据
summary(fit,censored = T)#包含删失数据 
 参数censored默认未FALSE，通过TRUE可以展示删失数据。包含各时间点的存活对象个数，生存率及置信区间，censored=T包括删失数据 
2. 寿命表法 
资料是分时间段收集，数据按时间区组分组
 对于每个时间段，假定死亡情况在每个时间段均匀出现，即生存率匀速下降
 假定删失的发生在时间段内也是均匀分布，区间中删失个体算半个人时，校正观察例数 
两种方法的选择： 
#寿命表法：适用于时间分段或无法得知准确生存时间的资料，区间内生存率呈线性的变化
#KM法：要求每个观测对象有明确的死亡时间和截尾时间 
3. 生存曲线 
生存曲线是一条横坐标以时间为参数，纵坐标以生存率为参数的统计图。 
（1）基础包plot绘图 
plot(fit,mark.time = T) 
 mark.time=T为删失数据添加标记"+"，不想按展示置信区间conf.int=F 
（2）用survival包函数ggsurvplot（）绘制生存曲线 
ggsurvplot(fit,pval=T,conf.int = T) #生存曲线 
（3）死亡曲线 
ggsurvplot(fit,pavl=T,conf.int = F,fun = "event",ylim=c(0,1)) 
 fun表示绘制不同类型的曲线，event绘制累积事件，cumhaz绘制累积危险函数，pct表示百分比表示生存率 
 Survfit(Surv(生存时间，结局)~分组变量，data)，survfit根据KM法创建生存曲线，如果没有分组变量，就用"1"代替
 ggsurvplot(fit,pval=TRUE,conf.int=FALSE),fit表示生存曲线函数，pval是否将p值显示到曲线上，conf.int是否显示置信区间条带
 surv_diff <- survdiff(Surv(生存时间，结局)~分组变量，data)，按照分组变量进行不同组别生存曲线的比较 
生存曲线比较 
案例：想要比较不同治疗方式下的生存率差异？ 
分析：由治疗方式作为分组，进行上述操作。 
surv.treat <- survfit(surv.obj~rx,data=ovarian)#创建生存对象
summary(surv.treat,censored = T)#计算生存率
plot(surv.treat,mark.time = T, #绘制生存曲线
     conf.int = F,
     lty=c(1,2),
     lwd=2,
     col = c("red","blue"))
legend(60,0.3,title = "KM cur",
       legend=c("环磷酰胺","环磷酰胺+阿霉素"),
       lty=c(1,2),
       lwd=2,
       col = c("red","blue")) 
#survmuner包中的ggsurplot以ggplot2的风格提供更丰富的图形输出
ggsurvplot(surv.treat,data=ovarian,pval=T)
survdiff(surv.obj~agegr,data = ovarian) 
 生存曲线的比较： 
Log-rank检验：时序检验，属于非参数的方法，比较整个生存时间的分布。可以通过函数survdiff（生存对象~分组，data=《MYDATA》）格式计算。 
surv_diff <- survdiff(Surv(futime,fustat)~rx,data=ovarian);surv_diff 
根据上述结果，两种治疗方法在生存率的差异没有统计学差异。但是两组中研究对象可能被其他影响生存的因素如年龄，疾病残留情况，患者的ECOG评分影响，需要使用回归分析，把这些因素作为协变量解释两组之间的差异。 
注意：当不同时间段内生存情况的优劣存在交叉时，log-rank检验不适用，，表现为生存曲线存在明显交叉 
COX比例风险回归模型 
风险率（hazard rate）或风险函数表示t时刻存活的个体在t时刻的瞬时死亡风险。风险函数h(t)用于度量在某个时刻t还存活的个体在极短时间内死亡的风险。h(t)=寿命分布函数f(t)/S(t) 
 COX回归模型的两个前提假设： 
各个危险因素的作用不随时间的变化而变化，即不会随着时间的变化而变化；
对数线性假定：模型中协变量应该与对数风险比成线性关系。 
COX回归分析的一般步骤： 
1. 分析前的准备——数据整理 
2. 参数估计，建立最佳模型 
3. 假设检验 
4. Cox模型的解释及应用 
5. Cox模型拟合优度的考察 
1. 建立COX回归模型 
案例分析： 
将所有协变量都包含进来建立cox回归模型： 
cox1 <- coxph(surv.obj~rx+resid.ds+agegr+ecog.ps,data=ovarian)
summary(cox1) 
 调整协变量后，rx组p值0.0322有统计学差异，因此两种治疗方法下死亡风险有统计学差异。 
与多重线性回归类似，多重cox回归也存在变量选择问题，对包含所有自变量的模型cox1使用函数drop1（）可以得到各个自变量的似然比检验结果。 
drop1(cox1,test = "Chisq") 
计算各个自变量似然比检验，ecog.ps对于的p值最大，从模型中剔除。也可以使用step（）函数基于AIC进行变量选择。 
cox2 <- coxph(surv.obj~rx+resid.ds+agegr,data=ovarian)
step.cox <- step(cox1) 
step的结果也表明cox2是最优模型 
2. 比例风险假定的检验 
验证cox回归的比例风险假定，使用函数cox.zph() 
cox.zph(cox2) 
结果表明各个自变量检验和全局检验p值都大于0.05，没有违背比例风险的假定 
summary(cox2) 
接受方法B治疗的死亡风险是方法A的27.78%，反过来，治疗A方法的死亡风险是B方法的3.599倍，同时不同残留情况以及不同年龄组之间的死亡风险差异均无统计学差异 
3. 生存预测 
newdata <- data.frame(rx = c("A" , "B"), 
                      resid.ds = c("no", "no"),
                      agegr = c(">50", ">50"))
newdata
hr <- predict(cox2, newdata = newdata, type = "risk")
hr[1]/hr[2] 
A治疗比B治疗的风险
 #大于50岁，都没有风险残留的病A治疗的死亡风险是B治疗的3.599倍
 #使用KM拟合生存曲线更直观展示两种治疗方式下的生存率变化 
绘制生存曲线： 
cox.fit <- survfit(cox2, newdata = newdata, type = "kaplan-meier")
plot(cox.fit, lty = c(1, 2), col = c(2, 4))
title(main="Cox survival curves by treatment
      for age > 50, no residual disease patients",
      xlab="Duration in days", 
      ylab="Survival probability",
      las = 1)
legend(5, 0.3, c("Treatment A", "Treatment B"), 
       lty = c(1, 2), col = c(2, 4)) 
生存描述-描述不同时间的总体生存率，计算中位生存时间，绘制生存函数曲线，一般用Kaplan-Meier方法和寿命表法；
生存曲线比较-比较不同处理组的生存率，一般用logrank检验；
生存相关因素的分析：回归模型；由于logrank检验仅能分析一个因素，因此两个或者两个以上因素的分析需要使用Cox比例风险模型；
2、生存分析使用的方法：
Kaplan-Meier plots
				生存情况信息统计分析：
案例题目：查看数据信息，求平均年龄。求男女生存率，未成年人获救的概率。
数据来源：https://www.kaggle.com/c/titanic/data
# coding=utf-8
import numpy as np
import pandas as pd
#中文的话这样打开，不会出现Initializing from file failed这种错误
f=open("./data/train.csv")
titanic=pd.read_csv(f)
titanic.head(
生存分析(Survival analysis)是指根据试验或调查得到的数据对生物或人的生存时间进行分析和推断，研究生存时间和结局与众多影响因素间关系及其程度大小的方法，也称生存率分析或存活率分析。临床医学数据生存分析是将临床医学中事件的现象与结果（终点事件）和出现临床现象与结果经历的时间结合起来分析的一种统计分析方法。
生存分析的目的
1.描述生存过程：预测不同时间的总体生存率、计算中位生存期、绘制生存函数曲线。统计方法包括Kaplan-Meier（K-M）法、寿命表法。
2.比较研究结果：比较不同处
				最近在读一篇医学信息处理相关的论文，由于在医学图像和医学信息相关的计算机处理方面还是个新手，遇到了一个生存分析预后的概念，提到KM生存曲线、COX回归模型这样的部分时读的一知半解，索性就上网学习了一下，留个记录吧。
什么是生存分析？
生存分析是对生存资料的分析。所谓生存资料是指描述寿命或者一个发生时间的数据。通俗来讲，“一个个体能活多久”这个问题扁鹊再世也无法回答，但通过对某一具有相同特质的群体进...
什么是生存分析
生存分析是对生存资料的分析。所谓生存资料是指描述寿命或者一个发生时间的数据。通过对某一具有相同特质的群体进行数据分析，我们可以得到这类人活过一定时间的概率。这就叫做生存分析。一个人的生存时间长短与许多因素有关，研究各个因素与生存时间有无关系以及关联程度大小，这也是生存分析。推广来说，疾病的复发可以看做“健康期”的生存，机器的故障可以看做“机器”的生存，甚至下岗职工再就业可以看做“失业期”的生存，生存分析可以被推广至很多问题的分析中去。
生存分析最早在生物医学中使用的最多，用来
MedCalc是一种常用的统计软件，用于进行医学和生物统计学分析。它提供了多种统计方法和工具，帮助研究人员和医生更容易地分析数据和做出决策。
MedCalc具有以下功能和应用：
1. 描述性统计分析：可以计算和展示数据的中心趋势（如均值、中位数）和离散程度（如标准差、范围），提供对所研究数据的直观理解。
2. 参数估计：通过应用不同的参数估计方法（如置信区间、假设检验）来进行数据分析和解释。
3. 数据可视化：可以生成各种图表（如直方图、散点图、箱线图）以及绘制ROC曲线和生存分析曲线，有助于展示数据的分布和趋势。
4. 假设检验：提供了多种统计检验方法（如t检验、方差分析、Chi-square检验），用于比较不同组别或类别之间的差异和关联。
5. 相关分析：可以计算和展示变量之间的相关性，包括Pearson相关系数和Spearman相关系数。
6. 回归分析：可以进行线性回归、多元回归和逻辑回归分析，用于预测和解释因变量与自变量之间的关系。
7. 生存分析：可用于生成Kaplan-Meier生存曲线，并计算生存率、风险比和生存时间的统计学参数。
总之，MedCalc是一个强大且易于使用的统计分析工具，广泛应用于医学研究、临床实践和生物统计学领域。它帮助用户进行数据分析、结果解释和相关决策，从而提高研究成果的质量和可靠性。  
### 回答2：
MedCalc是一种用于统计分析的软件工具，可用于支持医学和生物医学研究。它提供了各种统计方法和计算工具，帮助研究人员进行数据分析和解读结果。
MedCalc提供了常见的统计方法，如描述性统计、参数和非参数假设检验、回归分析、生存分析等。研究人员可以使用描述性统计方法来计算数据集的平均值、标准差、中位数等统计量，以了解数据的中心趋势和离散程度。参数和非参数假设检验方法可用于比较两组或多组数据是否存在显著差异，其中参数方法包括t检验、方差分析等，非参数方法包括Wilcoxon秩和检验、Kruskal-Wallis检验等。回归分析可以用于研究因变量与自变量之间的关系，生存分析则可以评估生存数据以及影响生存的因素。
MedCalc软件的应用广泛，适用于各种医学和生物医学研究领域。例如，临床试验研究可以使用MedCalc进行样本大小计算、随机化的生成和执行统计分析。此外，MedCalc还支持meta分析，可用于综合多个研究结果，以获得更精确的估计和结果。生物医学研究可以利用MedCalc来分析基因表达数据、多组比较、生物标志物评估等。
总之，MedCalc是一种功能强大的统计分析软件，可满足医学和生物医学领域的数据分析需求。它提供了丰富的统计方法和计算工具，帮助研究人员进行数据处理、假设检验、回归分析等，并支持多种应用领域，包括临床试验和生物医学研究。
                    Jing 我愿。: 
                    这个se计算有误，多了负号而且忘记除2了？ 赞同，而且这段代码过渡的时候，似乎是修改了列名result2$characteristic <- ifelse(is.na(result2$number),
                                 result2$characteristic,
                                 paste0("    ",result2$characteristic))
                实用技巧——绘制森林图
                    怎么控制forest生成图像的像素呢？
                6. 临床预测模型——可视化列线图
                    来点蛋挞778: 
                    请问二分类列线图中肺动脉栓塞的数据在哪可以获得，谢谢表情包
                7. 预测模型评价——DCA曲线
                    m0_67475320: 
                    family设置为二分类资料binomial，请问在一开始是否需要对结局进行因子化，在一开始就对结局进行因子化会发现报错