临床模型篇第三弹，如何用R语言对Logistic回归模型进行评价？

大家好，我是阿琛。在常用的临床模型构建中，主要分为两种，包括临床预测模型（Cox回归模型）和临床诊断模型（Logstic回归模型）。在之前的内容中，阿琛给大家介绍了如何使用Nomogram图将临床预测模型可视化，以及Cox回归模型的相关评价指标。

传送门

1. Nomogram图不会画？看了这篇，小白也能轻松看懂搞定；

2. 临床预测模型如何评价？看这篇就够了。

今天是临床模型的第三集临床诊断模型篇，即如何对Logistic回归模型进行评价。

1.建立Logistic预测模型

1.1 引用R包

#install.packages("foreign")#install.packages("rms")#install.packages("pROC")#install.packages("rmda")#install.packages("nricens")library(foreign)library(rms)   #构建Logstic模型library(pROC)   #绘制ROC曲线library(rmda)   #绘制DCA曲线library(nricens)   #计算NRI值

1.2 读取文件

setwd("C:\Users\000\Desktop\14_Logstic")   #设置工作目录rt <- read.table("clinical_2.txt",header=T,sep="  ")   #读取数据head(rt)   #查看前5行的数据

可以发现，该数据集主要包括5个不同的变量和1个结局变量；在此，我们主要通过构建临床诊断模型，来预测结局Status。

1.3 构建Logstic模型

首先，对数据集中的相关参数进行设置；

rt$Age <- factor(rt$Age,labels=c("<60",">=60"))rt$Gender <-  factor(rt$Gender,labels=c("No","Yes"))rt$Education <- factor(rt$Education,labels=c("Primary","Secondary","Higher"))rt$Alcohol <- factor(rt$Alcohol,labels=c("No","Yes"))str(rt)ddist <- datadist(rt)options(datadist = "ddist")   #使用函数datadist()将数据打包

随后，使用glm()函数构建Logstic模型，在R中glm()函数可以拟合广义线性模型，其中包括Logstic回归模型。其中，我们必须使用family=binomial这个参数来告诉R运行Logstic回归模型，而不是其他的广义线性模型；

fit <- glm(Status~Age + Gender + BMI + Education + Alcohol, data=rt,           family = binomial(link="logit"), x=T) #利用lrm(）函数对模型进行拟合 summary(fit)   #查看模型拟合结果

通过summary()函数，我们可以查看各个预测变量的系数及其P值；可以看到，只有两个特征的P值小于0.05（BMI和Education）；

2.计算C指数及绘制ROC曲线

2.1 C-Index的计算

Cindex <- rcorrcens(Status~predict(fit), data = rt)Cindex

通过这种方式，我们也可以获得该模型的C-Index为0.718，标准差为0.061；一般而言，我们将<0.5看成模型没有任何预测能力，0.51-0.7认为是较差的准确性，0.71-0.9为中等的准确性，> 0.9为高度的准确性。

2.2 ROC曲线的绘制

接下来，加个餐，我们一起来看下如何绘制预测模型的ROC曲线；

gfit <- roc(Status~predict(fit), data = rt)plot(gfit,    print.auc=TRUE,   #输出AUC值    print.thres=TRUE,   #输出cut-off值    main = "ROC CURVE",   #设置图形的标题    col= "red",   #曲线颜色    print.thres.col="black",   #cut-off值字体的颜色    identity.col="blue",   #对角线颜色    identity.lty=1,identity.lwd=1)

结果显示

该模型ROC曲线的曲线下面积（AUC值）为0.718；对于AUC值，其含义与C-Index较为类似，我们往往将AUC位于0.5~0.7时定义为模型的效果较低，位于0.7~0.85之间为效果一般，而位于0.85~0.95时效果很好。同时，该模型的最佳截断值（即cut-off值）为0.011，也就是说当以0.011进行分组时，两组之间具有最佳的区分度。

3.基于ggstatspot包的绘制

cal <- calibrate(fit, method="boot", B=1000)plot(cal,    xlab="Nomogram-predicted probability of nonadherence",    ylab="Actual diagnosed nonadherence (proportion)",    sub=F)

如图所示

X轴为模型预测得到的结局可能性，而Y轴为实际观察得到的值，并重复计算1000次，其中Bias-corrected为校正曲线，而对角线Ideal为理想的曲线。校正曲线与理想曲线之间越相近，说明模型的预测能力越好。在该结果中，模型具有良好的校准能力。

4.绘制DCA曲线

评估完模型的准确性，同时还需要进一步评估模型的获益率。因此，我们通过DCA曲线，来展示患者的净获益率。

modul<- decision_curve(data= rt,                      Status~Age + Gender + BMI + Education + Alcohol,                       family = binomial(link ='logit'),                      thresholds= seq(0,1, by = 0.01),                      confidence.intervals = 0.95)plot_decision_curve(modul,                   curve.names="Nonadherence prediction nomogram",   #曲线名称                   xlab="Threshold probability",   #x轴名称                   cost.benefit.axis =FALSE, col= "blue",                   confidence.intervals=FALSE,                   standardize = FALSE)

如图所示

蓝色的曲线为模型预测的获益情况，灰色的曲线为所有患者都接受干预的获益率，而横线为所有病人都不接受干预的获益率。取蓝色曲线与All的交点为起点，与None的交点为结束，在此范围内所对应的患者可以获益。

5.计算NRI指数

在此，我们对比一下BMI的纳入对模型的质量是否会产生影响。首先，构建两个不同的预测模型；

##构建模型fit_A <- glm(Status~Age + Gender + Education + Alcohol, data = rt, family = binomial(link="logit"),x=TRUE)fit_B <- glm(Status~Age + Gender + Education + Alcohol + BMI, data = rt, family = binomial(link="logit"),x=TRUE)

接着，我们分别使用了两种不同的方式，进一步对两个模型进行了比较；

#计算连续的NRI，取值为5%时就定义为重分类

NRI <- nribin(mdl.std = fit_A, mdl.new = fit_B,             updown = 'diff',             cut = 0.05, niter = 500, alpha = 0.05)

结果显示

两个模型之间比较的NRI值为0.419（95%CI：0.247-0.671），可以认为预测的风险在新旧模型中发生了重新分类；

#计算分类变量计算的NRI,只有概率超过0.011就定义为重分组了#0.011是根据之前模型C的ROC分析确定的切点NRI <- nribin(mdl.std = fit_A, mdl.new = fit_B,             updown = 'category',             cut = 0.011, niter = 500, alpha = 0.05)

结果显示

两个模型之间比较的NRI值为0.005（95%CI：0-0.021），可以认为预测的风险在新旧模型中并没有发生重新分类；

好了，今天的分享就到此结束了。大家可以使用示例数据或者自己的模型数据，进行相应的学习。