生存资料Cox回归校准曲线绘制（2）

weixin_49320263

已于 2023-08-04 23:31:24 修改

阅读量321

点赞数

分类专栏：模型评价文章标签： r语言数据分析

于 2023-08-04 11:29:43 首次发布

本文链接：https://blog.csdn.net/weixin_49320263/article/details/132099984

版权

模型评价专栏收录该内容

8 篇文章

订阅专栏

原始数据参考文章《生存资料Cox回归校准曲线绘制（1）》

方法3：

canshu<- val.surv(fit= model, # 模型
              newdata = veteran, # 测试集
              u=100, # 时间
              evaluate = 10,#设置点的数目以获取校准曲线
              S=Surv(veteran$time,veteran$status) 
)
plot(canshu,
     xlab="Predicted Probability of Surviving 100 Days",
     ylab="Actual Probability of Surviving 100 Days")

方法4：

#模型预测概率
train_pre <- c((summary(survfit(model, newdata=train), times=100)$surv))
head(train_pre)
#切分四个点
cut <- unique(quantile(c(0, 1, train_pre), seq(0, 1, length = 5), na.rm = TRUE))
cut
#用K-M法估计生存概率
library(rms)
km_surv<- groupkm(train_pre,
                   Srv = Surv(train_df$time,train_df$status),
                   u = 100,
                   cuts = cuts)
km_surv
#简单图
plot(km_surv[,1], km_surv[,4],xlim=c(0,1),ylim=c(0,1))
lines(km_surv[,1], km_surv[,4])
#美化
plot(km_surv[,1], km_surv[,4],
     xlim=c(0,1),ylim=c(0,1),
     xlab = 'Predicted 100-day Survival Probability',
     ylab = 'Observed 100-day Survival Probability',
     main="Calibration curves",
     col="red"
);lines(km_surv[,1], km_surv[,4],col="red")
# 计算误差线范围
errl <- ifelse(km_surv[,"KM"] == 0, 0,  
               km_surv[,"KM"] * exp(1.959964 * (-km_surv[,"std.err"])))
errh <- ifelse(km_surv[,"KM"] == 0, 0, 
               pmin(1, km_surv[,"KM"] * exp(1.959964 * km_surv[,"std.err"])))
# 添加误差线
errbar(x = km_surv[,"x"],
       y = km_surv[,"KM"],
       yminus = errl,yplus = errh,
       add = T,
       pch=16,cex=1,
       asp=1,xaxs='i',yaxs='i',col="red")
# 添加对角线
abline(a = 0,b = 1,col='grey')
#继续添加
#用K-M法估计生存概率
library(rms)
km_surv_1<- groupkm(train_pre,
                  Srv = Surv(train_df$time,train_df$status),
                  u = 200,
                  cuts = cuts)
km_surv_1
lines(km_surv_1[,1], km_surv_1[,4],col="blue")
# 计算误差线范围
errl <- ifelse(km_surv_1[,"KM"] == 0, 0,  
               km_surv_1[,"KM"] * exp(1.959964 * (-km_surv_1[,"std.err"])))
errh <- ifelse(km_surv_1[,"KM"] == 0, 0, 
               pmin(1, km_surv_1[,"KM"] * exp(1.959964 * km_surv_1[,"std.err"])))
# 添加误差线
errbar(x = km_surv_1[,"x"],
       y = km_surv_1[,"KM"],
       yminus = errl,yplus = errh,
       add = T,
       pch=16,cex=1,
       asp=1,xaxs='i',yaxs='i',col="blue")
#添加图例
legend("topleft",legend = c("100 days","200 days"),fill=c("red","blue"))