本文轉(zhuǎn)自易學統(tǒng)計R語言Logistic回歸模型深度驗證以及Nomogram繪制
研究背景
本章將常用的基于R語言實現(xiàn)二元Logistic回歸模型臨床預(yù)測模型的構(gòu)建和驗證，以及諾曼圖的繪制記錄下來，更為復(fù)雜的生存分析中的Cox回歸將在后續(xù)章節(jié)介紹。臨床預(yù)測模型的思路總結(jié)如下：①明確臨床問題，確定科學假設(shè)。②查找文獻，確定預(yù)測模型的研究思路。③確定模型中結(jié)局變量。④確定模型中的預(yù)測因子。⑤構(gòu)建模型，計算模型預(yù)測值。⑥模型區(qū)分度評估。⑦模型校準度評估。⑧臨床實用型DCA評估。
案例研究
本文采用的數(shù)據(jù)是上海交大出版<醫(yī)學統(tǒng)計學及SAS應(yīng)用>第十一章數(shù)據(jù)。預(yù)測因子有性別、年齡和高血壓等級，結(jié)局變量為是否患病。本文研究目的探討患病的危險因素構(gòu)建并驗證模型。因數(shù)據(jù)量少且只有一個數(shù)據(jù)集，故只用此數(shù)據(jù)集建模，并驗證，若有更多外部數(shù)據(jù)，最好拿外部數(shù)據(jù)來驗證模型。
臨床研究一般有提供多個危險因素，首先做單因素的篩選，具體篩選方法，見公眾號之前的文章。篩選完的危險因素用來構(gòu)建預(yù)測模型。
具體分析步驟是，①基于這些變量構(gòu)建模型。②繪制Nomogram圖。③計算模型ROC曲線面積(區(qū)分度)和繪制校準曲線并檢驗(校準度,U檢驗)，該步驟用神包rms一步實現(xiàn)。接下來直接上代碼。
R代碼及解讀
library(rms) ###加載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)
dt <- data.frame(y,age,sex,ECG) ##把數(shù)據(jù)集設(shè)置成數(shù)據(jù)框結(jié)構(gòu)
str(dt) ##查看每個變量結(jié)構(gòu)
'data.frame': 78 obs. of 4 variables:
$ y : num 0 1 0 0 0 1 1 1 0 0 ...
$ age: num 28 42 46 45 34 44 48 45 38 45 ...
$ sex: num 0 1 0 1 0 1 0 1 0 1 ...
$ ECG: num 0 0 1 1 0 0 1 1 0 0 ...
head(dt) ##查看數(shù)據(jù)框前幾行
y age sex ECG
1 0 28 0 0
2 1 42 1 0
3 0 46 0 1
4 0 45 1 1
構(gòu)建模型
設(shè)定環(huán)境參數(shù)#
ddist <- datadist(dt)
options(datadist='ddist')
f <- lrm(dt$y~.,data=dt) 注意此處使用lrm()函數(shù)構(gòu)建
summary(f) 也能用此函數(shù)看具體模型情況，模型的系數(shù)，置信區(qū)間等
繪制nomogram圖，
注意該函數(shù)里面的參數(shù)設(shè)置。
## nomogram
par(mgp=c(1.6,0.6,0),mar=c(2,2,2,2)) ##設(shè)置畫布
nomogram <- nomogram(f,fun=function(x)1/(1+exp(-x)), ##邏輯回歸計算公式
fun.at = c(0.001,0.01,0.05,seq(0.1,0.9,by=0.1),0.95,0.99,0.999),#風險軸刻度
funlabel = "Risk of Death", #風險軸便簽
lp=F, ##是否顯示系數(shù)軸
conf.int = F, ##每個得分的置信度區(qū)間，用橫線表示,橫線越長置信度越
abbrev = F#是否用簡稱代表因子變量
)
plot(nomogram)
該圖的使用，本質(zhì)上這是將邏輯回歸模型可視化展示，方便臨床快速判斷。假設(shè)有個病人年齡45歲，性別為男，高血壓正常，Nomogram用法是在age變量上找到其值為45的刻度，然后畫垂線投影到最上方的points刻度尺上，找到找到對應(yīng)的分值為50分，同理找到sex為1的分值約為37分，ECG為0對應(yīng)分值為0，將這三個因素的points值加起來總分87。下一步在下面的Total Points刻度尺上找到87，向下方的Risk of Death做垂線，87對應(yīng)的值在0.4和0.5之間，約為0.48，說明該患者患病風險預(yù)測概率值為48%。
利用rms包對該模型進行驗證。
##模型驗證
##以原數(shù)據(jù)集為驗證集
f.glm <- glm(y~.,data=dt,family = binomial(link = "logit"))
P1 <- predict(f.glm,type = 'response') ##獲得預(yù)測概率值
##關(guān)鍵的一步來了。
val.prob(P1,y) ##這個函數(shù)前面放概率值，后面芳結(jié)局變量
Dxy C (ROC) R2 D
5.675676e-01 7.837838e-01 3.164825e-01 2.578779e-01
D:Chi-sq D:p U U:Chi-sq
2.111448e+01 NA -2.564103e-02 -3.552714e-13
U:p Q Brier Intercept
1.000000e+00 2.835189e-01 1.885480e-01 -4.335689e-09
Slope Emax E90 Eavg
9.999998e-01 1.157412e-01 6.085456e-02 3.492462e-02
S:z S:p
-2.762507e-03 9.977958e-01
該函數(shù)可以一次性得到模型驗證的多個指標和P值，并繪制出校準曲線，功能很強大了。
首先看代碼中返回的結(jié)果，Emax是模型與理想模型的最大偏移量，Eavg是模型與理想模型的最小偏移量，這兩個值越小越好，越小則說明模型與理想模型越接近。U是指Unreliability test 即U檢驗，用來判斷構(gòu)建的模型是否能通過校準度檢驗，其對應(yīng)的P值在最下面，S:p，當S:p>0.05說明通過校準度檢驗。C(ROC)是ROC面積，該面積和C-index指數(shù)本質(zhì)上是一樣的，只不過一個對應(yīng)LR,一個對應(yīng)COX。
通過R計算的結(jié)果可看到，本模型通過校準度檢驗，p=0.998＞0.05，Roc面積為0.784具有良好的區(qū)分能力，總體來說，該模型的預(yù)測能力是很優(yōu)秀的~~。
總結(jié)
本文介紹了Logistic回歸模型的深度驗證和Nomogram的繪制及應(yīng)用。需要注意的是：一個預(yù)測模型的好壞除了內(nèi)部驗證，還要看外部驗證，即它的外推性是否好。本文由于數(shù)據(jù)量少，也沒有獲取外部驗證集，僅用原始數(shù)據(jù)集作為訓(xùn)練集和驗證集。
以上就是本次跟大家分享的內(nèi)容，覺得有用的話點贊、轉(zhuǎn)發(fā)哦~
更多閱讀
如何進行變量篩選和特征選擇(三)？交叉驗證
如何進行變量篩選和特征選擇(二)？最優(yōu)子集回歸
如何進行高維變量篩選和特征選擇(一)？Lasso回歸
