|
上回我们讲到通过年龄、性别和指标A来对高血压（hypertension）的发生进行预测并建模。复习一下，上一次我们讲了怎么建模以及做好看的ROC，接下我们要讲Calibration和DCA了。大家注意听喔。首先重申一下对于模型的Calibration，其主要目的是评估我们预测得到的预测值与实际的真实性的一致性。 通俗的说就是我们的到的预测值和真实值存在多少偏差，这个偏差能不能接受？ 这一步我们分为两个小步，第一步就是绘制我们熟悉的calibration图（我将给大家介绍两个方法），第二步Hosmer-Lemeshow good of fit test（拟合优度检验）来进行相关的评估得到P值。 Ok我们先开始，加载我们需要的包rms包的功能非常强大,里面的两个函数功能都能达到矫正的作用，但是它的方法是不一样的。我们先上数据。 ###data processingdataset = read.csv('disease.csv')dataset = dataset[2:5]dataset$hypertension = factor(dataset$hypertension, levels = c(0, 1))library(caTools)set.seed(123)split = sample.split(dataset$hypertension, SplitRatio = 0.7)training_set = subset(dataset, split == TRUE)test_set = subset(dataset, split == FALSE)###construct the modeloptions(datadist=NULL)classifier=lrm(hypertension~.,data=training_set,x=T,y=T)###method 1calib<-calibrate(classifier,method = "boot") ##calIbrate对模型进行评价 ###画图plot(calib,lwd=2,lty=1, xlim=c(0,1),ylim = c(0,1), xlab = "Predicted Probability of hypertension", ylab = "Actual probability", col=c(rgb(192,98,83,maxColorValue = 255)))
#### n=280 Mean absolute error=0.024 Mean squared error=0.00088## 0.9 Quantile of absolute error=0.053###method 2pred.logit = predict(classifier, newdata = test_set[-4])##得到预测值，由于是logistics回归，故需要进行logit转换y_pred <- 1/(1+exp(-pred.logit))#####val.prob进行校正val.prob(y_pred,as.numeric(test_set$hypertension),m=20,cex=0.5)
## Dxy C (ROC) R2 D D:Chi-sq## 8.205980e-01 9.102990e-01 5.492686e-01 5.030060e-01 6.136072e+01## D:p U U:Chi-sq U:p Q## NA 2.532954e+00 3.059545e+02 0.000000e+00 -2.029948e+00## Brier Intercept Slope Emax E90## 1.138589e+00 -2.709633e-02 7.785694e-01 1.135061e+00 1.094938e+00## Eavg S:z S:p## 1.000890e+00 2.083508e+01 2.081675e-96
好，我们使用两种方式都能得到calibration的图。提一句rms里面的这两种功能都需要使用rms的内置功能建模，也就是lrm进行建模才行。R本身自带的glm不能达到效果。好，接下来我们来看看两种方式的不同。 第一种 我们首先看calibration里面有一个boot，这个其实就是代表这个功能进行calibration的核心思想，也就是重抽样的思想。说起这个也是一个趣事，重抽样（bootstrap）本身的思想很简单，它和交叉验证（cross-validation）是机器学习非常重要的两大思想。重抽样可以理解为，你想知道池塘里面有多少条鱼，你先打捞上来100条，做好标记后放回去。然后然捞上来100条，根据这100条里有多少是你标记的鱼来推测总量，如此反复后不就知道这个池塘有多少条鱼了吗。 之前我说的趣事也就是当年这个思想的开发者写好论文向统计学的某某大牛期刊投稿时，期刊编辑给出了拒绝，并说“想法很好，但过于简单”。但是后面的实践发现bootstrap的作用远超大家的想象。现在想想，大牛发现这么有创造性的方法都能被拒，我们被拒几次岂不是稀松平常？哈哈哈哈，不吹牛我们继续。接下来我们看一看图的X轴和Y轴，X轴代表我们预测的可能性，而Y轴代表的是实际的可能性。理想状态我们的Calibration曲线应该是斜率为1从原点穿过的直线，也就是图中的虚线。我们可以看到我们的黑线，也就是实际的线和虚线差不太多，证明我们的结果还是可以的哈。并且平均的总的误插在0.024，还不错！ 第二种 val.prob，其实就是我们传统认知的calibration的方法，根据我们预测出来的预测概率与真实值进行矫正。我们看到这个val.prob和calibate功能不同的是它出现了很多其他奇奇怪怪的指标。我们主要关注四个指标Dxy, C(ROC), U: p, Brier.Dxy即是我们的预测值和真实值之间的相关性，C(ROC)即是ROC曲线下面积。U检验本身是unrealiability检验，即是假设预测值和真实值两者之间没有相关性。我们发现U:p很小，接近于0。那么则提示我们有很好的校正度。最后一个brier即是预测值和真实值的平均平方差，这个值当然越小越好。如果大家还想知道更多。像我一样help一下 help("calibrate")## starting httpd help server ... donehelp("val.prob")
第二步：Hosmer-Lemeshow检验得到校准P值 Hosmer-Lemeshow检验得到校准P值，其实HL检验是目前统计上用的比较多的检验校准度的方法。其思想是：1.首先根据预测模型来计算每个个体未来发生结局事件的预测概率； 2.根据预测概率从小到大进行排序，并按照十分位等分成10组； 3. 分别计算各组的实际观测数和模型预测数，其中模型预测数，即每个人的预测概率*人数，再求总和，这里人数即为1，最后总和就相当于每个个体预测概率的直接加和； 4. 根据每组实际观测数和模型预测数计算卡方值（自由度=8），再根据卡方分布得到对应的P值。 大家看懂了吗？没事就算没看懂也不耽误我们来用是吧。我们直接看结果。如果所得的统计量卡方值越小，对应的P值越大，则提示预测模型的校准度越好。若检验结果显示有统计学显著性（P<0.05），则表明模型预测值和实际观测值之间存在一定的差异，模型校准度差。好理论有了，上代码。###install.packages("RsourceSelection")library(ResourceSelection)## ResourceSelection 0.3-5 2019-07-22hoslem.test(as.numeric(test_set$hypertension),y_pred,g=10)#### Hosmer and Lemeshow goodness of fit (GOF) test#### data: as.numeric(test_set$hypertension), y_pred## X-squared = 6533.2, df = 8, p-value < 2.2e-16
好了，我们的Calibration弄完了，接下来我们要到最后的部分，也是最加分的部分，DCA。在上一篇我大概讲了什么是DCA，因为理论太多，我们先学会用，理论我们慢慢的学都可以。首先我们还是得加载包，rmda包是我们依赖得做出DCA得包，所以我们必须得加载它。 #install.packages(“rmda”)library(rmda)training_set[,4]=as.numeric(training_set[,4])-1dca=decision_curve(formula = hypertension ~Gender+Age+A, data = training_set,family = "binomial", confidence.intervals = F,bootstraps = 10,fitted.risk = F)## Note: The data provided is used to both fit a prediction model and to estimate the respective decision curve. This may cause bias in decision curve estimates leading to over-confidence in model performance.plot_decision_curve(dca,curve.names = "dca model")
红色得曲线即是我们DCA构建模型得到曲线，如果我们的红线在水平的黑线以及左侧斜向的灰线以上，那么我们可以认为这一段的红线取值是能够获益的。而我们的模型，红线从几乎0到1都在灰线和黑线之上，我们可以认为根据我们所构建模型所作出得决策能够给病人带来获益，因此这是一个相当完美得模型。 大家祝贺一下！希望我们自己文章的模型也能有这么完美！！！那么整体的内容我们已经弄完了，不过大家觉得是不是差点意思？大家觉得还可以，不差吗？但我觉得可能还差点意思。来，我们补上。 ddist <- datadist(training_set)options(datadist='ddist')nomogram1 <- nomogram(classifier, fun=function(x)1/(1+exp(-x)), # or fun=plogis fun.at=seq(0.3,0.8,by=0.1), funlabel="Risk of hypertension")plot(nomogram1)
这个图大家是不是很熟悉，对，就是我们得nomogram。我们每一个变量根据它的大小我们都能得到一个评分，并且根据每个变量进行相加。大家如果想知道更多，请看我们之前的文章——Nomogram图不会画？看了这篇，小白也能轻松搞定。 OK,我们的logistic模型的评价就做完了。谢谢大家！！
|
