用诺模图可视化你的模型
花花写于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
