logistic校准曲线（测试集）的6种实现方法

这个方法是基于。

医学和生信笔记

3707人浏览 · 2023-05-06 21:37:05

医学和生信笔记 · 2023-05-06 21:37:05 发布

关于临床预测模型的基础知识，小编之前已经写过非常详细的教程，包括了临床预测模型的定义、常用评价方法、列线图、ROC曲线、IDI、NRI、校准曲线、决策曲线等。

全都是免费获取的代码和数据：R语言临床预测模型合集

以上合集包括了临床预测模型绝大多数的内容，内容肯定比得上一个几千块的培训班！用以上教程完成一篇SCI绝对不是问题！

但是基础版合集留了几个问题尚未解决，主要集中在机器学习算法在临床预测模型中的使用以及临床预测模型细节解读和实现，并没有做详细的教程。所以后面的临床预测模型系列文章，会把重心放在以上问题中。

今天给大家展示的是测试集(或者叫验证集)的校准曲线如何实现(其实已经介绍过，不过没有单独说，有粉丝一直在后台问)。

本期目录：

文章目录

准备数据

数据来自于这篇推文：二分类资料校准曲线的绘制，数据获取方法也在上面的推文中给出了。

lowbirth <- read.csv("../000预测模型/lowbirth.csv")
lowbirth$black = ifelse(lowbirth$race == "black",1,0)
lowbirth$white = ifelse(lowbirth$race == "white",1,0)
lowbirth$other = ifelse(lowbirth$race %in% c("native American","oriental"),1,0)
lowbirth$delivery = factor(lowbirth$delivery)
lowbirth$sex <- factor(lowbirth$sex)
lowbirth$race <- NULL
str(lowbirth)
## 'data.frame':	565 obs. of  12 variables:
##  $ birth   : num  81.5 81.6 81.6 81.6 81.6 ...
##  $ lowph   : num  7.25 7.06 7.25 6.97 7.32 ...
##  $ pltct   : int  244 114 182 54 282 153 229 182 361 378 ...
##  $ bwt     : int  1370 620 1480 925 1255 1350 1310 1110 1180 970 ...
##  $ delivery: Factor w/ 2 levels "abdominal","vaginal": 1 2 2 1 2 1 2 2 1 2 ...
##  $ apg1    : int  7 1 8 5 9 4 6 6 6 2 ...
##  $ vent    : int  0 1 0 1 0 0 1 0 0 1 ...
##  $ sex     : Factor w/ 2 levels "female","male": 1 1 2 1 1 1 2 2 2 1 ...
##  $ dead    : int  0 1 0 1 0 0 0 0 0 1 ...
##  $ black   : num  0 1 1 1 1 1 0 1 0 0 ...
##  $ white   : num  1 0 0 0 0 0 1 0 1 1 ...
##  $ other   : num  0 0 0 0 0 0 0 0 0 0 ...

数据分割

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

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

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

训练集的校准曲线

在之前的推文中这种二分类资料训练集的校准曲线给大家介绍了非常多的方法：

二分类资料校准曲线的绘制

这里我们直接使用rms包实现，已在上面的推文中详细介绍过了，这里就不多解释了。

library(rms)

dd <- datadist(train_df)
options(datadist="dd")

fit1 <- lrm(dead ~ birth + lowph + pltct + bwt + vent + black + white,
            data = train_df,x=T,y=T)

cal1 <- calibrate(fit1, method='boot', B=100)

plot(cal1,
     xlim = c(0,1),
     ylim = c(0,1),
     xlab = "Prediced Probability",
     ylab = "Observed Probability"
     )

plot of chunk unnamed-chunk-3

## 
## n=395   Mean absolute error=0.01   Mean squared error=0.00016
## 0.9 Quantile of absolute error=0.021

测试集校准曲线方法1

测试集的校准曲线对于logistic回归很简单，任何可以计算概率的算法都可以轻松画出训练集、测试集的校准曲线，无非就是计算实际概率和预测概率而已。

二分类资料测试集的校准曲线在之前的推文中也做过很多次介绍，比如：

上面是3种实现方法，其实本质是一样的，前2种是手动计算，最后一种省去了自己计算的步骤，直接给你图形，并且完美继承yardstick的用法。

这里再给大家介绍3种方法，加上上面介绍的方法，logistic测试集的校准曲线一共给大家介绍了6种方法！

这个方法是基于rms包的。

# 首先获取测试集的预测结果
phat <- predict(fit1, test_df, type = 'fitted')

# 直接使用val.prob即可实现
val.prob(phat, test_df$dead,statloc = F,cex = 1)

plot of chunk unnamed-chunk-4

##         Dxy     C (ROC)          R2           D    D:Chi-sq         D:p 
##  0.74384874  0.87192437  0.42860880  0.28181902 48.90923309          NA 
##           U    U:Chi-sq         U:p           Q       Brier   Intercept 
## -0.01085663  0.15437207  0.92571762  0.29267565  0.08935692  0.12059928 
##       Slope        Emax         E90        Eavg         S:z         S:p 
##  1.08566597  0.29112563  0.07879941  0.03183303 -0.54088957  0.58858370

这个结果，你可以把它当做测试集的校准曲线用，但其实val.prob函数的真正作用是实现外部验证数据(external validation)的校准曲线，这一点在Harrell大神写的书：Regression Modeling Strategies 中写的很清楚，或者你可以看函数的帮助文档。

所以这个方法没有给你重抽样的选择，因为作者认为外部验证是对模型最后的检验，不需要重抽样。

你可能在文献看见过训练集和测试集的校准曲线都是上面那张图的样式，类似下面这张图展示的，训练集和测试集一样的图，实现方法也很简单。

# 获取训练集的预测结果
phat_train <- predict(fit1, train_df, type = 'fitted')

# 直接使用val.prob即可实现
val.prob(phat_train, train_df$dead,cex = 1)

          Dxy       C (ROC)            R2             D      D:Chi-sq           D:p 
 7.595104e-01  8.797552e-01  4.274140e-01  2.924665e-01  1.165243e+02            NA 
            U      U:Chi-sq           U:p             Q         Brier     Intercept 
-5.063291e-03  2.842171e-14  1.000000e+00  2.975298e-01  9.802204e-02 -8.039770e-10 
        Slope          Emax           E90          Eavg           S:z           S:p 
 1.000000e+00  3.476792e-02  3.030223e-02  1.012113e-02 -1.990865e-01  8.421951e-01

上面这张图可以用作训练集的校准曲线。

测试集校准曲线方法2

如果你非要对测试集的校准曲线进行重抽样，其实也很简单(除了rms还有很多手段可实现)。

这里还是用rms包实现。

二分类资料的校准曲线就是计算下实际概率和预测概率就好了，基于这个原理，我们可以自己实现，方法如下：

# 首先也是获取测试集的预测值
phat <- predict(fit1, test_df)
test_df$phat <- phat # 添加到测试集中

# 以预测值为自变量，结果变量为因变量，在测试集中建立逻辑回归
fit2 <- lrm(dead ~ phat, data = test_df,x=T,y=T)

# 对这个逻辑回归画校准曲线即可
cal2 <- calibrate(fit2, method='boot', B=100)

plot(cal2,
     xlim = c(0,1),
     ylim = c(0,1),
     xlab = "Prediced Probability",
     ylab = "Observed Probability"
     )

plot of chunk unnamed-chunk-5

## 
## n=170   Mean absolute error=0.031   Mean squared error=0.00183
## 0.9 Quantile of absolute error=0.078

这个图就是测试集的校准曲线，并且使用了bootstrap方法进行了重抽样。

可以看到其实两张图是一样的，唯一不同是我们手动实现的方法多了重抽样100次的矫正曲线，其余就都是一样的了！

测试集校准曲线方法3

使用riskRegression包。这是我推荐的方法，这个包真的好用！你可能还听说过另一个包：pec，riskRegression正是pec的升级版，目前大部分功能都已经转移到riskRegression了。

# 训练集建立模型
fit <- glm(dead ~ birth + lowph + pltct + bwt + vent + black + white,
            data = train_df, family = binomial)

library(riskRegression)
## riskRegression version 2022.03.22

# 测试集的表现
score <- Score(list("fit"=fit),
               formula = dead ~ 1,
               data = test_df, # 这里选择测试集即可
               metrics = c("auc","brier"),
               summary = c("risks","IPA","riskQuantile","ibs"),
               plots = "calibration",
               null.model = T,
               conf.int = T,
               B = 100,
               M = 50
               )

# 画图即可
plotCalibration(score,col="tomato",
                xlab = "Predicted Risk",
                ylab = "Observerd RISK",
                bars = F)