一个优秀的临床预测模型，不仅要有好的区分度(Discrimination)，同时还应有良好的校准度(Calibration)。区分度是模型拟合效果的指标，代表的是模型区分不同结局的能力，最常用的指标就是是ROC曲线的AUC和C-Statistic。校准度则是模型拟合优度的指标，衡量的是结局实际发生概率和模型预测概率之间的一致性，体现了模型对绝对风险预测的准确性，常用Hosmer-Lemeshow good of fit test来评估。校准曲线（Calibration Curves）是把Hosmer-Lemeshow拟合优度检验结果的可视化，是一种用于评估在不同分位数上预测概率和观测概率一致性的视觉工具，实际上就是实际发生概率与预测发生概率的散点图。
校准曲线的绘制原理并不难，根据模型计算每个观测的预测概率，按预测概率由小到大排列等分成10个组（常用10分位数），然后根据每个组发生结局的模型预测数和实际观测数来计算各组的发生率，分别对应横纵坐标绘制散点图即为校准曲线，Pearson卡方分析结果就是Hosmer-Lemeshow的检验结果。这里需要分组的意义在于，对某一个个体而言，其实际发生概率是无法确定的（结局要么发生要么不发生），需要通过分组实现组内实际发生概率和预测概率的计算。绘制二分类logistic回归模型校准曲线，一般需要两个核心的数据，一个是每个提提实际的观察结果（0/1型），另外一个就是通过模型获得的每个个体的预测概率。SPSS的logistic回归的【选项】中提供Hosmer-Lemeshow拟合优度检验，结果会输出检验结果，同时输出每个组的实际发生频数和预测发生频数，我们完全可以据此计算每个组的实际发生率和预测发生率然后绘制散点图。在R中有很多函数可以实现校准曲线的绘制和Hosmer-Lemeshow检验。绘制校准曲线比较经典的有rms包中的calibrate、val.prob函数。R中的其他一些相关程序包：PredictABEL、predtools、pmcalibration、CalibrationCurves、MachineShop、riskRegression、generalhoslem、gofcat、ResourceSelection。示例来自R语言bestglm包中的数据，该数据我们曾在《R笔记：二分类资料的logistic回归》使用过。数据来自一项南非西开普省心脏病高危地区男性的回顾性样本。响应变量chd表示冠心病发生与否，解释变量有sbp（收缩压）、tobacco（累积烟草量(kg)）、ldl（低密度脂蛋白胆固醇）、adiposity（肥胖）、famhist（家族史）、typea（A型行为）、obesity（肥胖）、alcohol（当前饮酒量）、age（发病年龄）。数据处理是一个复杂的过程，包括缺失值处理、变量转换等很多步骤，此处仅演示有这样一个步骤。library(bestglm);library(leaps)
SAheart<-SAheart[-4] #考虑到adiposity和obesity的高度相关性，本例去掉与obesity高度相关的变量adiposity。我们后面用到了bestglm函数来筛选建模变量，这里需要注意结局变量须在最后一列SAheart$chd<-factor(SAheart$chd,levels=c(0,1),labels=c("NO","YES")) #设置因子。有很多函数其实只支持数值型的观测值，对结局变量赋值标签化反而不一定方便。【2】将数据随机划分为训练集和测试集，并利用训练集筛选建模变量
ss<-sample(nrow(SAheart),nrow(SAheart)*0.7)bestglm(trainset,family=binomial,IC="AIC") #需要首先载入bestglm包：library(bestglm)，数据导入部分已经载入。此次笔记的重点在于绘制校准曲线，变量的筛选不再按部就班，直接采用了最优子集方法来筛选。结果筛选到5个建模变量，我们将以这5个变量与预测变量建立二分类logistic回归模型。
【3】利用基础函数glm构建二分类logistic回归模型，获得训练集和测试集的预测概率
CHD.glm<-glm(chd~age+typea+tobacco+famhist+ldl,family=binomial(link = logit), data=trainset)P.glm.train<- predict(CHD.glm,trainset,type ="response") #结合type参数，获得训练集中每个观测值的预测概率。当然也可以直接通过fitted(CHD.glm)获取P.glm.test<- predict(CHD.glm,testset,type ="response") #测试集观测值的预测概率rms包是回归模型策略（Regression Modeling Strategies）的函数集合，可以绘制列线图、校准曲线的绘制等。该包中各种模型有专用的构建函数，二分类logistic回归用lrm函数。不过rms包中未提供Hosmer-Lemeshow检验，多少有点遗憾。library(rms);library(Hmis)
NM<-datadist(trainset)
options(datadist='NM')
label(trainset$typea)<-"type A" #变量在列线图上显示的名称CHD.pred<-lrm(chd~age+typea+tobacco+famhist+ldl,data=trainset,x=T,y=T)Nom.CHD<-nomogram(CHD.pred,fun=plogis,fun.at=c(.01,.05,seq(.1,.9,by=.2),.95,.99),lp=F,funlabel="CHD Rate")plot(Nom.CHD,xfrac=.25,col.grid=c("red","green"))结果图略。
绘制训练集校准曲线，建议使用calibrate函数
Cal.CHD<-calibrate(CHD.pred,method="boot",B=1000)plot(Cal.CHD,xlab="Predicted Probability",ylab="Observed Probability")#calibrate函数是基于非参数平滑（生存分析中是基于将生存时间的间隔子集）采用的重抽样（bootstrap）或交叉验证（cross-validation）来获得校正偏差（校正过拟合）的观测值和预测值。函数的设置参数非常多，具体可参加帮助文件。在logistic回归中核心参数有三个：fit、method和B。fit指定lrm建立的logistic回归模型；method默认boot可用方法有"crossvalidation"、"boot"、".632"、 "randomization"，可简写为 "cross"、"b"、".6"、"rand"。B参数表示最大抽样次数或者交叉验证折数#plot函数中可以通过xlab、ylab指定横坐标和纵坐标名称；xlim、ylim指定横坐标和纵坐标的范围，比如xlim=c(0,1.0)表示很坐标取值范围在0到1；subtitles用于指定是否做作图中显示校正偏差数据的获得方法及一些基本信息n=323 Mean absolute error=0.016 Mean squared error=0.00027
0.9 Quantile of absolute error=0.021结果显示：模型的校准曲线与参考线相接近，说明预测概率与实际概率的一致性较好。你可能觉得绘制的校准图不够美观，想换一些颜色。calibrate的返回值有很多，其中包含了预测概率（predy）、calibrated.orig值（Apparent）和calibrated.corrected值（Bias-corrected）。可借助lines和abline函数对初始的数据线进行覆盖，这样就可以对结果图进行美化了，也可以将不同模型的校准曲线绘于一图。最后别忘了再用legend对图例做修改。plot(Cal.CHD,xlab="Predicted Probability",ylab="Observed Probability",legend=FALSE)lines(Cal.CHD[,c("predy","calibrated.orig")],type="l",lwd=2,pch=20,col="green")
#利用Cal.CHD返回值中的predy列和calibrated.orig列在既有的校准曲线上绘制进行绘制未进行偏差校正的校准曲线。type指定曲线类型，可用l、p、b、c、o、s、h、n；lwd指定曲线的粗细；pch是点的符号；col指定颜色
lines(Cal.CHD[,c("predy","calibrated.corrected")],type="l",lwd=3,pch=20,col="blue")
abline(0,1,lwd=2,pch=20,col="red")#绘制对角线
legend("bottomright",c("Apparent","Bias-corrected","Ideal"),col=c("green","blue","red"),lty=1,lwd=2,bty = "n",cex=0.8)
绘制测试集校准曲线，建议使用val.prob函数
P.test<-predict(CHD.pred,newdata=testset,type="fitted")CHD.test<- ifelse(testset$chd=="NO", 0, 1)val.prob(P.test,CHD.test,pl=TRUE,xlab="Predicted Probability", ylab="Actual Probability", lim=c(0, 1))#predict.lrm函数可获取二分类或有序多分类logistic回归（lrm和orm）的预测值，object指定模型，newdata指定需要预测的数据集，type指定预测的结果，logistic回归对应的type=c("lp", "fitted", "fitted.ind", "mean", "x", "data.frame", "terms", "cterms", "ccterms", "adjto","adjto.data.frame", "model.frame")#val.prob函数用于验证二分类结局的预测概率，提供了Dxy、C、Brier、R2、D、U、Q、Intercept、Slope、Emax、E90、Eavg等众多指标，也可绘制校准曲线。常用的参数有p指定预测概率；y指定结局指标，需要是数值型变量，本例在导入数据时将结局设置为字符型的因子，因此在先用ifelse对数据集中的字符进行了替换；pl是否绘制校准曲线及可选的一些统计信息，默认为TRUE；logistic.cal是否绘制linear logistic calibration，这个线性逻辑校准曲线具体指的什么笔者尚未弄清楚，从名字上猜测可能是利用logitP做的校准曲线；xlab、,ylab指定横纵坐标的名称；lim指定横纵坐标的范围，默认(0,1)；legendloc和statloc分别指定在校准图中图例和统计信息的位置，若为FALSE则不显示#也可用plot函数可对val.prob函数的对象做图结果中除了显示校准曲线与参考线之间的关系，校准曲线和返回值中还提供了其他一些统计指标。如Dxy（预测概率与观测值的秩相关大小，等于2C-0.5）；C（ROC）即C指数，也就是ROC曲线下面积，等于1+Dxy/2；R2为模型的复相关系数（Nagelkerke-Cox-Snell-Maddala-Magee R-squared index,）；D为区分指数（Discrimination index），等于 (Logistic model LR χ2 - 1)/n，值越大越好；LR χ2 及P值；U为不可靠指数（Unreliability index），越小越好；U指数检验χ2和P值 (H0: intercept=0, slope=1))；Q为品质指数（Quality index），越大越好；Brier是预测值和真实值的均方误差，值越小越好；Intercept为截距；Slope为斜率；Emax为预测值和实际值的最大绝对差值；Eavg为预测值和实际值的平均差值；E90为预测值和真实值差值的90%分位数；S:z和S:p是Spiegelhalter Z-test的Z值和P值，本例P=0.910，表明拟合效果较好。当然，val.prob函数也可以绘制训练集的校准曲线，我们只需要把数据集指向训练集就可以了：
P.train<-predict(CHD.pred,newdata=trainset,type="fitted")
CHD.train<- ifelse(trainset$chd=="NO", 0, 1)val.prob(P.train,CHD.train)【4.2】predtools包：Prediction Model Tools
校正曲线绘制函数：
calibration_plot(data,obs,follow_up=NULL,pred,group=NULL,nTiles=10,legendPosition="right",title=NULL,x_lim=NULL,y_lim=NULL,xlab="Prediction",ylab="Observation",points_col_list=NULL,data_summary=FALSE)#需要特别注意的是obs和pred需要是data数据集的变量，直接引用非data内的数据会报错，所以我们需要首先将pred值写入到数据集中trainset$P.pred<-predict(CHD.glm,type="response")#模型对训练集的预测概率写入训练集testset$P.pred<-predict(CHD.glm,testset,type="response") #模型对测试集的预测概率写入测试集trainset$chd<- ifelse(trainset$chd=="NO", 0, 1) #因变量需要是0/1型的数值型testset$chd<- ifelse(testset$chd=="NO", 0, 1)calibration_plot(data=trainset,obs="chd",pred="P.pred",title="Calibration Curves for Trainset",x_lim=c(0,1),y_lim=c(0,1),xlab="Prediction Probability",ylab="Observation Probability")calibration_plot(data=testset,obs="chd",pred="P.pred",title="Calibration Curves for Testset",xlab="Prediction Probability",ylab="Observation Probability")logistic回归模型可以看做在某个时间点上无删失的二分类结局的生存分析。riskRegression包可以实现竞争风险生存数据的风险回归模型和预测评分分析，其中plotCalibratio函数可绘制生存模型和logistic回归模型的校准曲线。上接第【3】步
在第【3】步我们构建了含有5个预测变量的logistic模型CHD.glm，为演示在同一幅图中绘制多条校准曲线，我们再构建一个含有两个变量的预测模型：CHD.2var<-glm(chd~age+tobacco,family=binomial(link = logit), data=trainset)
scores.train<-Score(object=list("CHD.5var"=CHD.glm,"CHD.2var"=CHD.2var),formula=chd~1,data=trainset,metrics=c("AUC","Brier"),plots=c("ROC", "Calibration", "boxplot"),B=1000) #模型对训练集预测性能评分。formula:Y ~ 1 for binary and Hist(time,status) ~ 1 for time-to-event outcome
scores.test<-Score(object=list("CHD.pred"=CHD.glm),formula=chd~1,data=testset,metrics=c("AUC","Brier"),plots=c("ROC", "Calibration", "boxplot")) #模型对测试集预测性能评分单独绘制训练集的校准曲线
plotCalibration(x=scores.train,models="CHD.5var",percent=FALSE)
将训练集和测试集的校准曲线绘制在同一幅图中，并进行美化
plotCalibration(x=scores.train,models="CHD.5var",method="quantile",q=10,xlim=c(0,1),ylim=c(0,1),xlab="Predicted Risk",ylab="Observed Frequency",auc.in.legend=FALSE,brier.in.legend=FALSE,col="blue",lwd=2,cex=0.8)
#col指定颜色，lwd指定线条宽度，lty指定线条的类型，pch指定点和符号的类型，type控制散点图的类型(p/l/b/o/h/s/n),cex符号大小
plotCalibration(x=scores.test,models="CHD.pred",method="quantile",q=10,add=TRUE,col="green",lwd=2,cex=0.8)
绘制多个模型的校准曲线
plotCalibration(x=scores.train,models=list("CHD.5var","CHD.2var"),col=c("blue","red"),lwd=2,cex=0.8)plotCalibration(x=scores.test,models="CHD.pred",method="nne",bar=TRUE,show.frequencies=TRUE,col=c("tomato","blueviolet"))
横坐标表示预测概率，按大小分为了10个风险段（q=10），纵坐标表示发生频率。第一个风险段组风险大小1.9%-6.8%，预测概率7%，实际频率4%；第十个风险段组风险大小68.9%-84.6%，预测概率79%，实际频率77%。整体来说，预测概率是实际概率的一致性还是不错的。【4.4】PredictABEL包：Assessment of Risk Prediction Models上接第【2】步
CHD.lgr<-fitLogRegModel(data=trainset,cOutcome=9,cNonGenPreds=c(2:3,5,8),cNonGenPredsCat=4,cGenPreds=c(0),cGenPredsCat=c(0))#此函数实际上就是利用标准的GLM函数拟合了一个二分类logistic回归，结果同glm函数的结果【CHD.glm<-glm(chd~age+typea+tobacco+famhist+ldl,,family=binomial(link = logit), data=trainset)】是一致的。data需含有结局变量和预测变量；cOutcome指定结局变量的列数，本例第9列为结局变量。另外, 结局变量须是0/1型，1为感兴趣的结局；cNonGenPreds指定连续型非基因预测变量的列数，本例为第2、3、5、8列；cNonGenPredsCat指定非基因的分类变量，本例是第4列；cGenPreds和cGenPredsCat分别指定连续型和分类型的基因变量，示例数据均不是基因的变量，使用c(0)表示没有P.risk.train<-predRisk(riskModel=CHD.lgr, data=trainset)#获得每个个体的预测风险。其实就是每个个体的预测概率【P.glm.train<- predict(CHD.glm,trainset,type ="response")或fitted(CHD.glm)】P.risk.test<-predRisk(riskModel=CHD.lgr, data=testset)trainset$chd<- ifelse(trainset$chd=="NO", 0, 1) #因变量需要是0/1型的数值型testset$chd<- ifelse(testset$chd=="NO", 0, 1)plotCalibration(data=trainset, cOutcome=9, predRisk=P.risk.train, groups=10, rangeaxis=c(0,1), plottitle="Calibration plot", xlabel="Predicted risk", ylabel="Observed risk",filename="calibration.txt",fileplot="plotname", plottype="jpg")#data指定数据集；cOutcome指定结局变量的列数；predRisk为预测概率, groups分位数分组组数；rangeaxis横纵坐标范围；plottitle绘图名称；xlabel、ylabel横纵坐标名称；filename指定在工作目录中对校准数据的表格进行保存的名称，不指定则不保存。将该文件重新导入到R中，可以在既有的校准图中通过abline函数绘制拟合曲线，但目前保存的text文件将数据的小数点保存成了逗号，导致数据无法使用，期望后续版本的修订。当然一般分10组，全部的观测概率和预测概率也不过20个，手动输入数据也并不是很麻烦；fileplot指定在工作目录中保存校准图的名称，不指定则不保存；plottype指定保存校准图的理性，默认格式是epsplotCalibration(data=testset, cOutcome=9, predRisk=P.risk.test) 执行命令后不仅有校准曲线，还输出了Hosmer-Lemeshow拟合优度检验结果。结果显示各个数据点均在理想的校准曲线附近波动，Hosmer-Lemeshow校验（训练集：chi2=5.978，P=0.6497；测试集：chi2=8.744，P=0.3644）模型对绝对风险的估计与实际发生风险一致性较好。
除了PredictABEL包，其他程序包诸如generalhoslem、gofcat、ResourceSelection等也提供一些函数可以完成HL的检验。【5】Hosmer-Lemeshow拟合优度检验
如果Hosmer-Lemeshow拟合优度检验P>0.05，则表明预测概率与实际概率一致性较好。【5.1】logitgof{generalhoslem}library(reshape);library(MASS)训练集HL检验
logitgof(obs=trainset$chd,exp=fitted(CHD.glm),g=10,ord=FALSE) #obs为观测值；exp为期望值,此处指模型的拟合值。在logistic回归中，拟合值为概率，可以利用predict结合type参数，训练集中也可以直接使用fitted值；g为分数数分组，默认为10；ord是否运行有序版本，当obs值为有序多分类时可改用TRUE# logitgof(obs=trainset$chd,exp=P.train,g=10,ord=FALSE) #利用rms包中predict.lrm函数获的预测概率Hosmer and Lemeshow test (binary model)data: trainset$chd, P.trainX-squared = 5.8379, df = 8, p-value = 0.6654测试集HL检验
logitgof(obs=testset$chd,exp=P.glm.test,g=10,ord=FALSE) #logitgof(obs=testset$chd,exp=P.test,g=10,ord=FALSE) #rms包中predict.lrm函数获的预测概率Hosmer and Lemeshow test (binary model)
data: testset$chd, P.testX-squared = 7.3966, df = 8, p-value = 0.4945In logitgof(obs = testset$chd, exp = P.test, g = 10, ord = FALSE) :At least one cell in the expected frequencies table is < 1. Chi-square approximation may be incorrect.结果中一个警示信息：至少有一个单元格的期望频数<1，卡方分析进行了校正。你可以减少分组数（默认10）来去掉这个警示。
【5.2】hoslem.test{ResourceSelection}CHD.train<- ifelse(trainset$chd=="NO", 0, 1)CHD.test<- ifelse(testset$chd=="NO", 0, 1)library(ResourceSelection)ResourceSelection 0.3-6 2023-06-27hoslem.test(CHD.train,P.glm.train,g=10) #该函数观测值要求为0/1型二分类变量，且需要为数值型。本例在导入数据时将结局设置为字符型的因子需要在转换一下@_@
Hosmer and Lemeshow goodness of fit (GOF) testdata: CHD.train, P.glm.trainX-squared = 5.8379, df = 8, p-value = 0.6654hoslem.test(CHD.test,P.glm.test,g=10)
Hosmer and Lemeshow goodness of fit (GOF) testdata: CHD.test, P.glm.testX-squared = 7.3966, df = 8, p-value = 0.4945关注下方公众号，分享更多更好玩的R语言知识。
觉得有帮助的请点赞、分享、在看走起！
点个在看，SCI马上发表。