- 模拟数据
n <- 1000 # sample size
set.seed(88)
age <- rnorm(n,65,11)
lac <- round(abs(rnorm(n,3,1)),1)
sex <- factor(sample(1:2,n,prob = c(0.6,0.4),TRUE),
labels = c("male","female"))
shock <- factor(sample(1:4,n,prob = c(0.3,0.3,0.25,0.15),TRUE),
labels = c("no","mild","moderate","severe"))
z <- 0.2*age +3*lac*as.numeric(sex)+5*as.numeric(shock)-rnorm(n,36,15) # linear combination with a bias
y <- ifelse(runif(n)<=plogis(z),1,0)
# plogis()为logit的反函数，根据线性函数计算出概率probability
Y <- ifelse(y == 0,0,sample(1:3,length(y),TRUE))
data <- data.frame(age = age, lac = lac,sex=sex,shock=shock,y =y,Y =Y)
var.labels = c(age = "Age in Years",
lac = "lactate",
sex = "Sex of the participant",
shock = "shock",
y = "outcome",
Y = "ordinal")
library(rms)
label(data) = lapply(names(var.labels),
function(x) label(data[,x])=var.labels[x]) #label()为rms里的函数
head(data)
image.png
- 绘制nomogram
library(rms)
ddist <- datadist(data)
options(datadist="ddist")
mod.bi <- lrm(y~shock+lac*sex+age,data)
nom.bi <- nomogram(mod.bi,
lp.at = seq(-3,4,by = 0.5),#轴上标注点
fun = function(x)1/(1+exp(-x)), #线性轴转换成概率，显示在另一条轴上
fun.at = c(.001,.01,.05,seq(.1,.9,by=.1),.95,.99,.999),
funlabel = "Risk of Death",
conf.int = c(0.1,0.7),
abbrev = TRUE, # 对于交互作用的factor各个水平进行缩写
minlength=1,lp =F)
plot(nom.bi,lplabel ="Linear Predictor",
fun.side = c(3,3,1,1,3,1,3,1,1,1,1,1,3), #标度在上或在下
col.conf = c("red","green"),
conf.space = c(0.1,0.5),
label.every = 3, # 每3个tick标注一次
col.grid = gray(c(0.8,0.95)),
which = "shock")
legend.nomabbrev(nom.bi, which = "shock",x =.5,y = .5) #简写
image.png
- 结局变量具有多个水平(ordinal logistic regression model)
mod.ord <- lrm(Y~age+rcs(lac,4)*sex)
fun2 <- function(x) plogis(x-mod.ord$coef[1]+mod.ord$coef[2])
fun3 <- function(x) plogis(x-mod.ord$coef[1]+mod.ord$coef[3])
image.png
f <- Newlabels(mod.ord,c(age="Age(years)")) #对label重新定义
nom.ord <- nomogram(f,fun = list("Prob Y >= 1" = plogis,
"Prob Y >= 2" = fun2,
"Prob Y >= 3" = fun3),
lp = F,
fun.at = c(.01,.05,seq(.1,.9,by = .1),.95,.99))
plot(nom.ord,lmp =.2,cex.axis = .6)
image.png
- 生存资料
library(survival)
str(lung)
image.png
参数模型
lung$sex <- factor(lung$sex,labels = c("male","female"))
mod.sur <- psm(Surv(time,status)~ph.ecog+sex+age,
lung,dist = "weibull")
med <- Quantile(mod.sur) #将线性预测转化为中位生存时间
surv <- Survival(mod.sur) #将线性预测转化为某个时间点的生存概率
ddist <- datadist(lung)
options(datadist = "ddist")
nom.sur1 <- nomogram(mod.sur,
fun = list(function(x) med(lp=x,q=0.5),
function(x) med(lp=x,q=0.25)),
funlabel = c("Median Survival Time",
"1Q Survival Time"),
lp = F)
plot(nom.sur1,
fun.side = list(c(rep(1,7),3,1,3,1,3),rep(1,7)),
col.grid = c("red","green"))
image.png
在某个时间点的生存概率
nom.sur2 <- nomogram(mod.sur,fun = list(function(x)
surv(200,x),
function(x) surv(400,x)),
funlabel = c("200-Day Survival Probability",
"400-Day Survival Probability"),
lp = F)
plot(nom.sur2,
fun.side = list(c(rep(c(1,3),5),1,1,1,1),
c(1,1,1,rep(c(3,1),6))),
xfrac=.7,
col.grid = c("red","green"))
image.png
半参数模型
ddist <- datadist(lung)
options(datadist="ddist")
mod.cox <- cph(Surv(time,status)~ph.ecog+sex+age,
lung,surv=TRUE)
surv.cox <- Survival(mod.cox)
nom.cox <- nomogram(mod.cox,fun = list(function(x) surv.cox(200,x),
function(x) surv.cox(400,x)),
funlabel = c("200-Day Sur.Prob.","400-Day Sur.Prob."),
lp = F)
plot(nom.cox,
fun.side = list(c(rep(c(1,3),5),1,1,1,1),c(1,1,1,rep(c(3,1),6))))
image.png
参考资料
文中代码及部分截图来自章仲恒教授的丁香园课程：回归模型可视化：列线图制作
