歡迎關(guān)注”生信修煉手冊(cè)”!

在之前關(guān)于列線圖的文章中，我們介紹了利用列線圖來(lái)可視化預(yù)后模型，同時(shí)也提到了模型性能的幾種評(píng)估方式，校準(zhǔn)度以及校準(zhǔn)曲線就是其中一種方式。

校準(zhǔn)度，用來(lái)描述一個(gè)模型預(yù)測(cè)個(gè)體發(fā)生臨床結(jié)局的概率的準(zhǔn)確性。在實(shí)際應(yīng)用中，通常用校準(zhǔn)曲線來(lái)表征。校準(zhǔn)曲線展示了模型預(yù)測(cè)值與實(shí)際值之間的偏差，一個(gè)典型的校準(zhǔn)曲線示例如下

橫軸表示模型預(yù)測(cè)的不同臨床結(jié)局概率，縱軸表示實(shí)際觀察到的患者的臨床結(jié)局的概率，用中位數(shù)加均值的errorbar 形式表征，并繪制了一條斜率為1的理想曲線作為參照，實(shí)際曲線越接近理想曲線，表明模型預(yù)測(cè)結(jié)果與實(shí)際結(jié)果的偏差越小，模型效果高好。

在數(shù)據(jù)分析過(guò)程中，我們可以通過(guò)rms包中的calibrate函數(shù)來(lái)創(chuàng)建校準(zhǔn)曲線，首先來(lái)運(yùn)行下官方示例

> 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)

效果圖如下

參數(shù)u指定了我們想要分析的時(shí)間節(jié)點(diǎn)，m指定了樣本分組個(gè)數(shù)，該參數(shù)決定了圖中errorbar的個(gè)數(shù)，示例數(shù)據(jù)有200個(gè)樣本，m取50時(shí)，group的個(gè)數(shù)為4。該函數(shù)通過(guò)有放回的抽樣方法對(duì)模型效能進(jìn)行評(píng)估，利用函數(shù)返回值可以查看具體的繪圖數(shù)據(jù)，示例如下

> cal calibrate.cph(fit?=?f,?cmethod?=?"KM",?u?=?1.5,?m?=?50,?B?=?20) n=200??B=20??u=1.5?Dayindex.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 20mean.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對(duì)應(yīng)的x軸坐標(biāo)，KM.corrected列表示圖中黑色原形散點(diǎn)的縱坐標(biāo)，星形散點(diǎn)的縱坐標(biāo)為KM列，errobar的上下區(qū)間則通過(guò)如下公式計(jì)算

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列的數(shù)據(jù)進(jìn)行計(jì)算，?我們可以提取其中的數(shù)據(jù)，自己來(lái)畫(huà)圖，?代碼如下

> 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"]

效果圖如下

可以看到和直接用函數(shù)繪制出來(lái)的圖是完全一致的，掌握這個(gè)用法之后，我們就可以實(shí)現(xiàn)文獻(xiàn)中所示的個(gè)性化校準(zhǔn)曲線，比如下圖

只需要提取4個(gè)時(shí)間點(diǎn)的校準(zhǔn)曲線數(shù)據(jù)，然后自己繪圖賦予不同顏色即可實(shí)現(xiàn)。

·end·

—如果喜歡，快分享給你的朋友們吧—

原創(chuàng)不易，歡迎收藏，點(diǎn)贊，轉(zhuǎn)發(fā)！生信知識(shí)浩瀚如海，在生信學(xué)習(xí)的道路上，讓我們一起并肩作戰(zhàn)！

本公眾號(hào)深耕耘生信領(lǐng)域多年，具有豐富的數(shù)據(jù)分析經(jīng)驗(yàn)，致力于提供真正有價(jià)值的數(shù)據(jù)分析服務(wù)，擅長(zhǎng)個(gè)性化分析，歡迎有需要的老師和同學(xué)前來(lái)咨詢。

? 更多精彩

• KEGG數(shù)據(jù)庫(kù)，除了pathway你還知道哪些

• 全網(wǎng)最完整的circos中文教程

• DNA甲基化數(shù)據(jù)分析專題

• 突變檢測(cè)數(shù)據(jù)分析專題

• mRNA數(shù)據(jù)分析專題

• lncRNA數(shù)據(jù)分析專題

• circRNA數(shù)據(jù)分析專題

• miRNA數(shù)據(jù)分析專題

• 單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析專題

• chip_seq數(shù)據(jù)分析專題

• Hi-C數(shù)據(jù)分析專題

• HLA數(shù)據(jù)分析專題

• TCGA腫瘤數(shù)據(jù)分析專題

• 基因組組裝數(shù)據(jù)分析專題

• CNV數(shù)據(jù)分析專題

• GWAS數(shù)據(jù)分析專題

• 機(jī)器學(xué)習(xí)專題

• 2018年推文合集

• 2019年推文合集

• 2020推文合集

? 寫(xiě)在最后

轉(zhuǎn)發(fā)本文至朋友圈，后臺(tái)私信截圖即可加入生信交流群，和小伙伴一起學(xué)習(xí)交流。

掃描下方二維碼，關(guān)注我們，解鎖更多精彩內(nèi)容！

一個(gè)只分享干貨的

生信公眾號(hào)

總結(jié)

以上是生活随笔為你收集整理的校准曲线的绘制的小技巧的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。