在上一篇广义线性模型一（Generalized Linear Models，GLM） - 简书 (jianshu.com)，我们大致了解了glm的应用范围，并从三个方面探索模型构建：

1、如何使用glm构建logistics回归
2、如何提取模型中的参数
3、不同模型之间如何比较

接下来，我们继续从四个方面谈一谈logistics回归：

1、模型可视化(列线图)
2、使用构建的模型预测结局
3、评价模型的效能
区分度(discrimination ability)：C指数、ROC曲线、NRI、IDI
一致性（calibration ability）：校准曲线
4、临床效能分析（Decision Curve Analysis）：DCA, 生存曲线

【参考】rms.pdf (r-project.org)

以上四个方面并不能完全分开来讲，因此我将可以使用同一R代码的部分一起讲解，如下：列线图的构建和校准曲线、DCA放在一起，预测结局变量及ROC放在一起。

一、列线图及校准曲线

在上一篇文章中，构建logistic回归我们使用了glm()函数，构建方式如下:
fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, data=Affairs, family=binomial())
当我们要对其进行可视化时，就像我在cox回归模型的展示中说的，想要使用rms绘制列线图，就只能使用rms来构建模型

survival包和rms包都生成原材料比如都养猪，同时rms还做香肠，rms还只用自己家的猪。害，所以没办法，你家猪再优质，想吃香肠就只能去找rms包
survival包学习笔记——cox回归（一） - 简书 (jianshu.com)

列线图

继续上一篇构建的模型，我们来绘制列线图

library(rms)
ddist <- datadist(Affairs)
options(datadist='ddist')
fit.reduced1 <- lrm(ynaffair ~ age + yearsmarried + religiousness + 
                      rating, data=Affairs,x = T,y = T)
summary(fit.reduced1 )
nomo <- nomogram(fit.reduced1,
                 lp = F,                                #是否为线性
                 fun = plogis,                          #系数转换的具体过程
                 fun.at = c(0.1,0.2,0.3,0.5,0.7,0.9),   #最后一行的分段
                 funlabel = "离婚率" )                   #最后一行的名字
plot(nomo)

image.png

这里我们每次使用只需要更改数据集的名字即可啦！

校准曲线

一致性分析：校准曲线

首先绘制校准曲线，校准曲线所依赖的模型是上步rms::lrm()构建的模型，因此绘制完列线图可以直接绘制校准曲线。

library(rms)
ddist <- datadist(Affairs)
options(datadist='ddist')
fit.reduced1 <- lrm(ynaffair ~ age + yearsmarried + religiousness + 
                      rating, data=Affairs,x = T,y = T)
summary(fit.reduced1 )


calb <- calibrate(fit.reduced1, 
                  cmethod='hare',    #logistic回归使用的方法
                  method='boot',    #使用抽样的办法，抽样1000次
                  B=1000)
plot(calb,xlim=c(0,1.0),ylim=c(0,1.0))

image.png

二、模型预测及ROC曲线绘制

模型已经拟合好了，那么他是否可以根据我们的变量来预测结局呢？

结局预测

那是当然，使用函数predict

S3 method for class 'glm'

predict(object, newdata = NULL, type = c("link", "response", "terms"), se.fit = FALSE, dispersion = NULL, terms = NULL, na.action = na.pass, ...)

S3 method for class 'lm'

predict(object, newdata, se.fit = FALSE, scale = NULL, df = Inf, interval = c("none", "confidence", "prediction"), level = 0.95, type = c("response", "terms"), terms = NULL, na.action = na.pass, pred.var = res.var/weights, weights = 1, ...)

参数介绍

fit: 这里fit就是刚刚拟合的函数，
newdata: 可以是多种多样的
1、可以是原来的数据，预测后可以用来绘制ROC曲线
2、可以是新的数据，来计算敏感性、特异性
3、可以是自己生成的数据，来展示某一变量的影响情况
type: 预测数据的类型
1、response，fit是线性模型时，返回值是规范化后的连续值; fit是二分类模型时，返回值是发生结局变量的可能性（probabilities）
2、term: 返回一个矩阵，里面包含这个模型中的每个参数项的预测值，每一项之和或者其他运算得到返回值（formular中的y值）

实战

这三个参数，已经够我们使用的了，接下来我们实战看一下

library(rms)
library(xlsx)
library(tidyverse)

fit <- glm(HPD与否~.,data = data,family = binomial())
summary(fit)

image

a <- predict(fit,newdata = data,type = "terms",se.fit = T)
b <- predict(fit,newdata = data,type = "response")

type="terms"

term="response"

计算C指数及ROC曲线绘制

区分度——C指数

rms包计算C指数

library(rms)
ddist <- datadist(Affairs)  #将数据组装
options(datadist='ddist')
fit.reduced1 <- lrm(ynaffair ~ age + yearsmarried + religiousness + 
                      rating, data=Affairs,x = T,y = T)
 #相对于基础包中的glm，这里多了 x 和 y 参数
fit.reduced1 

image.png

ROCR包绘制ROC曲线

这里使用ROCR必须依赖rms得到的fit，同列线图

library(rms)
ddist <- datadist(Affairs)  #将数据组装
options(datadist='ddist')
fit.reduced1 <- lrm(ynaffair ~ age + yearsmarried + religiousness + 
                      rating, data=Affairs,x = T,y = T)


Affairs$predvalue<-predict(fit.reduced1)
#这里已经在前文中包装了相关的数据，所以无需赋值newdata，type
library(ROCR)
pred <- prediction(Affairs$predvalue, Affairs$ynaffair) 
 #结局变量必须是二分类变量（0，1），否则结果会报错
perf<- performance(pred,"tpr","fpr")
plot(perf)
abline(0,1)
auc <- performance(pred,"auc")
auc@y.values #auc即是C-statistics

perf

AUC

ROC

ROCR得到perf后，使用plot(perf)绘制的ROC曲线较简单，接下来我们使用pROC优化ROC曲线

pROC包

输入变量只需要predict得到的预测值和准确值

library(pROC)
pdf(file = 'roc.pdf',width = 5,height = 5)
plot.roc(Affairs$ynaffair~Affairs$predvalue,    #输入变量，实际变量~预测变量
         data = Affairs,
         axes=T, ## 是否显示xy轴
         legacy.axes=T, ## TRUE为1-specificity,FLASE为specificity
         main="ROC for HPD与否", ## Title
         
         col= "black", ## 曲线颜色
         lty=1, ## 曲线形状
         lwd=3, ## 曲线粗细
         identity=T, ## 是否显示对角线
         identity.col="black", ## 对角线颜色
         identity.lty=2, ## 对角线形状
         identity.lwd=2, ## 对角线粗细
         
         print.thres=TRUE, ## 是否输出cut-off值
         print.thres.best.method="youden",
         print.thres.pch=20, ## cut-off点的形状
         print.thres.col="black", ## cut-off点和文本的颜色
         print.thres.cex=1.2, ## 文本放大的倍数（比起原始值） 
         #print.thres.pattern="%.3f", ## cut-off文本的格式
         #print.thres.pattern.cex=1.2,
         
         print.auc=T, ## 是否显示AUC
         ci=T, ## 是否显示CI
         print.auc.pattern="AUC = 0.704", ## 展示AUC的格式
         print.auc.x=0.4, ## AUC值的X位置
         print.auc.y=0.15, ## AUC值的Y位置
         print.auc.cex=1.2, ## AUC值的放大倍数
         print.auc.col='black', ## ACU值的颜色
         auc.polygon=TRUE, ## 是否将ROC下面积上色
         auc.polygon.col='white', 
         auc.polygon.density=NULL,
         auc.polygon.lty=2,
         auc.polygon.angle=45,
         auc.polygon.border='black',
         
         max.auc.polygon=TRUE,
         max.auc.polygon.col='white',
         max.auc.polygon.lty=2
)
dev.off()

ROC曲线

最后

由于篇幅的问题，剩余的内容继续在下一篇中分享