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

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

之前给大家介绍了逻辑回归(logistic)测试集的校准曲线画法，给大家介绍了6种方法！

• logistic回归测试集校准曲线的6种实现方法

不过Cox回归测试集的校准曲线就没有那么简单了。

目前最简单的COX回归测试集的实现方法是通过riskRegression包。

本期目录：

使用 survival包的自带的 lung数据集做演示。

加载数据

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

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

dim(lung)
## [1] 167  10
str(lung)
## 'data.frame':	167 obs. of  10 variables:
##  $ inst     : num  3 5 12 7 11 1 7 6 12 22 ...
##  $ time     : num  455 210 1022 310 361 ...
##  $ status   : num  1 1 0 1 1 1 1 1 1 1 ...
##  $ age      : num  68 57 74 68 71 53 61 57 57 70 ...
##  $ sex      : num  1 1 1 2 2 1 1 1 1 1 ...
##  $ ph.ecog  : num  0 1 1 2 2 1 2 1 1 1 ...
##  $ ph.karno : num  90 90 50 70 60 70 70 80 80 90 ...
##  $ pat.karno: num  90 60 80 60 80 80 70 80 70 100 ...
##  $ meal.cal : num  1225 1150 513 384 538 ...
##  $ wt.loss  : num  15 11 0 10 1 16 34 27 60 -5 ...
##  - attr(*, "na.action")= 'omit' Named int [1:61] 1 3 5 12 13 14 16 20 23 25 ...
##   ..- attr(*, "names")= chr [1:61] "1" "3" "5" "12" ...

数据分割

把数据随机划分为训练集、测试集，划分比例为7:3

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

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

dim(train_df)
## [1] 116  10
dim(test_df)
## [1] 51 10

测试集校准曲线方法1

使用riskRegression轻松实现，训练集的校准曲线就不画了，直接画测试集的校准曲线。

library(riskRegression)
## riskRegression version 2022.03.22

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

# 画出测试集的校准曲线
set.seed(1)
cox_fit_t <- Score(list("fit1" = coxph_fit),
               formula = Surv(time, status) ~ 1,
               data = test_df, # 这里写测试集即可
               plots = "calibration",
               conf.int = T,
               B = 500, # 重抽样500次
               M = 50,
               times=c(100) # 时间
               )

# 画图即可
plotCalibration(cox_fit_t,
                cens.method="local",
                xlab = "Predicted Risk",
                ylab = "Observerd RISK")

这个就是测试集的校准曲线了，非常简单，而且支持返回数据自己用ggplot2画，你想要的样子它都可以有！ 可参考推文：生存资料校准曲线的绘制

测试集校准曲线方法2

回归建模神包rms当然也是支持COX回归测试集校准曲线的，不过隐藏的比较深。

在上次介绍logistic回归测试集的校准曲线时，我们介绍过val.prob函数，这里给大家介绍下val.surv函数。

首先还是要使用使用cph函数建立cox回归模型，这里我们指定时间为100天：

suppressMessages(library(rms))

# 一定要记得打包数据
dd <- datadist(train_df)
options(datadist = "dd")
units(train_df$time) <- "days"

cph_fit <- cph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno + pat.karno,
               data = train_df,
               x = T, y = T,
               time.inc = 100
               )

然后计算模型在测试集中的各种指标：

v <- val.surv(fit= cph_fit, # 模型
         newdata = test_df, # 测试集
         u=100, # 时间
         #evaluate = 10,
         S=Surv(test_df$time,test_df$status) # 测试集的生存对象
         )
v
## 
## Validation of Predicted Survival at Time= 100 	n= 51 , events= 40 
## 
## hare fit:
## 
## dim A/D   loglik       AIC        penalty 
##                                 min    max 
##   1 Add   -284.61    573.15   11.98     Inf
##   2 Add   -278.62    565.11    7.05   11.98
##   3 Add   -275.81    563.41      NA      NA
##   4 Add   -271.57    558.87    3.15    7.05
##   5 Add   -270.00    559.65    0.00    3.15
## 
## the present optimal number of dimensions is 4.
## penalty(AIC) was the default: BIC=log(samplesize): log(51)=3.93
## 
##   dim1           dim2           beta        SE         Wald
## Constant                            -5.6          1   -5.65
## Time   3.6e+02                      0.01     0.0052    1.95
## Co-1  linear                      -0.091       0.52   -0.17
## Time   3.6e+02 Co-1  linear       0.0094     0.0033    2.89
## 
## Function used to transform predictions:
## function (p)  log(-log(p))
## 
## Mean absolute error in predicted probabilities: 0.0894 
## 0.9 Quantile of absolute errors               : 0.1025

之后通过调用plot实现测试集的校准曲线，可以看到这个校准曲线并不是通过分箱的方法获得的，这里使用了一种hare方法(通过polspline包实现，感兴趣的可以自己看看)，且在测试集中没得选，只能用这种方法(训练集的校准曲线可通过method参数选择)。

但是这里也要注意，仔细阅读Regression Modeling Strategies 就会发现，val.surv计算的是外部验证集的各种指标，画出来的图也是外部验证集的，这一点在logistic回归测试集校准曲线的6种实现方法中也说过。

plot(v, 
     scat1d.opts=list(nhistSpike=200, side=3)
     )

这里由于数据集样本量太少，结果并不好看，大家可以使用自己的数据(样本量尽量大一点)尝试一下。

2种方法，我还是比较推荐riskRegression，因为简单又强大。

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

目前好像并没有包可以直接实现，不过也不是什么难题，下次介绍。

如果大家在文献中发现一些比较有意思的校准曲线图，欢迎通过后台、微信、QQ群等方法发给我，我会复现后推送。