校准曲线的绘制的小技巧
欢迎关注”生信修炼手册”!
在之前关于列线图的文章中，我们介绍了利用列线图来可视化预后模型，同时也提到了模型性能的几种评估方式，校准度以及校准曲线就是其中一种方式。
校准度，用来描述一个模型预测个体发生临床结局的概率的准确性。在实际应用中，通常用校准曲线来表征。校准曲线展示了模型预测值与实际值之间的偏差，一个典型的校准曲线示例如下
横轴表示模型预测的不同临床结局概率，纵轴表示实际观察到的患者的临床结局的概率，用中位数加均值的errorbar 形式表征，并绘制了一条斜率为1的理想曲线作为参照，实际曲线越接近理想曲线，表明模型预测结果与实际结果的偏差越小，模型效果高好。
在数据分析过程中，我们可以通过rms包中的calibrate函数来创建校准曲线，首先来运行下官方示例
> set.seed(1)
> n <- 200
> d.time <- rexp(n)
> x1 <- runif(n)
> x2 <- factor(sample(c('a', 'b', 'c'), n, TRUE))
> f <- cph(Surv(d.time) ~ pol(x1,2) * x2, x=TRUE, y=TRUE, surv=TRUE, time.inc=1.5)
> cal <- calibrate(f, u=1.5, cmethod='KM', m=50, B=20)
> plot(cal)效果图如下
参数u指定了我们想要分析的时间节点，m指定了样本分组个数，该参数决定了图中errorbar的个数，示例数据有200个样本，m取50时，group的个数为4。该函数通过有放回的抽样方法对模型效能进行评估，利用函数返回值可以查看具体的绘图数据，示例如下
> cal
calibrate.cph(fit = f, cmethod = "KM", u = 1.5, m = 50, B = 20)
n=200 B=20 u=1.5 Day
index.orig training test mean.optimism mean.corrected n
[1,] -0.02180909 -0.006492867 0.053098128 -0.05959099 0.03778191 20
[2,] 0.01161824 0.013463692 0.031802035 -0.01833834 0.02995658 20
[3,] 0.07007320 -0.064043654 -0.007650977 -0.05639268 0.12646588 14
[4,] -0.07103626 -0.015150576 -0.055302350 0.04015177 -0.11118804 20
mean.predicted KM KM.corrected std.err
[1,] 0.1418091 0.12 0.1795910 0.3829708
[2,] 0.1883818 0.20 0.2183383 0.2828427
[3,] 0.2299268 0.30 0.3563927 0.2160247
[4,] 0.3110363 0.24 0.1998482 0.2516611其中，mean.predicted列代表图中4处errorbar对应的x轴坐标，KM.corrected列表示图中黑色原形散点的纵坐标，星形散点的纵坐标为KM列，errobar的上下区间则通过如下公式计算
cal <- x[,"KM"]
se <- x[,"std.err"]
ciupper <- function(surv, d) ifelse(surv==0, 0, pmin(1, surv*exp(d)))
cilower <- function(surv, d) ifelse(surv==0, 0, surv*exp(-d))
cilower(cal, 1.959964*se)
ciupper(cal, 1.959964*se)利用KM列和std.err列的数据进行计算， 我们可以提取其中的数据，自己来画图， 代码如下
> x <- cal
> plot(x = x[,"mean.predicted"], y = x[,"KM"], pch = 20, xlab = "", ylab = "")
> errbar(x[,"mean.predicted"], x[,"KM"] , cilower(x[,"KM"], 1.959964 * x[,"std.err"]), ciupper(x[,"KM"], 1.959964 * x[,"std.err"]))
> points(x = x[,"mean.predicted"], y = x[,"KM.corrected"], pch = 4)
> lines(x = x[,"mean.predicted"], y = x[,"KM"])
> plot(x = x[,"mean.predicted"], y = x[,"KM"], pch = 20, xlab = "", ylab = "")
> errbar(x[,"mean.predicted"], x[,"KM"] , cilower(x[,"KM"], 1.959964 * x[,"std.err"]), ciupper(x[,"KM"], 1.959964 * x[,"std.err"]), xlab = "", ylab = "")
> points(x = x[,"mean.predicted"], y = x[,"KM.corrected"], pch = 4)
> lines(x = x[,"mean.predicted"], y = x[,"KM"]效果图如下
可以看到和直接用函数绘制出来的图是完全一致的，掌握这个用法之后，我们就可以实现文献中所示的个性化校准曲线，比如下图
只需要提取4个时间点的校准曲线数据，然后自己绘图赋予不同颜色即可实现。
·end·
—如果喜欢，快分享给你的朋友们吧—
原创不易，欢迎收藏，点赞，转发！生信知识浩瀚如海，在生信学习的道路上，让我们一起并肩作战！
本公众号深耕耘生信领域多年，具有丰富的数据分析经验，致力于提供真正有价值的数据分析服务，擅长个性化分析，欢迎有需要的老师和同学前来咨询。
更多精彩
- KEGG数据库，除了pathway你还知道哪些
- 全网最完整的circos中文教程
- DNA甲基化数据分析专题
- 突变检测数据分析专题
- mRNA数据分析专题
- lncRNA数据分析专题
- circRNA数据分析专题
- miRNA数据分析专题
- 单细胞转录组数据分析专题
- chip_seq数据分析专题
- Hi-C数据分析专题
- HLA数据分析专题
- TCGA肿瘤数据分析专题
- 基因组组装数据分析专题
- CNV数据分析专题
- GWAS数据分析专题
- 机器学习专题
- 2018年推文合集
- 2019年推文合集
- 2020推文合集
写在最后
转发本文至朋友圈，后台私信截图即可加入生信交流群，和小伙伴一起学习交流。
扫描下方二维码，关注我们，解锁更多精彩内容！
一个只分享干货的
生信公众号