使用 ggplot2 自定义美化你的校准曲线

1引言

来自粉丝(弱小，可怜，又被嫌弃的路易李四)投稿,感谢粉丝分享。(如果你也想投稿分享知识,欢迎联系老俊俊!)

2正文

在两个月前，俊俊的 QQ 群里有人分享一个推送，想问一个校准曲线的结果如何绘制。当时比较忙，拖了一个星期之后才想着复现。现在分享代码给大家。

这个校准曲线还是比较漂亮的，让人感觉眼前一亮。想要复原的话，得要先理清绘图思路。首先，这是一个双 Y 轴图，主图是点线图，也就是图中的红色部分。然后就是直方图部分，画的是啥？那当然得要看 x 轴是啥，看到这里，相信对模型建立和校准曲线有些了解的冰雪聪明的同学们应该不难猜出，其实就是预测概率。不同颜色分别对应 0 和 1，也就是正例和反例。这个直方图就是模型对正例和反例们的预测概率。

所以，绘图思路明确之后，开始操作：

使用基于 rms 包的结果进行校准曲线的美化：首先加载包及数据:

library(tidyverse)
library(rms)
library(survival)
library(ggsci)
library(ggthemes)

head(colon)
colon <- na.omit(colon)

（为了便于同学们使用，我坚持使用 R 的内置数据集，因为我不想让大家后台留言甚至付费去换取示例数据，麻烦。我当然知道使用生存数据建分类模型不合理，这点水平我还是有的，不用抬杠。）

封装数据，建立模型及校准曲线对象：

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

fit1 <- lrm(status ~ rx + sex + age + obstruct + adhere + nodes + etype,
            data = colon,x=T,y=T)

cal1 <- calibrate(fit1, method='boot', B=500,smoother = 'none')

然后，就是最常规的绘制校准曲线的方式：

pdf('jzqx_yuanshi.pdf')
plot(cal1,
     xlim = c(0,1),
     ylim = c(0,1),
     xlab = "Prediced Probability",
     ylab = "Observed Probability",
     cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=0.8,
     #subtitles = FALSE,
     legend = FALSE
)
lines(cal1[,c("predy","calibrated.corrected")],
      type = 'l', #连线的类型，可以是"p","b","o"
      lwd = 3, #连线的粗细
      pch = 16, #点的形状，可以是0-20
      col = "#2166AC") #连线的颜色
lines(cal1[,c("predy","calibrated.orig")],type="l",pch=16,lwd=3,col="tomato")
abline(0,1,
       lty = 2, #对角线为虚线
       lwd = 2, #对角线的粗细
       col = "#224444") #对角线的颜色
legend(0.6,0.2,
       c("Apparent","Bias-corrected","Ideal"),
       lty = c(2,1,1),
       lwd = c(2,3,3),
       col = c("black","#2166AC","tomato"),
       bty = "n"
)
dev.off()

解螺旋的代码就是这样绘制的，包括实际中见到的大多也是这样的，效果如下：

红色的为原始的校准曲线，蓝色的是经过 500 次自助法抽样后矫正了的校准曲线，因此，在下面复现中，使用自助法矫正的校准曲线。

下面，开始美化工作。使用 ggplot2 重绘图，首先，就要考虑数据从哪里进行提取。 在这里，自然是在校准曲线对象中进行提取（不然也没其他地方去提取数据）。

cal1

直接查看一下对象，输出的东西较多，主要就是两个：

1：一个 matrix：

这里储存的是曲线的绘制数据。

2：预测概率：

这里是待会需要绘制直方图的数据。

首先，绘制曲线部分。横坐标为 predy，纵坐标则为 calibrated.corrected。因此，首先提取出数据，保存为一个数据框：

df <- cal1[,c("predy","calibrated.corrected")]%>%
  as.data.frame()
head(df)

先尝试绘制一下校准曲线部分：

ggplot()+
  geom_line(data = df,aes(predy,calibrated.corrected),col='red3',linewidth = 1)+
  geom_point(data = df[c(seq(1,50,5),50),],      #将数据分为10段，一共11个点
             aes(predy,calibrated.corrected),col='red3',size = 2)

然后，绘制直方图部分，得要先提取数据：

df2 <- data.frame(prob = attr(cal1,"predicted"), #预测概率
                  group = colon$status) #从源数据中提取实际分类情况
df2$group <- as.factor(df2$group) #将分类变量设为因子型
head(df2)

在尝试绘制一下直方图部分：

#用这个看一下直方图分成多少份合适，切决定接下来左侧y轴大小
ggplot(data = df2,aes(x=prob,fill=group))+
  geom_histogram(binwidth = 1/20,alpha = 0.5)+
  scale_fill_aaas()
#分成20份后，坐标轴以400为上限比较合适。这个值需要根据情况调整。

在这一步需要确定将直方图分成多少份以及 Y 轴高度，根据需要设定。colon 数据集样本量大，因此接下来的绘图中设为 400。这一步很重要，决定了下一步 Y 轴高度使用多大的数值。

最后，将二者合二为一，画作双 Y 轴图。

ggplot()+
  geom_line(data = df,aes(predy,calibrated.corrected),col='red3',linewidth = 1)+
  geom_point(data = df[c(seq(1,50,5),50),],aes(predy,calibrated.corrected),col='red3',size = 2)+
  labs(x = "Prediced Probability",y = "Observed Probability")+
  scale_x_continuous(limits = c(0,1),
                     expand = c(0,0))+
  scale_y_continuous(limits = c(0,1),
                     expand = c(0,0),
                     sec.axis = sec_axis(~.*400,name = 'Samples', #一定要注意这里的400是Y轴高度！！！Samples是右侧坐标轴名
                                         breaks = c(50,100,150,200,250,300,350)))+ #加右边的y轴就是这样
  geom_histogram(data = df2,aes(x=prob,after_stat(count/400),fill=group),#一定要注意这里的400是Y轴高度！！！记得改
                 binwidth = 0.05,alpha = .5,show.legend = F)+ #关闭图例
  scale_fill_d3()+
  annotate("segment", x = 0, xend = 1, y = 0, yend = 1, #绘制对角线，以虚线绘制
           colour = "black",lty = 2)+
  annotate('text',x=0.3,y=0.9,label = 'Hosmer-Lemeshow: p =0.8',size = 5)+ #随便编的数值，根据实际数值添加
  theme_bw()

如果使用的话，一定要记得把两个“400”改为自己实际需要的 Y 轴高度。最后效果：

还可以根据实际需要再对主题、颜色等进行修改。

3结尾

路漫漫其修远兮,吾将上下而求索。

要使用ggplot2自定义绘制各种统计图形，提取模型中的绘图数据很重要，学习下。

本文转载自老俊俊的生信笔记。

关注下方公众号，分享更多更好玩的R语言知识。

点个在看，SCI马上发表。