R语言|生存分析：Weibull回归

小易学统计

生物统计师|自立，爱家人

各位小伙伴们，大家好！最近很多同学问到关于 weibull回归 在生存分析中如何应用的问题，整理一些资料和R代码，分享一下自己 浅薄的经验 ，其中不乏有个人见解，仅供大家参考。

提起生存分析，很多人在做生存预后因素的时候，几乎想都不想就选择cox回归。为什么呢？个人看法是，主要因为cox回归简单易用，几乎不大用考虑什么前提条件，而weibull分布要检验数据是否符合分布。

本文采用的是survival包中函数来拟合weibull分布，并和KM拟合对比，同时分析它在预后中危险因素的分析应用。

案例实操

本案例数据来R加载数据集lung，因变量为生存资料，自变量是患者的ecog评分，研究不同ecog评分的患者的生存率随时间的变化趋势，以及对预后的影响。

1 加载数据

library(survival)
dt <- read.csv('lung.csv')

2 无变量拟合

2.1 构建模型

data <- dt[dt$status==1,]
fit0 <- survfit(Surv(time, status)~ 1, data=data)
weibull估计
fit1 <- survreg(Surv(time, status) ~ 1, data=data)

2.2 计算拟合值

由于不能直接对fit1作图，因此需要利用其他的函数，提取拟合曲线坐标信息，即拟合的生存率。提取拟合数据的方法有很多种，这里经过前后对比，给大家演示2种简单好操作的方法。

2.2.1 方法一

pct <- 1:98/100  ##定义死亡率
ptime <- predict(fit1,newdata = data.frame(1),type='quantile',p=pct, se=T)
temp0 <- data.frame(time=ptime$fit,
                    ymin=ptime$fit -1.96*ptime$se.fit,
                    ymax=ptime$fit +1.96*ptime$se.fit,
                    p=1-pct)

这个ptime计算的就是pct对应的生存时间，返回的结果有生存时间以及生存时间的标准误，可以计算出生存时间的95%置信区间。提取ptime中的数据，并组装成一个datafram。

绘制曲线

plot(fit0,conf.int=F,col = 4,lwd = 2,
     main="Kaplan-Meier estimate and Weibull",
     xlab = "months", ylab = "survival function",
     cex.main=0.8)
lines(temp0$time, temp0$p, type = "l",lwd = 2,col=9)
lines(temp0$ymin, temp0$p, type = "l",lwd = 2,col=9,lty=3)
lines(temp0$ymax, temp0$p, type = "l",lwd = 2,col=9,lty=3)
legend("topright",pch = c(15,15),legend = c("KM-estimator","Weibull"), 
       col =c(4,9),bty="n")

图形解读：

上图是KM估计和Weibull估计的生存率曲线拟合图，横坐标是生存时间(月)，纵坐标是生存率，中间的蓝色实线是KM拟合曲线，黑色的实线是Weibul拟合生存率的，两边的虚线是95%可信区间。可以看到Weibull拟合的和KM非常的接近。

2.2.2 方法二

利用Weibull分布的参数，计算概率

# survreg's scale = 1/(rweibull shape)
# survreg's intercept = log(rweibull scale)
## reponse rweibull
shape <- 1/fit1$scale   ##weibull分布形状参数
scale <- exp(fit1$coefficients[1]) ##weibull分布尺度参数
p <- 1-pweibull(sort(data$time),shape = shape,scale = scale)