【Rcode】生存分析： KM & COX回归& 随机森林& nomogram
本文使用了R——survival包下的lung数据进行方法测试：
1）KM
2) Cox回归：多因素
3）随机森林因子（根据cox回归结果）
4）nomogram
（ 根据cox回归结果，建立了中位生存时间，1年5年生存率的概率计算）
-----🧚🏽♀️本文只有干货，非常干！🤣-----
一、数据加载
#生存分析
library("survival")
library("survminer")
data("lung")
head(lung)
结果：
<head(lung)
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1 3 306 2 74 1 1 90 100 1175 NA
2 3 455 2 68 1 0 90 90 1225 15
3 3 1010 1 56 1 0 90 90 NA 15
4 5 210 2 57 1 1 90 60 1150 11
5 1 883 2 60 1 0 100 90 NA 0
6 12 1022 1 74 1 1 50 80 513 0
二、KM方法和非参数检验
#---------------KM方法-------------------------
fit<-survfit(Surv(time,status)~sex,data=lung) #拟合生存函数
print(fit)
res.sum<-surv_summary(fit)#看整体结果
head(res.sum)
#绘制KM曲线
ggsurvplot(fit,
pval=TRUE,#添加P值
conf.int = TRUE,#增加置信区间
risk.tabel.col="strata", #change risk tabel color by groups
linetype = "strata", #change line type by groups
surv.median.line = "hv", #specify median survival 增加中位生存时间
ggtheme = theme_bw(), #change ggplots theme
palette = c("hue"),
#palette = c("#E7B800","#2E9FDF"),#自定义调色板
risk.table = TRUE,#add risk tabel绘制累计风险曲线
xlab = "Follow up time(d)", # 指定x轴标签
legend = c(0.8,0.75), # 指定图例位置
legend.title = "", # 设置图例标题
legend.labs = c("Male", "Female"), # 指定图例分组标签
break.x.by = 100
#fun="event"
)
结果：
Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
n events median 0.95LCL 0.95UCL
sex=1 138 112 270 212 310
sex=2 90 53 426 348 550
> res.sum<-surv_summary(fit)#看整体结果
三、cox方法和结果
#----------------cox方法-------------------------
data("lung")
head(lung)
res.cox<-coxph(Surv(time,status)~age+sex+ph.ecog,data=lung) #多因素COX回归：age+sex+ph.ecog对生存结局的影响
#这里先做单因素分析纳入P<0.05的，或者临床认为显著相关的
summary(res.cox)
ggsurvplot(survfit(res.cox),
data=lung,
palette = "#2E9FDF",
ggtheme = theme_minimal(),
legend="none"
)
结果：
Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)
n= 227, number of events= 164
(1 observation deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.011067 1.011128 0.009267 1.194 0.232416
sex -0.552612 0.575445 0.167739 -3.294 0.000986 ***
ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0111 0.9890 0.9929 1.0297
sex 0.5754 1.7378 0.4142 0.7994
ph.ecog 1.5900 0.6289 1.2727 1.9864
Concordance= 0.637 (se = 0.025 )
Likelihood ratio test= 30.5 on 3 df, p=1e-06
Wald test = 29.93 on 3 df, p=1e-06
Score (logrank) test = 30.5 on 3 df, p=1e-06
#根据cox回归结果，绘制基础森林图
ggforest(res.cox, data = lung,
main = 'Hazard ratio of LIHC', #标题
cpositions = c(0.05, 0.15, 0.35), #前三列距离
fontsize = 1, #字体大小
refLabel = 'reference', #相对变量的数值标签，也可改为1
noDigits = 3 #保留HR值以及95%CI的小数位数
)
四、诺莫图：根据cox方法
先安装和加载包
#install.packages("rms")
library(rms)
构建COX比例风险模型
dd=datadist(lung)
options(datadist="dd")
f2 <- psm(Surv(time,status) ~ age+sex+ph.ecog,data = lung, dist='lognormal')
med <- Quantile(f2) # 计算中位生存时间
surv <- Survival(f2) # 构建生存概率函数
绘制COX回归中位生存时间的Nomogram图
nom <- nomogram(f2, fun=function(x) med(lp=x),funlabel="Median Survival Time")
plot(nom)
绘制COX回归生存概率的Nomogram图
nom2 <- nomogram(f2, fun=list(function(x) surv(365, x),
function(x) surv(1825, x),
function(x) med(lp=x)),
funlabel=c("1-year Survival Probability", "5-year Survival Probability","Median Survival Time"))
plot(nom2, xfrac=.2)
