- 模擬數(shù)據(jù)
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的反函數(shù)姑宽，根據(jù)線性函數(shù)計算出概率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里的函數(shù)
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)), #線性軸轉(zhuǎn)換成概率腋颠，顯示在另一條軸上
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
- 結(jié)局變量具有多個水平(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
參數(shù)模型
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) #將線性預(yù)測轉(zhuǎn)化為中位生存時間
surv <- Survival(mod.sur) #將線性預(yù)測轉(zhuǎn)化為某個時間點的生存概率
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
半?yún)?shù)模型
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
參考資料
文中代碼及部分截圖來自章仲恒教授的丁香園課程：回歸模型可視化：列線圖制作
