以下內容來自B站up主：大鵬統(tǒng)計工作室 的系列教學視頻《R語言Logistic回歸臨床預測模型》
第七節(jié) calibration校準曲線(4種款式）

假定已經擬合好了下面的模型：

image.png

dd = datadist(mydata)
options(datadist = “dd”)

Method 1

fit1 <- lrm(formula1, data = mydata, x=TRUE,y=TRUE)
cal1 <- calibrate(fit1, method = "boot", B=1000)

plot(cal1,
     xlim = c(0,1),
     xlab = "Predicted Probability",
     ylab = "Observed Probability",
     legend = FALSE,
     subtitles = FALSE)
abline(0,1,col = "black",lty = 2,lwd = 2)
lines(cal1[,c("predy","calibrated.orig")], type = "l",lwd = 2,col="red",pch =16)
lines(cal1[,c("predy","calibrated.corrected")], type = "l",lwd = 2,col="green",pch =16)
legend(0.55,0.35,
       c("Apparent","Ideal","Bias-corrected"),
       lty = c(2,1,1),
       lwd = c(2,1,1),
       col = c("black","red","green"),
       bty = "n") # "o"為加邊框

image.png

Method 2

library(riskRegression)

fit2 <- glm(formula2,data = mydata,family = binomial())
xb <- Score(list("fit"=fit2),formula= Group~1,
            null.model = FALSE,
            conf.int = TRUE,
            plots = c("calibration","ROC"),
            metrics = c("auc","brier"),
            B=1000,M = 50, # bootstrap；每次抽樣50
            data= mydata)
plotCalibration(xb,col="red")

image.png

Method 3(多個模型結果匯總1)

fit1 = glm(formula1, data=mydata,family = binomial())
fit2 = glm(formula2, data=mydata,family = binomial())
fit3 = glm(formula3, data=mydata,family = binomial())
fit4 = glm(formula4, data=mydata,family = binomial())
xb <- Score(list("fit1"=fit1,
                 "fit2"=fit2,
                 "fit3"=fit3,
                 "fit4"=fit4),
            formula=Group~1,
            null.model = FALSE,
            conf.int =TRUE,
            plots =c("calibration","ROC"),
            metrics = c("auc","brier"),
            B=1000,M=50,
            data=mydata)
plotCalibration(xb)

image.png

Method 4(多個模型結果匯總2)

fit1 <- lrm(formula1, data = mydata, x=TRUE,y=TRUE)
cal1 <- calibrate(fit1, method = "boot", B=1000)

fit2 <- lrm(formula2, data = mydata, x=TRUE,y=TRUE)
cal2 <- calibrate(fit2, method = "boot", B=1000)

fit3 <- lrm(formula3, data = mydata, x=TRUE,y=TRUE)
cal3 <- calibrate(fit3, method = "boot", B=1000)

fit4 <- lrm(formula4, data = mydata, x=TRUE,y=TRUE)
cal4 <- calibrate(fit4, method = "boot", B=1000)

plot(1,type ="n",
     xlim =c(0,1),
     ylim =c(0,1),
     xlab = "Predicted Probability",
     ylab = "Observed Probability",
     legend = FALSE,
     subtitles = FALSE)
abline(0,1,col = "black",lty = 2,lwd = 2)
lines(cal1[,c("predy","calibrated.corrected")], type = "l",lwd = 2,col="red",pch =16)
lines(cal2[,c("predy","calibrated.corrected")], type = "l",lwd = 2,col="green",pch =16)
lines(cal3[,c("predy","calibrated.corrected")], type = "l",lwd = 2,col="blue",pch =16)
lines(cal4[,c("predy","calibrated.corrected")], type = "l",lwd = 2,col="yellow",pch =16)

legend(0.55,0.35,
       c("fit1","fit2","fit3","fit4"),
       lty = c(2,2,2,2),
       lwd = c(2,2,2,2),
       col = c("red","green","blue","yellow"), #和上面的保持一致
       bty = "n") # "o"為加邊框

image.png