在上一篇广义线性模型一（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, 生存曲线
以上四个方面并不能完全分开来讲，因此我将可以使用同一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曲线
最后
由于篇幅的问题，剩余的内容继续在下一篇中分享
网友评论