用诺模图可视化你的模型

花花写于2020.6.2
诺模图出镜率很高，用于多因素cox或logstic模型的可视化展示。

0.输入数据

需要病人临床信息和生存信息表格，如下

rm(list = ls())
load("ph.Rdata")
head(ph)

##                         sample_id gender age T N  M stage event time
## TCGA-MP-A4T4-01A TCGA-MP-A4T4-01A female  68 2 1  0     2     1 2617
## TCGA-05-4250-01A TCGA-05-4250-01A female  79 3 1  0     3     1  121
## TCGA-64-5774-01A TCGA-64-5774-01A   male  60 2 0  0     1     0 2676
## TCGA-97-7937-01A TCGA-97-7937-01A   male  65 2 0 NA     1     0  564
## TCGA-97-A4M6-01A TCGA-97-A4M6-01A female  45 1 0  0     1     0  568
## TCGA-44-A4SU-01A TCGA-44-A4SU-01A female  67 1 0 NA     1     1  409
##                            fp
## TCGA-MP-A4T4-01A -0.005061199
## TCGA-05-4250-01A  0.286351686
## TCGA-64-5774-01A  0.389468347
## TCGA-97-7937-01A  0.098764940
## TCGA-97-A4M6-01A -0.325363900
## TCGA-44-A4SU-01A -0.183674628

1.简约诺模图

使用rms包里的cph函数建模，nomogram函数画图

library(rms)
dd<-datadist(ph)
options(datadist="dd")
mod <- cph(formula =  as.formula(paste("Surv(time, event) ~ ",paste(colnames(ph)[2:7],collapse = "+"))),
           data=ph,x=T,y=T,surv = T)

surv<-Survival(mod) 
surv3<-function(x) surv(1095,x)
surv5<-function(x) surv(1825,x)

x<-nomogram(mod,
            fun = list(surv3,surv5),
            funlabel = c('3-year survival Probability',
                         '5-year survival Probability'))

plot(x)

2.美化版本

regplot画出的图和上面的简约版意义一致，只是可修改的细节多一些。

library(regplot)
library(survival)

mod2 <- coxph(formula =  as.formula(paste("Surv(time, event) ~ ",paste(colnames(ph)[2:7],collapse = "+"))),data=ph)

2.1 简约风

regplot(mod2,
        failtime = c(1095,1825), 
        plots = c("no plot","no plot"),
        points = T,
        prfail = T)

2.2 在图上标出某个病人

regplot(mod2,
        observation=ph[1,], 
        obscol = "#326db1",
        failtime = c(1095,1825), 
        plots = c("no plot","no plot"),
        points = T,
        prfail = T)

2.3. 增加表示分布情况的图形

regplot(mod2,
        observation=ph[1,], 
        failtime = c(1095,1825), 
        plots = c("bars","boxes"),
        points = T,
        prfail = T)

帮助文档中对plots参数的解读：The plots parameter specifies initial plot types. It is length 2. The first item specifies a plot type for non-factor variables as one of: “no plot”, “density”, “boxes”, “spikes”, “ecdf”, “bars”, “boxplot”, “violin” or “bean”. The second item, is for factors and is one of: “no plot”, “boxes”, “bars” or “spikes”.

其他图形也可以试试哦。我觉得简单的就挺好的，已经可以说明问题了。

3.校准曲线

校准曲线的作用是展示模型质量，曲线越贴近0-1对角线越好。

3.1 建模和完成计算

cph的time.inc参数和calibrate的u参数，后面都是写天数，三年就是1095天，5年就是1825天

f3 <- cph(formula =  as.formula(paste("Surv(time, event) ~ ",paste(colnames(ph)[2:7],collapse = "+"))),
           data=ph,x=T,y=T,surv = T, time.inc=1095)
cal3 <- calibrate(f3, cmethod="KM", method="boot", u=1095, m=50, B=1000)

## Using Cox survival estimates at 1095 Days

f5 <- cph(formula =  as.formula(paste("Surv(time, event) ~ ",paste(colnames(ph)[2:7],collapse = "+"))),
           data=ph,x=T,y=T,surv = T,  time.inc=1825)
cal5 <- calibrate(f5, cmethod="KM", method="boot", u=1825, m=50, B=1000)

## Using Cox survival estimates at 1825 Days

3.2 画校准曲线图

plot(cal5)就可以画，但出来的图一言难尽，主要是参考线画的并不是0-1，所以需要下面的代码手动去画咯。

plot(cal3,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(cal3[,c('mean.predicted',"KM")],
      type = 'b', lwd = 1, col = c("#2166AC"), pch = 16)
mtext("")

plot(cal5,lwd = 2,lty = 0,errbar.col = c("#B2182B"),
     xlim = c(0,1),ylim= c(0,1),col = c("#B2182B"),add = T)
lines(cal5[,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("3-year","5-year"), #图例文字
       col =c("#2166AC","#B2182B"), #图例线的颜色，与文字对应
       lwd = 2,#图例中线的粗细
       cex = 1.2,#图例字体大小
       bty = "n")#不显示图例边框

参考：https://blog.csdn.net/weixin_30592887/article/details/112865824 https://www.jianshu.com/p/9085e4e13843 https://www.dxy.cn/bbs/newweb/pc/post/33323120