1. 使用
包绘制 nomogram图
1.1 基本用法
# 安装并加载所需的R包
# install.packages("rms")
library(rms)
# 加载数据
pbc <- pbc[pbc$status %in% c(0, 1), ]
head(pbc)
## id time status trt age sex ascites hepato spiders edema bili chol albumin copper alk.phos ast trig platelet protime stage
## 2 2 4500 0 1 56.44627 f 0 1 1 0.0 1.1 302 4.14 54 7394.8 113.52 88 221 10.6 3
## 5 5 1504 1 2 38.10541 f 0 1 1 0.0 3.4 279 3.53 143 671.0 113.15 72 136 10.9 3
## 7 7 1832 0 2 55.53457 f 0 1 0 0.0 1.0 322 4.09 52 824.0 60.45 213 204 9.7 3
## 13 13 3577 0 2 45.68925 f 0 0 0 0.0 0.7 281 3.85 40 1181.0 88.35 130 244 10.6 3
## 16 16 3672 0 2 40.44353 f 0 0 0 0.0 0.7 204 3.66 28 685.0 72.85 58 198 10.8 3
## 19 19 4232 0 1 49.56057 f 0 1 0 0.5 0.7 235 3.56 39 1881.0 93.00 123 209 11.0 3
ddist <- datadist(pbc)
options(datadist = 'ddist')
# 拟合-逻辑回归模型
f <- lrm(status ~ age + bili + copper + stage, data = pbc)
f # ★ 模型的C指数为0.864，预测能力良好
## Frequencies of Missing Values Due to Each Variable
## status age bili copper stage
## 0 0 0 71 2
##
## Logistic Regression Model
##
## lrm(formula = status ~ age + bili + copper + stage, data = pbc)
##
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 186 LR chi2 26.67 R2 0.277 C 0.864
## 0 167 d.f. 4 R2(4,186)0.115 Dxy 0.727
## 1 19 Pr(> chi2) <0.0001 R2(4,51.2)0.358 gamma 0.728
## max |deriv| 2e-06 Brier 0.080 tau-a 0.134
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -1.4859 1.5504 -0.96 0.3378
## age -0.0842 0.0288 -2.92 0.0035
## bili 0.2208 0.1146 1.93 0.0541
## copper 0.0064 0.0033 1.93 0.0539
## stage 0.6617 0.3577 1.85 0.0643
nom <- nomogram(f, fun = plogis, funlabel = "Risk of Death")
plot(nom)
rms-1
1.2 构建生存模型，预测患者的生存时间和生存概率
f <- psm(Surv(time,status) ~ age + bili + copper + stage, data = pbc, dist='lognormal')
# 计算中位生存时间
med <- Quantile(f)
# 构建生存概率
surv <- Survival(f)
nom <- nomogram(f,
fun = list(function(x) med(lp = x, q = 0.5),
function(x) surv(365, x),
function(x) surv(365 * 3, x),
function(x) surv(365 * 5, x)),
funlabel = c("Median Survival Time", # 添加轴的标签
"1-year Survival Probability",
"3-year Survival Probability",
"5-year Survival Probability"),
fun.at = list(c(10000, 20000, 40000, 90000), # 指定要在轴上标记的函数值
c(0.8, 0.90, 0.95, 0.99),
c(0.5, 0.7, 0.9),
c(0.2, 0.4, 0.6, 0.8, 0.9)),
lp = F) # 不显示“Linear Predictor”这个轴
plot(nom)
rms-2
2. 使用
包绘制 nomogram图
2.1 基本用法
# 安装并加载所需的R包
# install.packages("regplot")
library(regplot)
library(survival)
mod2 <- coxph(formula = as.formula(paste("Surv(time,status) ~ ", paste(colnames(pbc)[c(5, 11, 14, 20)],collapse = "+"))), data = pbc)
## Call:
## coxph(formula = as.formula(paste("Surv(time,status) ~ ", paste(colnames(pbc)[c(5,
## 11, 14, 20)], collapse = "+"))), data = pbc)
##
## coef exp(coef) se(coef) z p
## age -0.069487 0.932872 0.025029 -2.776 0.00550
## bili 0.219639 1.245627 0.073620 2.983 0.00285
## copper 0.005304 1.005318 0.001977 2.683 0.00729
## stage 0.746101 2.108763 0.324986 2.296 0.02169
##
## Likelihood ratio test=31.07 on 4 df, p=2.956e-06
## n= 186, number of events= 19
## (71 observations deleted due to missingness)
regplot(mod2,
failtime = c(1095, 1825), # 定义生存的两个概率标度
plots = c("violin","boxes"), # 连续性变量形状，可选"no plot" "density" "boxes" "ecdf" "bars" "boxplot" "violin" "bean" "spikes"；分类变量的形状，可选"no plot" "boxes" "bars" "spikes"
points = T, # 截距项显示为0-100
prfail = T)
regplot-1
2.2 其他参数
regplot(mod2,
observation=pbc[1,], #用哪行观测
obscol = "#326db1",
failtime = c(1095,1825),
plots = c("bars","boxes"),
droplines = T, # 是否画竖线
points = T,
title = "nomogram", # 更换标题
# odds = T, # 是否显示OR值
showP = T, # 是否显示变量的显著性标记（默认：T）
rank = "sd", # 根据sd给变量排序
# interval="confidence", # 展示可信区间
clickable = F, # 是否可以交互（默认：F）
prfail = T)
regplot-2
3. 使用包绘制校准曲线图
校准曲线图表示的是预测值和实际值的差距，作为预测模型的重要部分
# 安装并加载所需的R包
# install.packages("survival")
library("survival")
Datdata("lung")
head(lung)
## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 0
# 建模并完成计算
f300 <- cph(formula = as.formula(paste("Surv(time, status) ~ ",paste(colnames(lung)[4:7],collapse = "+"))),
data = lung, x = T,y = T,surv = T, time.inc = 300) # time.inc参数和calibrate的u参数后接天数
cal300 <- calibrate(f300, cmethod = "KM", method = "boot", u = 300, m = 75, B = 1000) # m，分组到平均包含 m 个受试者的区间；
f500 <- cph(formula = as.formula(paste("Surv(time, status) ~ ",paste(colnames(lung)[4:7],collapse = "+"))),
data = lung, x = T,y = T,surv = T, time.inc = 500)
cal500 <- calibrate(f500, cmethod="KM", method="boot", u = 500, m = 75, B = 1000)
# 绘制校准曲线图
plot(cal300, lwd = 2, # error bar的粗细
lty = 0, # error bar的类型，可以是0-6
errbar.col = c("#2166AC"), # error bar的颜色
# bty = "l", # 只画左边和下边框
xlim = c(0.3,0.7),ylim= c(0.3,0.8),
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(cal300[,c('mean.predicted',"KM")],
type = 'b', lwd = 1, col = c("#2166AC"), pch = 16)
mtext("")
plot(cal500, lwd = 2, lty = 0, errbar.col = c("#B2182B"),
xlim = c(0.1,0.7), ylim = c(0.1,0.7), col = c("#B2182B"), add = T)
lines(cal500[,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("300-Day","500-Day"), # 图例文字
col =c("#2166AC","#B2182B"), # 图例线的颜色，与文字对应
lwd = 2, # 图例中线的粗细
cex = 1.2, # 图例字体大小
bty = "n") # 不显示图例边框
rms-3
参考：
1.用诺模图可视化你的模型(https://www.jianshu.com/p/bb847e41dd92)
网友评论