Cox回归校准曲线(测试集)的实现方法（下）

临床预测模型基础入门必看合集：R语言临床预测模型合集

之前的推文给大家介绍了Cox回归各种校准曲线的实现方法，包括训练集和测试集：

• Cox回归校准曲线(测试集)的实现方法（上）

在最后个大家留了一个疑问，今天继续：

大家经常读文献就会发现这种COX回归测试集的校准曲线↓：

目前好像并没有包可以直接实现，不过也不是非常困难，下面给大家介绍实现方法。

本文目录：

准备数据

首先把数据集划分为训练集、测试集，然后在训练集建立cox模型

rm(list = ls())
library(survival)

lung$status <- ifelse(lung$status == 2,1,0)
lung <- na.omit(lung)

set.seed(123)
ind <- sample(1:nrow(lung),nrow(lung)*0.7)

train_df <- lung[ind,]
test_df <- lung[- ind, ]

# 在训练集建立cox模型
coxph_fit <- coxph(Surv(time, status) ~ age + sex + ph.ecog,
                   data = train_df,
                   x = T, y = T)

训练集的校准曲线

对于生存分析来说，校准曲线的横纵坐标分别是模型预测的生存概率(predicted survival probability)和K-M法估计的生存概率(也可以称为Observed survival probability)，其中模型预测的生存概率是模型算出来的(生存分析一般是根据COX算出来的)，K-M法估计的生存概率则是根据乘积极限法(K-M法)计算的。

模型计算的生存概率很简单可以得到，比如训练集100天的校准曲线，首先计算cox模型预测的100天的生存概率：

# 模型预测的概率
set.seed(123)
train_p <- c((summary(survfit(coxph_fit, newdata=train_df), times=100)$surv))
head(train_p)

0.6958190 0.9347173 0.8681793 0.7350013 0.8694693 0.8213824

以上是通过模型计算出的生存概率，为了画出校准曲线，我们还需要真实的生存概率，怎么计算？通过K-M法。

如何通过K-M法计算真实的生存概率？

首先，分组。根据什么东西分组？就根据我们通过模型计算的概率分组。

为了保证每组都有人，我们就分为4组(可以多试几次)：

# 根据这几个点进行切分
cuts <- unique(quantile(c(0, 1, train_p), seq(0, 1, length = 5), na.rm = TRUE))
cuts

0.0000000 0.8097154 0.8526484 0.8892975 1.0000000

然后计算K-M法估计的生存概率：

suppressMessages(library(rms))

km_surv <- groupkm(train_p,
                   Srv = Surv(train_df$time,train_df$status),
                   u = 100,
                   cuts = cuts)

km_surv

             x  n events        KM    std.err
[1,] 0.7597781 28     24 0.7500000 0.10910895
[2,] 0.8291634 29     20 0.8620690 0.07427814
[3,] 0.8695150 30     23 0.8333333 0.08164966
[4,] 0.9131132 29     13 0.9310345 0.05053987

其中x就是模型预测的概率，KM是K-M法估计的真实概率，所以就可以画图了：

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'
     )
lines(km_surv[,1], km_surv[,4])

# 计算误差线范围
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'
       )
# 添加对角线
abline(a = 0,b = 1,col='grey')

这个图就基本和上面一样了，和开头的文献里的图一模一样！

测试集的校准曲线

方法一模一样，很简单！

测试集100天的校准曲线，首先也是计算概率：

set.seed(123)

test_p <- c((summary(survfit(coxph_fit, newdata=test_df), times=100)$surv))

然后分组：

cuts_t <- unique(quantile(c(0, 1, test_p), seq(0, 1, length = 5), na.rm = TRUE))
cuts_t

0.0000000 0.7887829 0.8588841 0.8829986 1.0000000

然后计算K-M法估计的生存概率：

suppressMessages(library(rms))

km_surv_t <- groupkm(test_p,
                     Srv = Surv(test_df$time,test_df$status),
                     u = 100,
                     cuts = cuts_t)

km_surv_t

             x  n events        KM   std.err
[1,] 0.7379403 12     10 0.6666667 0.2041241
[2,] 0.8263202 13     12 1.0000000 0.0000000
[3,] 0.8693478 13     11 0.8461538 0.1182625
[4,] 0.9073257 13      7 1.0000000 0.0000000

然后就可以画图了：

plot(km_surv_t[,1], km_surv_t[,4],
     xlim=c(0,1),ylim=c(0,1),
     xlab = 'Predicted 100-day Survival Probability',
     ylab = 'Observed 100-day Survival Probability'
     )
lines(km_surv_t[,1], km_surv_t[,4])

# 计算误差线范围
errl <- ifelse(km_surv_t[,"KM"] == 0, 0,  
               km_surv_t[,"KM"] * exp(1.959964 * (-km_surv_t[,"std.err"])))
errh <- ifelse(km_surv_t[,"KM"] == 0, 0, 
               pmin(1, km_surv_t[,"KM"] * exp(1.959964 * km_surv_t[,"std.err"])))
# 添加误差线
errbar(x = km_surv_t[,"x"],
       y = km_surv_t[,"KM"],
       yminus = errl,yplus = errh,
       add = T,
       pch=16,cex=1,
       asp=1,xaxs='i',yaxs='i'
       )
# 添加对角线
abline(a = 0,b = 1,col='grey')

这个就是测试集的校准曲线了！

两张图拼在一起就是大家文献中常见的图了！后期我会把这些函数打包，然后变为长数据，支持使用ggplot2画图！