rm(list = ls())
lowbirth <- read.csv("./datasets/lowbirth.csv")31 logistic回归校准曲线绘制
评价模型的好坏主要看区分度和校准度，校准度方面目前最推荐的还是校准曲线（calibration curve），可用于评价模型预测概率和实际概率一致性。
除了我介绍的这几种方法外，还有predtools、caret等可以用于分类变量的校准曲线绘制。如果你明白了校准曲线就是真实概率和预测概率的分箱平均值散点图，你其实可以自己画，并不局限于logistic和cox，像随机森林、lasso、SVM等很多模型都可以画出校准曲线。
拟合优度检验（Hosmer-Lemeshow goodness-of-fit test）可以用来比较预测概率和实际概率是否有显著性差异，但是这个检验也只是能说明两者有没有统计学意义，并不能说明好多少、差多少。
31.1 加载数据
使用lowbirth数据集，这个数据集是关于低出生体重儿是否会死亡的数据集，其中dead这一列是结果变量，0代表死亡，1代表存活，其余列都是预测变量。
注意：需要把分类变量因子化，对于无序分类变量，需要设置哑变量！
查看一下数据：
dim(lowbirth) # 565行，10列
## [1] 565 10
str(lowbirth)
## 'data.frame': 565 obs. of 10 variables:
## $ birth : num 81.5 81.6 81.6 81.6 81.6 ...
## $ lowph : num 7.25 7.06 7.25 6.97 7.32 ...
## $ pltct : int 244 114 182 54 282 153 229 182 361 378 ...
## $ race : chr "white" "black" "black" "black" ...
## $ bwt : int 1370 620 1480 925 1255 1350 1310 1110 1180 970 ...
## $ delivery: chr "abdominal" "vaginal" "vaginal" "abdominal" ...
## $ apg1 : int 7 1 8 5 9 4 6 6 6 2 ...
## $ vent : int 0 1 0 1 0 0 1 0 0 1 ...
## $ sex : chr "female" "female" "male" "female" ...
## $ dead : int 0 1 0 1 0 0 0 0 0 1 ...简单的把人种分为白色和黑色人种（无序分类变量需要设置哑变量），再去掉race这一列，然后其余分类变量因子化。
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
tmp <- lowbirth %>%
mutate(across(where(is.character),as.factor),
vent = factor(vent),
black = ifelse(race == "black",1,0),
white = ifelse(race == "white",1,0),
other = ifelse(race %in% c("native American","oriental"),1,0)
) %>%
select(- race)
glimpse(tmp)
## Rows: 565
## Columns: 12
## $ birth <dbl> 81.514, 81.552, 81.558, 81.593, 81.610, 81.624, 81.626, 81.68…
## $ lowph <dbl> 7.250000, 7.059998, 7.250000, 6.969997, 7.320000, 7.160000, 7…
## $ pltct <int> 244, 114, 182, 54, 282, 153, 229, 182, 361, 378, 255, 186, 26…
## $ bwt <int> 1370, 620, 1480, 925, 1255, 1350, 1310, 1110, 1180, 970, 770,…
## $ delivery <fct> abdominal, vaginal, vaginal, abdominal, vaginal, abdominal, v…
## $ apg1 <int> 7, 1, 8, 5, 9, 4, 6, 6, 6, 2, 4, 8, 1, 8, 5, 9, 9, 9, 6, 2, 1…
## $ vent <fct> 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1…
## $ sex <fct> female, female, male, female, female, female, male, male, mal…
## $ dead <int> 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0…
## $ black <dbl> 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0…
## $ white <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1…
## $ other <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…31.2 方法1：rms
使用rms包。
library(rms)
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "ndiMatrix" of class "replValueSp"; definition not updated
dd <- datadist(tmp)
options(datadist="dd")
fit2 <- lrm(dead ~ birth + lowph + pltct + bwt + vent + black + white,
data = tmp,x=T,y=T)
cal2 <- calibrate(fit2, method='boot', B=500)接下来就是画图：
plot(cal2,
xlim = c(0,1),
ylim = c(0,1),
xlab = "Prediced Probability",
ylab = "Observed Probability",
cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=0.8,
#subtitles = FALSE,
legend = FALSE
)
##
## n=565 Mean absolute error=0.012 Mean squared error=0.00029
## 0.9 Quantile of absolute error=0.032
lines(cal2[,c("predy","calibrated.corrected")],
type = 'l', #连线的类型，可以是"p","b","o"
lwd = 3, #连线的粗细
pch = 16, #点的形状，可以是0-20
col = "#2166AC") #连线的颜色
lines(cal2[,c("predy","calibrated.orig")],type="l",pch=16,lwd=3,col="tomato")
abline(0,1,
lty = 2, #对角线为虚线
lwd = 2, #对角线的粗细
col = "#224444") #对角线的颜色
legend(0.6,0.2,
c("Apparent","Bias-corrected","Ideal"),
lty = c(2,1,1),
lwd = c(2,3,3),
col = c("black","#2166AC","tomato"),
bty = "n"
)
还可以计算hosmer-lemeshow的P值。
# 进行hosmer-lemeshow 检验
library(ResourceSelection)
## ResourceSelection 0.3-6 2023-06-27
model_glm <- glm(dead ~ birth + lowph + pltct + bwt + vent + black + white,
data = tmp, family = binomial)
# hosmer-lemeshow 检验
p.hoslem <- hoslem.test(model_glm$y, fitted(model_glm), g=10)$p.value
p.hoslem
## [1] 0.2340365
# 构建 calibration
fit2 <- lrm(dead ~ birth + lowph + pltct + bwt + vent + black + white,
data = tmp,x=T,y=T)
cal2 <- calibrate(fit2, method='boot', B=500)画图还是和上面一样，就是多了一个添加 hosmer-lemeshow P值的步骤。
plot(cal2,
xlim = c(0,1),
ylim = c(0,1),
xlab = "Prediced Probability",
ylab = "Observed Probability",
cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=0.8,
#subtitles = FALSE,
legend = FALSE
)
##
## n=565 Mean absolute error=0.012 Mean squared error=0.00028
## 0.9 Quantile of absolute error=0.032
lines(cal2[,c("predy","calibrated.corrected")],
type = 'l', #连线的类型，可以是"p","b","o"
lwd = 3, #连线的粗细
pch = 16, #点的形状，可以是0-20
col = "#2166AC") #连线的颜色
lines(cal2[,c("predy","calibrated.orig")],type="l",pch=16,lwd=3,col="tomato")
abline(0,1,
lty = 2, #对角线为虚线
lwd = 2, #对角线的粗细
col = "#224444")#对角线的颜色
legend(0.6,0.2,
c("Apparent","Bias-corrected","Ideal"),
lty = c(2,1,1),
lwd = c(2,3,3),
col = c("black","#2166AC","tomato"),
bty = "n"
)
text(0,0,bquote("Hosmer-Lemeshow "~italic(P)~" = "~.(round(p.hoslem,3))),adj = 0)
31.3 方法2：riskRegression
使用riskRegression包。
fit2 <- glm(dead ~ birth + lowph + pltct + bwt + vent + black + white,
data = tmp, family = binomial)library(riskRegression)
## riskRegression version 2023.09.08
fit22 <- Score(list("fit"=fit2),
formula = dead ~ 1,
data = tmp,
metrics = c("auc","brier"),
summary = c("risks","IPA","riskQuantile","ibs"),
plots = "calibration",
null.model = T,
conf.int = T,
B = 500,
M = 50
)
plotCalibration(fit22,col="tomato",
xlab = "Predicted Risk",
ylab = "Observerd RISK",
bars = F)
非常神奇的是，还可以用ggplot2来画！
plotdata <- plotCalibration(fit22,plot = F,method = "nne"
#bandwidth = 0.1
)
library(ggplot2)
ggplot(plotdata$plotFrames$fit, aes(x=Pred,y=Obs))+
geom_line(color="tomato",linewidth=1.5)+
scale_x_continuous(limits = c(0,1),name = "Predicted Risk")+
scale_y_continuous(limits = c(0,1),name = "Observerd Risk")+
geom_abline(slope = 1,intercept = 0,lty=2)+
geom_rug(color="grey")+
theme_bw()
31.4 多个calibration画在一起
31.4.1 方法1
首先构建多个画图需要的结果：
fit2 <- lrm(dead ~ birth + lowph + pltct + bwt + vent + black + white,
data = tmp,x=T,y=T)
cal2 <- calibrate(fit2, method='boot', B=500)
fit1 <- lrm(dead ~ lowph + bwt + vent + black,
data = tmp,x=T,y=T)
cal1 <- calibrate(fit1, method='boot', B=500)
fit3 <- lrm(dead ~ birth + lowph + pltct + bwt + vent,
data = tmp,x=T,y=T)
cal3 <- calibrate(fit3, method='boot', B=500)然后使用基础函数画在一起即可：
plot(1,
type = "n",
xlim = c(0,1),
ylim = c(0,1),
xlab = "Nomogram-prediced Probability",
ylab = "Observed Probability",
cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=0.8
)
lines(cal1[,c("predy","calibrated.corrected")],
type = 'l', #连线的类型，可以是"p","b","o"
lwd = 3, #连线的粗细
pch = 16, #点的形状，可以是0-20
col = "#2166AC") #连线的颜色
lines(cal2[,c("predy","calibrated.corrected")],type="l",pch=16,lwd=3,col="tomato")
lines(cal3[,c("predy","calibrated.corrected")],type="l",pch=16,lwd=3,col="skyblue")
abline(0,1,
lty = 2, #对角线为虚线
lwd = 2, #对角线的粗细
col = "#224444")#对角线的颜色
legend(0.6,0.2,
c("fit1","fit2","fit3"),
lty = c(1,1,1),
lwd = c(3,3,3),
col = c("#2166AC","tomato","skyblue"),
bty = "n"
)
31.4.2 方法2
library(riskRegression)
fit1 <- glm(dead ~ lowph + bwt + vent + black,data = tmp,family=binomial)
fit2 <- glm(dead ~ birth + lowph + pltct + bwt + vent + black + white,
data = tmp, family = binomial)
fit3 <- glm(dead ~ birth + lowph + pltct + bwt + vent,data=tmp,
family=binomial)
fit_all <- Score(list("Fit1"=fit1,
"Fit2"=fit2,
"Fit3"=fit3
),
formula = dead ~ 1,
data = tmp,
metrics = c("auc","brier"),
summary = c("risks","IPA","riskQuantile","ibs"),
plots = "calibration",
null.model = T,
conf.int = T,
B = 500,
M = 50
)
data_all <- plotCalibration(fit_all,plot = F)接下来数据稍作变换，就可以画图了！
plot_df <- bind_rows(data_all$plotFrames) %>%
mutate(fits = rep(c("fit1","fit2","fit3"),each=76))
ggplot(plot_df, aes(Pred,Obs))+
geom_line(aes(group=fits,color=fits),linewidth=1.2)+
scale_color_manual(values = c("#2166AC","tomato","skyblue"),name=NULL)+
scale_x_continuous(limits = c(0,1),name = "Predicted Risk")+
scale_y_continuous(limits = c(0,1),name = "Observerd Risk")+
geom_abline(slope = 1,intercept = 0,lty=2)+
geom_rug(aes(color=fits))+
theme_bw()
获取lowbirth数据请在公众号后台回复20220520。