校準(zhǔn)曲線圖表示的是預(yù)測值和實際值的差距，作為預(yù)測模型的重要部分，目前很多函數(shù)能繪制校準(zhǔn)曲線。
一般分為兩種，一種是通過Hosmer-Lemeshow檢驗，把P值分為10等分，求出每等分的預(yù)測值和實際值的差距


另外一種是calibration函數(shù)重抽樣繪制連續(xù)的校準(zhǔn)圖

我們既往文章《手動繪制R語言Logistic回歸模型的外部驗證校準(zhǔn)曲線（Calibration curve）（2）》已經(jīng)介紹了如何繪制外部驗證模型的校準(zhǔn)曲線，今天我們來介紹一下如何繪制分類的校準(zhǔn)曲線，如下面的圖

其實如果我們了解了繪圖原理就會明白，這4個曲線就是4個概率，只要求出概率，繪制這樣的圖形非常輕松，我們今天來演示一下
我們先導(dǎo)入數(shù)據(jù)，繼續(xù)使用我們的早產(chǎn)數(shù)據(jù)，

bc<-read.csv("E:/r/test/zaochan.csv",sep=',',header=TRUE) library(ggplot2) library(rms) library(tidyverse)

這是一個關(guān)于早產(chǎn)低體重兒的數(shù)據(jù)（公眾號回復(fù)：早產(chǎn)數(shù)據(jù)，可以獲得該數(shù)據(jù)），低于2500g被認(rèn)為是低體重兒。數(shù)據(jù)解釋如下：low 是否是小于2500g早產(chǎn)低體重兒，age 母親的年齡，lwt 末次月經(jīng)體重，race 種族，smoke 孕期抽煙，ptl 早產(chǎn)史（計數(shù)），ht 有高血壓病史，ui 子宮過敏，ftv 早孕時看醫(yī)生的次數(shù)，bwt 新生兒體重數(shù)值。
我們先把分類變量轉(zhuǎn)成因子

bc$race<-ifelse(bc$race=="black",1,ifelse(bc$race=="white",2,3)) bc$smoke<-ifelse(bc$smoke=="nonsmoker",0,1) bc$race<-factor(bc$race) bc$ht<-factor(bc$ht) bc$ui<-factor(bc$ui) bc$smoke<-factor(bc$smoke)

假設(shè)我們想了解吸煙人群和不吸煙人群比較，模型的預(yù)測能力有什么不同，可以把原數(shù)據(jù)分成2個模型，分別做成校準(zhǔn)曲線，然后進(jìn)行比較，
先分成吸煙組和不吸煙組兩個數(shù)據(jù)

dat0<-subset(bc,bc$smoke==0) dat00<-dat0[,-6] dat1<-subset(bc,bc$smoke==1) dat11<-dat1[,-6]

建立兩個回歸方程

fit0<-glm(low ~ age + lwt + race + ptl + ht + ui + ftv,family = binomial("logit"),data = dat00) fit1<-glm(low ~ age + lwt + race + ptl + ht + ui + ftv,family = binomial("logit"),data = dat11)

我們做校準(zhǔn)曲線還需求出方程的概率和Y值，兩個方程都要求，

pr0<- predict(fit0,type = c("response"))#得出預(yù)測概率 y0<-dat00[, "low"] pr1<- predict(fit1,type = c("response"))#得出預(yù)測概率 y1<-dat11[, "low"]

得出了數(shù)據(jù)，概率和Y值后就可以按我們上一篇的方法做出校準(zhǔn)曲線了，我這里為了節(jié)省時間直接用我自己編寫的gg2程序代替了，就是就是把原來的步驟整合在一起，gg2程序可以做出單獨和分類的校準(zhǔn)曲線，

smoke0<-gg2(dat00,pr0,y0,group = 2,leb = "nosmoke")

做分類的時候有5個參數(shù)，前面3個是數(shù)據(jù)，概率和Y值，group = 2是固定的，leb = "nosmoke"是你想給這個分類變量取的名字，生成如下數(shù)據(jù)

接下來做吸煙組的數(shù)據(jù)

smoke1<-gg2(dat11,pr1,y1,group = 2,leb = "smoke")

把兩個數(shù)據(jù)合并最后生成繪圖數(shù)據(jù)

plotdat<-rbind(smoke0,smoke1)

生成了繪圖數(shù)據(jù)后就可以繪圖了，只需把plotdat放進(jìn)去其他不用改，當(dāng)然你想自己調(diào)整也是可以的

ggplot(plotdat, aes(x=meanpred, y=meanobs, color=gro,fill=gro,shape=gro)) + geom_line() +geom_point(size=4)+annotate(geom = "segment", x = 0, y = 0, xend =1, yend = 1)+expand_limits(x = 0, y = 0)

美化一下圖形，這樣一個用于發(fā)表的圖就做好了

ggplot(plotdat, aes(x=meanpred, y=meanobs, color=gro,fill=gro,shape=gro)) + geom_line() +geom_point(size=4)+annotate(geom = "segment", x = 0, y = 0, xend =1, yend = 1)+expand_limits(x = 0, y = 0)+scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0))+xlab("predicted probability")+ylab("actual probability")+theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+theme(legend.justification=c(1,0), legend.position=c(1,0))

我們還可以做出帶可信區(qū)間的分類校準(zhǔn)曲線，我們把概率區(qū)間調(diào)小一點，10個太多了，畫圖不夠美觀，我看很多函數(shù)都是做成5個。（當(dāng)然你不調(diào)也是可以的）

smoke0<-gg2(dat00,pr0,y0,group = 2,leb = "nosmoke",g=5) smoke1<-gg2(dat11,pr1,y1,group = 2,leb = "smoke",g=5) plotdat<-rbind(smoke0,smoke1)

得出數(shù)據(jù)后就可以繼續(xù)繪圖了

ggplot(plotdat, aes(x=meanpred, y=meanobs, color=gro,fill=gro)) + geom_errorbar(aes(ymin=meanobs-1.96*se, ymax=meanobs+1.96*se,), width=.02)+geom_point(size=4)+annotate(geom = "segment", x = 0, y = 0, xend =1, yend = 1)+expand_limits(x = 0, y = 0)+scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0))+xlab("predicted probability")+ylab("actual probability")+theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+theme(legend.justification=c(1,0),legend.position=c(1,0))

也可以加入連線，不過我這個數(shù)據(jù)加入連線感覺不是很美觀

ggplot(plotdat, aes(x=meanpred, y=meanobs, color=gro,fill=gro)) + geom_errorbar(aes(ymin=meanobs-1.96*se, ymax=meanobs+1.96*se,), width=.02)+geom_point(size=4)+annotate(geom = "segment", x = 0, y = 0, xend =1, yend = 1)+expand_limits(x = 0, y = 0)+scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0))+xlab("predicted probability")+ylab("actual probability")+theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+theme(legend.justification=c(1,0), legend.position=c(1,0)) +geom_line()

OK，本期文章結(jié)束，覺得有用的話多多分享喲。

總結(jié)

以上是生活随笔為你收集整理的R语言手动绘制分类Logistic回归模型的校准曲线（Calibration curve）（3）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。