R语言Logistic回归模型深度验证以及Nomogram绘制

小易学统计 互联网医疗统计师|自立，爱家人 15 人 赞同了该文章

01 研究背景

本章将常用的基于R语言实现二元Logistic回归模型临床预测模型的构建和验证，以及诺曼图的绘制记录下来，更为复杂的生存分析中的Cox回归将在后续章节介绍。临床预测模型的思路总结如下：①明确临床问题，确定科学假设。②查找文献，确定预测模型的研究思路。③确定模型中结局变量。④确定模型中的预测因子。⑤构建模型，计算模型预测值。⑥模型区分度评估。⑦模型校准度评估。⑧临床实用型DCA评估。

02 案例研究

本文采用的数据是上海交大出版<医学统计学及SAS应用>第十一章数据。预测因子有性别、年龄和高血压等级，结局变量为是否患病。本文研究目的探讨患病的危险因素构建并验证模型。因数据量少且只有一个数据集，故只用此数据集建模，并验证，若有更多外部数据，最好拿外部数据来验证模型。

临床研究一般有提供多个危险因素，首先做单因素的筛选，具体筛选方法，见公众号之前的文章。筛选完的危险因素用来构建预测模型。

具体分析步骤是，①基于这些变量构建模型。②绘制Nomogram图。③计算模型ROC曲线面积(区分度)和绘制校准曲线并检验(校准度,U检验)，该步骤用神包rms一步实现。接下来直接上代码。

03 R代码及解读

library(rms)   ###加载rms包#
##建立数据集
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)  ##把数据集设置成数据框结构

str(dt)  ##查看每个变量结构
'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) ##查看数据框前几行
  y age sex ECG
1 0  28   0   0
2 1  42   1   0
3 0  46   0   1
4 0  45   1   1

第一步，构建模型。

#设定环境参数#
ddist <- datadist(dt)
options(datadist='ddist')
###logistic 
f <- lrm(dt$y~.,data=dt)   ##注意此处使用lrm()函数构建二元LR
summary(f)   ##也能用此函数看具体模型情况，模型的系数，置信区间等

第二步：绘制nomogram图，注意该函数里面的参数设置。

###  nomogram
par(mgp=c(1.6,0.6,0),mar=c(2,2,2,2))  ##设置画布
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,  ##是否显示系数轴
                     conf.int = F, ##每个得分的置信度区间，用横线表示,横线越长置信度越
                     abbrev = F#是否用简称代表因子变量
                    )
plot(nomogram)

该图的使用，本质上这是将逻辑回归模型可视化展示，方便临床快速判断。假设有个病人年龄45岁，性别为男，高血压正常，Nomogram用法是在age变量上找到其值为45的刻度，然后画垂线投影到最上方的points刻度尺上，找到找到对应的分值为50分，同理找到sex为1的分值约为37分，ECG为0对应分值为0，将这三个因素的points值加起来总分87。下一步在下面的Total Points刻度尺上找到87，向下方的Risk of Death做垂线，87对应的值在0.4和0.5之间，约为0.48，说明该患者患病风险预测概率值为48%。

重头戏来了，第三步：利用rms包对该模型进行验证。

##模型验证
##以原数据集为验证集
f.glm <- glm(y~.,data=dt,family = binomial(link = "logit"))
P1 <- predict(f.glm,type = 'response')  ##获得预测概率值
##关键的一步来了。
val.prob(P1,y)     ##这个函数前面放概率值，后面芳结局变量
          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 

该函数可以一次性得到模型验证的多个指标和P值，并绘制出校准曲线，功能很强大了。

首先看代码中返回的结果，Emax是模型与理想模型的最大偏移量，Eavg是模型与理想模型的最小偏移量，这两个值越小越好，越小则说明模型与理想模型越接近。U是指Unreliability test 即U检验，用来判断构建的模型是否能通过校准度检验，其对应的P值在最下面，S:p，当S:p>0.05说明通过校准度检验。C(ROC)是ROC面积，该面积和C-index指数本质上是一样的，只不过一个对应LR,一个对应COX。

通过R计算的结果可看到，本模型通过校准度检验，p=0.998＞0.05，Roc面积为0.784具有良好的区分能力，总体来说，该模型的预测能力是很优秀的~~。

04 总结

本文介绍了Logistic回归模型的深度验证和Nomogram的绘制及应用。需要注意的是：一个预测模型的好坏除了内部验证，还要看外部验证，即它的外推性是否好。本文由于数据量少，也没有获取外部验证集，仅用原始数据集作为训练集和验证集。

05 更多阅读

文章在公粽号：易学统计

文章里的干货更多哟

欢迎交流，欢迎提问