原文地址:https://blog.csdn.net/anshiquanshu/article/details/53444800
首先什么是Nomogram鲤屡?簡單的說這是一種將Logistic回歸或Cox回歸圖形化呈現(xiàn)的方法福侈,可以讓讀者從圖中很簡便地根據(jù)預測變量的值得到因變量的大致概率數(shù)值肪凛。其對于Logistic回歸或Cox回歸的意義辽社,大概相當于散點圖對于簡單線性回歸的意義衡奥。具體的介紹以及作圖原理,這里就不詳述了,有興趣的請參照附件中SAS公司的一份文檔档址。
下面簡單說下Nomogram怎么看邻梆。如下圖。欲知年齡50歲的女性(sex=1)的患病風險浦妄,只需要將age=45歲向points軸投射剂娄,則points=50;同理sex=1時和二,points≈37耳胎。兩者相加則Total
points=87;將此數(shù)值在Total points軸上向Risk概率軸投射废登,則可知風險大概在0.4和0.5之間郁惜。(參見圖中紅線)對于單個變量,只需要令Total points = points進行投射即可吏颖。
接下來講如何用R語言做出上面的這張圖恨樟。簡單起見劝术,此帖僅討論Logistic回歸呆奕,Cox回歸的方法類似衬吆,但相對更復雜。本帖所用數(shù)據(jù)引用自上海交大出版《醫(yī)學統(tǒng)計學及SAS應用(修訂版)》的例11.4姆泻,詳細的內容請參見我的另一個帖子的附件(http://www.dxy.cn/bbs/topic/26880076)冒嫡。
-------------------------------------代碼開始---------------------------------------------
require(rms)
# 建立數(shù)據(jù)集
y = c(0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0,
1, 0, 0, 0, 1, 1, 0, 1,
1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1,
0, 1, 0, 0, 0, 1,
1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1,
0, 0, 1, 1, 1, 0,
0, 0, 1, 0, 1, 0, 1, 0, 1)
age = c(28, 42, 46, 45, 34, 44, 48, 45, 38, 45, 49,
45, 41, 46, 49, 46, 44, 48,
52, 48, 45, 50, 53, 57, 46, 52, 54, 57, 47, 52, 55,
59, 50, 54, 57, 60,
51, 55, 46, 63, 51, 59, 48, 35, 53, 59, 57, 37, 55,
32, 60, 43, 59, 37,
30, 47, 60, 38, 34, 48, 32, 38, 36, 49, 33, 42, 38,
58, 35, 43, 39, 59,
39, 43, 42, 60, 40, 44)
sex = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0,
1, 0, 1, 0, 1, 0, 1, 0, 1,
0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0,
1, 0, 1, 0, 1, 0, 1,
0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0,
1, 1, 1, 0, 1, 1, 1,
0, 1, 1, 1, 0, 1)
ECG = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1,
1, 0, 0, 1, 1, 0, 0, 1, 1,
0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 1, 0,
0, 2, 2, 0, 0, 2, 2,
0, 1, 2, 2, 0, 1, 0, 2, 0, 1, 0, 2, 1, 1, 0, 2, 1,
1, 0, 2, 1, 1, 0, 2,
1, 1, 0, 2, 1, 1)
# 設定nomogram的參數(shù)
ddist <- datadist(age, sex, ECG)
options(datadist='ddist')
# logistic回歸
f <- lrm(y ~ age + sex + ECG)
# nomogram
nom <- nomogram(f, fun=plogis,
fun.at=c(.001, .01, .05, seq(.1,.9, by=.1), .95, .99,
.999),
lp=F, funlabel="Risk")
plot(nom)
-------------------------------------代碼結束---------------------------------------------
END方咆。歡迎交流蟀架。
Cox回歸模型會復雜一些,因為可能涉及到不同時間點(3年煌集、5年)生存概率的計算穆碎。下面討論最簡單的概率軸為中位生存時間的情況。若需要在映射軸呈現(xiàn)生存概率方面,或需要比較不同模型之間的優(yōu)劣(比如c-index)色徘,另外再留言討論。
require(rms)
require(Hmisc)
# 建立數(shù)據(jù)集(使用rms包example的代碼横腿,未改動)
n <- 1000
set.seed(731)
age <- 50 + 12*rnorm(n)
label(age) <- "Age"
sex <- factor(sample(c('Male','Female'), n,
rep=TRUE, prob=c(.6, .4)))
cens <- 15*runif(n)
h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
dt <- -log(runif(n))/h
label(dt) <- 'Follow-up Time'
e <- ifelse(dt <= cens,1,0)
dt <- pmin(dt, cens)
units(dt) <- "Year"
# 設定nomogram的參數(shù)
ddist <- datadist(age, sex)
options(datadist='ddist')
# Cox回歸
S <- Surv(dt,e)
f <- cph(S ~ rcs(age,4)
+ sex, x=T, y=T)
med <- Quantile(f)
# nomogram
nom <- nomogram(f, fun=function(x) med(x),
fun.at=c(13,12,11,9,8,7,6,5),
lp=F, funlabel="Median Survival Time")
plot(nom)