生存分析
K-M分析
- 定义：KM分析即乘积极限法（product-limit method），是现在生存分析最常用的方法。是由Kaplan和Meier于1958年提出，因此称Kaplan-Meier法，简称KM法。KM法这样估计生存曲线：首先计算出活过一定时期的病人再活过下一时期的概率（即生存概率），然后将逐个生存概率相乘，即为相应时段的生存率。
- R包实现：survival和survminer是其中最常用的两个包，其中survival负责分析，survimner负责可视化。survival提供了Surv对象和survfit()函数用于进行KM生存分析；survminer提供了ggsurvplot()函数**用于KM生存分析结果的可视化。根据性别这个二分类变量，采用KM算法来估计生存曲线
- 示例
library("survival")
library("survminer")
#载入并查看数据集
data("lung")
attach(lung)
fit <- survfit(Surv(time, status) ~ sex, data = lung)
## 作图
ggsurvplot(fit,
pval = TRUE, # 在图上添加log rank检验的p值
conf.int = TRUE,# 添加置信区间
risk.table = TRUE, # 在图下方添加风险表
risk.table.col = "strata", # 根据数据分组为风险表添加颜色
linetype = "strata", # 改变不同组别的生存曲线的线型
surv.median.line = "hv", # 标注出中位生存时间
ggtheme = theme_bw(), # 改变图形风格
palette = c("#E7B800", "#2E9FDF")) # 图形颜色风格
Log-rank检验
- 当两组或多组生存曲线进行比较时,常用的假设检验方法是对数秩检验(log-rank test ) : Log-rank又称时序检验,属于非参数检验,用于比较两组或多组生存曲线或生存时间是否相同，检验统计量为卡方。（注意：这里的 log 表示 count、register 或 record，与对数毫无关系，也称为时序检验/对数秩检验）
- 应用： 检验两条或多条生存曲线的发病率是否相同。与普通卡方检验不同之处是：log-rank 检验能充分利用生存时间（包括删失数据），而且能对各组的生存曲线作整体比较。log-rank 检验是比较两条或多条生存曲线的非参数方法，属于单因素分析方法
- 应用条件：除比较因素外，影响生存率的各协变量组间均衡可比。如果需要调整协变量或同时分析众多因素对生存结局和生存时间的影响，应采用 Cox 比例风险回归模型。
- R代码示例
## log-rank test
survdiff(Surv(time, status) ~ sex, data = lung)
## add p-value in fig
ggsurvplot(fit,pval = TRUE)
# 提取（Log-rank）检验的P值
sd = survdiff(Surv(time, status) ~ sex, data = lung)
1 - pchisq(sd$chisq, length(sd$n) - 1)
# 估测在某一段时间内的生存概率
summary(survfit(Surv(time, status) ~ sex, data = lung), times = 180)
Cox比例风险模型
- Cox 比例风险模型（cox proportional-hazards model）简称Cox 模型，是英国统计学家 D.R.Cox（1972）提出的一种半参数回归模型，可以用来预测一个或多个不同变量在某一时间对死亡率的影响。它同时适用于数值变量和类别变量，可以同时评估几种风险因素对生存时间的影响，检验特定因素如何影响特定时间点特定事件（例如，感染，死亡）的发生率。
- Cox比例风险模型以生存结局和生存时间为应变量，可同时分析众多因素对生存期的影响，能分析带有截尾生存时间（又称删失值或终检值。生存时间观察过程的截止不是由于死亡事件，而是由于其他原因引起的，称为截尾）的资料，且不要求估计资料的生存分布类型。
- Cox模型是根据多种特征（即协变量xi），在一个随时间变化的基准风险函数基础上建立的多元线性回归模型，可以用来评估具有特定特征值的患者在某一时刻的瞬时风险概率。
- COX比例风险回归模型(Proportional hazards model)常用于多因素生存分析，探索影响多因素对生存期的影响。该模型是一种半参数回归模型，对数据分布要求较低。
- R示例
library("survival")
library("survminer")
data("lung")
head(lung)
# 单变量Cox分析可以如下计算：
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox
summary(res.cox)
# 要将单变量coxph函数一次应用于多个协变量，请输入以下命令：
covariates <- c("age", "sex", "ph.karno", "ph.ecog", "wt.loss")
univ_formulas <- sapply(covariates,
function(x) as.formula(paste('Surv(time, status)~', x)))
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
## 提取数据，并制作数据表格
univ_results <- lapply(univ_models,
function(x){
x <- summary(x)
p.value<-signif(x$wald["pvalue"], digits=2)
wald.test<-signif(x$wald["test"], digits=2)
beta<-signif(x$coef[1], digits=2);#coeficient beta
HR <-signif(x$coef[2], digits=2);#exp(beta)
HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
HR <- paste0(HR, " (",
HR.confint.lower, "-", HR.confint.upper, ")")
res<-c(beta, HR, wald.test, p.value)
names(res)<-c("beta", "HR (95% CI for HR)", "wald.test",
"p.value")
return(res)
#return(exp(cbind(coef(x),confint(x))))
})
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)
- COX 模型结果整理
## cox regression
coxph(Surv(time, status) ~ sex, data = lung)
# 变成表格形式
coxph(Surv(time, status) ~ sex, data = lung) %>%
gtsummary::tbl_regression(exp = TRUE)
# 风险比：HR = exp(β)
## cox regression
f4<-coxph(Surv(days,status==1)~sex)
summary(f4)
#summary(coxph(Surv(days,status==1)~sex))
f5<-coxph(Surv(days,status==1)~sex+log(thick)+ulc)
summary(f5)
#summary(coxph(Surv(days,status==1)~sex+log(thick)+strata(ulc)))
plot(survfit(coxph(Surv(days,status==1)~
log(thick)+sex+ulc)),col=c("red","blue"))
#plot(survfit(coxph(Surv(days,status==1)~
#log(thick)+sex+strata(ulc))),col=c("red","blue"))
legend(locator(n=1),legend=c("ulceration present","ulceration absent"),lty=1,col=c("red","blue"))
竞争风险
定义：竞争风险是指在观察队列中，存在某种已知事件可能会影响另一事件发生的概率或者完全阻碍其发生，则可认为前者与后者存在竞争风险。竞争风险属于临床研究关注的临床结局（称之为兴趣终点）之外的另外一个结局，竞争风险的发生率与兴趣终点相等，甚至高于兴趣终点，同时它可以影响兴趣终点的发生的可能性。因此，竞争风险有三个属性，A，不是兴趣终点；B，临床常见；C，与兴趣终点的发生不相互独立。三者同时具备才能称之为竞争风险。
当存在竞争事件时，传统的做法，通常是将竞争事件的发生视作刪失，来计算结局事件的生存概率。如果采用cox模型分析，会导致结局事件的生存概率估计有偏差，一般会高估累计风险发生率，因此在存在竞争风险时，我们优先使用竞争风险模型。
竞争风险的单因素分析常用来估计关心终点事件的发生率，多因素分析常用来探索预后影响因素及效应值。
竞争风险模型包括单因素和多因素的分析：
- 单因素竞争风险模型 (Single-factor competition risk model)：估计目标事件与竞争事件的累积发生率；组间比较时采用 Gray检验。仅考虑单个事件时，类似经典生存分析的Kaplan-Meier 法估计事件的累积发生率。
- 多因素竞争风险模型 (Multi-factor competitive risk model): 多因素分析时,运用Fine-Gray模型探讨影响目标事件累积发生率的因素。类似经典COX模型探讨影响目标事件累积发生率的因素。
- 竞争风险模型适用于多个终点的生存数据。其单因素分析常为估计关心终点发生率（对应生存分析中常用Kaplan-Meier边际回归法（KM法）估计死亡率），多因素分析常为估计预后影响因素及效应值（对应生存分析中常用cox比例风险模型中估计死亡影响因素及效应）。
- R中竞争风险
##fine-gray test
#library(foreign)
bmt <-read.csv('./ch07/bmtcrr.csv')
head(bmt)
bmt$D <- as.factor(bmt$D)
# install.packages("cmprsk")
library(survival)
library(cmprsk)
library(splines)
attach(bmt)
crmod <- cuminc(ftime,Status,D)
crmod
plot(crmod,xlab = 'Month', ylab = 'CIF',lwd=2,lty=1,
col = c('red','blue','black','forestgreen'))
##cmprsk model
cov1 <- data.frame(age = bmt$Age,
sex_F = ifelse(bmt$Sex=='F',1,0),
dis_AML = ifelse(bmt$D=='AML',1,0),
phase_cr1 = ifelse(bmt$Phase=='CR1',1,0),
phase_cr2 = ifelse(bmt$Phase=='CR2',1,0),
phase_cr3 = ifelse(bmt$Phase=='CR3',1,0),
source_PB = ifelse(bmt$Source=='PB',1,0)) ## 手动设置哑变量
cov1
mod1 <- crr(bmt$ftime, bmt$Status, cov1, failcode=1, cencode=0)
summary(mod1)
# install.packages("aod")
library(aod)
wald.test(mod1$var,mod1$coef,Terms = 4:6)
