R语言代码-绘制SCI发表级的nomogram及calibration图

需求描述

画出paper里的nomogram图和校准曲线


demo1.png
demo2.png

出自https://linkinghub.elsevier.com/retrieve/pii/S1470-2045(13)70491-1

画出新型nomogram


demo3.png

应用场景

列线图(nomogram,诺莫图)是在平面直角坐标系中用一簇互不相交的线段表示多个独立变量之间函数关系的图。

将Logistic回归或Cox回归的结果进行可视化呈现,给出其发病风险或比例风险。

输入数据

至少要有临床信息easy_input.csv,还可以加入FigureYa31lasso输出的lasso_output.txt

pbc<-read.table("easy_input.txt")

pbc$catbili <- cut(pbc$bili,breaks=c(-Inf, 2, 4, Inf),
                   labels=c("low","medium","high"))
pbc$died <- pbc$status==2

head(pbc)

经典版

#install.packages("rms")
library(rms)

dd<-datadist(pbc)
options(datadist="dd")
options(na.action="na.delete")
summary(pbc$time)
coxpbc<-cph(formula = Surv(time,died) ~  age + catbili + sex + copper + stage + trt ,data=pbc,x=T,y=T,surv = T,na.action=na.delete)  #,time.inc =2920

print(coxpbc)
surv<-Survival(coxpbc) 
surv3<-function(x) surv(1825,x)
surv4<-function(x) surv(2920,x)

x<-nomogram(coxpbc,fun = list(surv3,surv4),lp=T,
            funlabel = c('5-year survival Probability','8-year survival Probability'),
            maxscale = 100,fun.at = c(0.95,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1))

pdf("nomogram_classical.pdf",width = 12,height = 10)
plot(x, lplabel="Linear Predictor",
     xfrac=.35,varname.label=TRUE, varname.label.sep="=", ia.space=.2, 
     tck=NA, tcl=-0.20, lmgp=0.3,
     points.label='Points', total.points.label='Total Points',
     total.sep.page=FALSE, 
     cap.labels=FALSE,cex.var = 1.6,cex.axis = 1.05,lwd=5,
     label.every = 1,col.grid = gray(c(0.8, 0.95)))
dev.off()
nomogram0.png

输出公式

#print(x)
#install.packages("nomogramEx")
library(nomogramEx)
nomogramEx(nomo=x,np=2,digit = 9)
$RESULT
[1] "The equation of each variable as follows:"

[[2]]
[1] "points = 0 * age ^2 + 1.720070183 * age + -43.001754587"

[[3]]

[[4]]

[[5]]
[1] "points = 0 * copper ^3 + 0 * copper ^2 + 0.166666667 * copper + 0"

[[6]]

[[7]]

[[8]]
[1] "5-year survival Probability = 5.9e-08 * points ^3 + -5.1054e-05 * points ^2 + 0.008126192 * points + 0.583874866"

[[9]]
[1] "8-year survival Probability = 5.9e-08 * points ^3 + -4.3812e-05 * points ^2 + 0.004242289 * points + 0.835039968"

新型nomogram

regplot可以交互式调整nomogram,代码更简单,输出的图形更漂亮。

regplot那段代码粘贴到Console里,回车,就会在Plots窗口出现图。

鼠标点击栏目上方蓝色按钮,可以选择变量展示的类型:箱图、小提琴图、密度曲线图等;需要什么图,导出PDF进行AI修改即可。

鼠标点击Esc或按键盘上的Esc可以退出交互模式

#不喜欢默认的颜色,先设置几个颜色
mycol<-c("#A6CEE3","#1F78B4","#33adff","#2166AC")
names(mycol) = c("dencol","boxcocl","obscol","spkcol")
mycol<- as.list(mycol)

#install.packages("regplot","vioplot","sm","beanplot")
library(survival)
pbccox <- coxph(formula = Surv(time,died) ~  age + catbili + sex + 
                   copper + stage + trt , data = pbc)

library(regplot)

pdf("nomogram_new.pdf")
regplot(pbccox,
        #对观测2的六个指标在列线图上进行计分展示
        observation=pbc[2,], #也可以不展示
        
        #预测3年和5年的死亡风险,此处单位是day
        failtime = c(1095,1825), 
        prfail = TRUE, #cox回归中需要TRUE
        showP = T, #是否展示统计学差异
        droplines = F,#观测2示例计分是否画线
        colors = mycol, #用前面自己定义的颜色
        rank="decreasing", #根据统计学差异的显著性进行变量的排序
        interval="confidence") #展示观测的可信区间
dev.off()
nomogram.png

展示逻辑回归,支持"lm", "glm", "coxph", "survreg" "negbin"

## Plot a logistic regression, showing odds scale and confidence interval
pbcglm <- glm(died ~  age + catbili + sex + copper, family = "binomial", data=pbc )

regplot(pbcglm, 
        observation=pbc[1,], 
        odds=TRUE, 
        interval="confidence")
nomo.png

绘制calibration curve进行验证

采用和nomogram一样的变量进行多因素cox回归

5年

f5<-cph(formula = Surv(time,died) ~  age + catbili + sex + copper +stage + trt,data=pbc,x=T,y=T,surv = T,na.action=na.delete,time.inc = 1825) 

#参数m=50表示每组50个样本进行重复计算
cal5<-calibrate(f5, cmethod="KM", method="boot",u=1825,m=50,B=1000) 

pdf("calibration_5y.pdf",width = 8,height = 8)
plot(cal5,
     lwd = 2,#error bar的粗细
     lty = 1,#error bar的类型,可以是0-6
     errbar.col = c("#2166AC"),#error bar的颜色
     xlim = c(0,1),ylim= c(0,1),
     xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)",
     cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=0.6) #字的大小
lines(cal5[,c('mean.predicted',"KM")], 
      type = 'b', #连线的类型,可以是"p","b","o"
      lwd = 2, #连线的粗细
      pch = 16, #点的形状,可以是0-20
      col = c("#2166AC")) #连线的颜色
mtext("")
box(lwd = 1) #边框粗细
abline(0,1,lty = 3, #对角线为虚线
       lwd = 2, #对角线的粗细
       col = c("#224444")#对角线的颜色
       ) 
dev.off()
cali5.png

8年

f8<-cph(formula = Surv(time,died) ~  age + catbili + sex + copper +stage + trt,data=pbc,x=T,y=T,surv = T,na.action=na.delete,time.inc = 2920) 
cal8<-calibrate(f8, cmethod="KM", method="boot",u=2920,m=50,B=1000)

plot(cal8,
     lwd = 2,
     lty = 1,
     errbar.col = c("#B2182B"),
     xlim = c(0,1),ylim= c(0,1),
     xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)",
     col = c("#B2182B"),
     cex.lab=1.2,cex.axis=1, cex.main=1.2, cex.sub=0.6)
lines(cal8[,c('mean.predicted',"KM")],
      type= 'b',
      lwd = 2,
      col = c("#B2182B"),
      pch = 16)
mtext("")
box(lwd = 1)
abline(0,1,lty= 3,
       lwd = 2,
       col =c("#224444"))

同时展示两条curve

pdf("calibration_compare.pdf",width = 8,height = 8)
plot(cal5,lwd = 2,lty = 0,errbar.col = c("#2166AC"),
     bty = "l", #只画左边和下边框
     xlim = c(0,1),ylim= c(0,1),
     xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)",
     col = c("#2166AC"),
     cex.lab=1.2,cex.axis=1, cex.main=1.2, cex.sub=0.6)
lines(cal5[,c('mean.predicted',"KM")],
      type = 'b', lwd = 1, col = c("#2166AC"), pch = 16)
mtext("")

plot(cal8,lwd = 2,lty = 0,errbar.col = c("#B2182B"),
     xlim = c(0,1),ylim= c(0,1),col = c("#B2182B"),add = T)
lines(cal8[,c('mean.predicted',"KM")],
      type = 'b', lwd = 1, col = c("#B2182B"), pch = 16)

abline(0,1, lwd = 2, lty = 3, col = c("#224444"))

legend("topleft", #图例的位置
       legend = c("5-year","8-year"), #图例文字
       col =c("#2166AC","#B2182B"), #图例线的颜色,与文字对应
       lwd = 2,#图例中线的粗细
       cex = 1.2,#图例字体大小
       bty = "n")#不显示图例边框
dev.off()
cali.png
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,684评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,143评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,214评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,788评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,796评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,665评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,027评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,679评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 41,346评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,664评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,766评论 1 331
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,412评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,015评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,974评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,073评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,501评论 2 343