從大家好，我是風(fēng)間琉璃。好久不見(jiàn)甚是想念，之前的坑拖著一直沒(méi)補(bǔ)上十分愧疚。今天我們繼續(xù)好吧。 之前的坑大家請(qǐng)看“生信高分文章delogistic臨床模型如何評(píng)價(jià)？看完這篇就明白了（上）”。

上回我們講到通過(guò)年齡、性別和指標(biāo)A來(lái)對(duì)高血壓（hypertension）的發(fā)生進(jìn)行預(yù)測(cè)并建模。復(fù)習(xí)一下，上一次我們講了怎么建模以及做好看的ROC，接下我們要講Calibration和DCA了。大家注意聽(tīng)喔。首先重申一下對(duì)于模型的Calibration，其主要目的是評(píng)估我們預(yù)測(cè)得到的預(yù)測(cè)值與實(shí)際的真實(shí)性的一致性。

通俗的說(shuō)就是我們的到的預(yù)測(cè)值和真實(shí)值存在多少偏差，這個(gè)偏差能不能接受？

這一步我們分為兩個(gè)小步，第一步就是繪制我們熟悉的calibration圖（我將給大家介紹兩個(gè)方法），第二步Hosmer-Lemeshow good of fit test（擬合優(yōu)度檢驗(yàn)）來(lái)進(jìn)行相關(guān)的評(píng)估得到P值。

Ok我們先開(kāi)始，加載我們需要的包rms包的功能非常強(qiáng)大,里面的兩個(gè)函數(shù)功能都能達(dá)到矯正的作用，但是它的方法是不一樣的。我們先上數(shù)據(jù)。

###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對(duì)模型進(jìn)行評(píng)價(jià)   ###畫(huà)圖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])##得到預(yù)測(cè)值，由于是logistics回歸，故需要進(jìn)行l(wèi)ogit轉(zhuǎn)換y_pred <- 1/(1+exp(-pred.logit))#####val.prob進(jìn)行校正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的內(nèi)置功能建模，也就是lrm進(jìn)行建模才行。R本身自帶的glm不能達(dá)到效果。好，接下來(lái)我們來(lái)看看兩種方式的不同。

第一種

我們首先看calibration里面有一個(gè)boot，這個(gè)其實(shí)就是代表這個(gè)功能進(jìn)行calibration的核心思想，也就是重抽樣的思想。說(shuō)起這個(gè)也是一個(gè)趣事，重抽樣（bootstrap）本身的思想很簡(jiǎn)單，它和交叉驗(yàn)證（cross-validation）是機(jī)器學(xué)習(xí)非常重要的兩大思想。重抽樣可以理解為，你想知道池塘里面有多少條魚(yú)，你先打撈上來(lái)100條，做好標(biāo)記后放回去。然后然撈上來(lái)100條，根據(jù)這100條里有多少是你標(biāo)記的魚(yú)來(lái)推測(cè)總量，如此反復(fù)后不就知道這個(gè)池塘有多少條魚(yú)了嗎。

之前我說(shuō)的趣事也就是當(dāng)年這個(gè)思想的開(kāi)發(fā)者寫(xiě)好論文向統(tǒng)計(jì)學(xué)的某某大牛期刊投稿時(shí)，期刊編輯給出了拒絕，并說(shuō)“想法很好，但過(guò)于簡(jiǎn)單”。但是后面的實(shí)踐發(fā)現(xiàn)bootstrap的作用遠(yuǎn)超大家的想象?，F(xiàn)在想想，大牛發(fā)現(xiàn)這么有創(chuàng)造性的方法都能被拒，我們被拒幾次豈不是稀松平常？哈哈哈哈，不吹牛我們繼續(xù)。接下來(lái)我們看一看圖的X軸和Y軸，X軸代表我們預(yù)測(cè)的可能性，而Y軸代表的是實(shí)際的可能性。理想狀態(tài)我們的Calibration曲線應(yīng)該是斜率為1從原點(diǎn)穿過(guò)的直線，也就是圖中的虛線。我們可以看到我們的黑線，也就是實(shí)際的線和虛線差不太多，證明我們的結(jié)果還是可以的哈。并且平均的總的誤插在0.024，還不錯(cuò)！

第二種

val.prob，其實(shí)就是我們傳統(tǒng)認(rèn)知的calibration的方法，根據(jù)我們預(yù)測(cè)出來(lái)的預(yù)測(cè)概率與真實(shí)值進(jìn)行矯正。我們看到這個(gè)val.prob和calibate功能不同的是它出現(xiàn)了很多其他奇奇怪怪的指標(biāo)。我們主要關(guān)注四個(gè)指標(biāo)Dxy, C(ROC), U: p, Brier.Dxy即是我們的預(yù)測(cè)值和真實(shí)值之間的相關(guān)性，C(ROC)即是ROC曲線下面積。U檢驗(yàn)本身是unrealiability檢驗(yàn)，即是假設(shè)預(yù)測(cè)值和真實(shí)值兩者之間沒(méi)有相關(guān)性。我們發(fā)現(xiàn)U:p很小，接近于0。那么則提示我們有很好的校正度。最后一個(gè)brier即是預(yù)測(cè)值和真實(shí)值的平均平方差，這個(gè)值當(dāng)然越小越好。如果大家還想知道更多。像我一樣help一下

help("calibrate")## starting httpd help server ... donehelp("val.prob")

###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弄完了，接下來(lái)我們要到最后的部分，也是最加分的部分，DCA。在上一篇我大概講了什么是DCA，因?yàn)槔碚撎啵覀兿葘W(xué)會(huì)用，理論我們慢慢的學(xué)都可以。首先我們還是得加載包，rmda包是我們依賴(lài)得做出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構(gòu)建模型得到曲線，如果我們的紅線在水平的黑線以及左側(cè)斜向的灰線以上，那么我們可以認(rèn)為這一段的紅線取值是能夠獲益的。而我們的模型，紅線從幾乎0到1都在灰線和黑線之上，我們可以認(rèn)為根據(jù)我們所構(gòu)建模型所作出得決策能夠給病人帶來(lái)獲益，因此這是一個(gè)相當(dāng)完美得模型。

大家祝賀一下！希望我們自己文章的模型也能有這么完美?。?！那么整體的內(nèi)容我們已經(jīng)弄完了，不過(guò)大家覺(jué)得是不是差點(diǎn)意思？大家覺(jué)得還可以，不差嗎？但我覺(jué)得可能還差點(diǎn)意思。來(lái)，我們補(bǔ)上。

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)

這個(gè)圖大家是不是很熟悉，對(duì)，就是我們得nomogram。我們每一個(gè)變量根據(jù)它的大小我們都能得到一個(gè)評(píng)分，并且根據(jù)每個(gè)變量進(jìn)行相加。大家如果想知道更多，請(qǐng)看我們之前的文章——Nomogram圖不會(huì)畫(huà)？看了這篇，小白也能輕松搞定。

OK,我們的logistic模型的評(píng)價(jià)就做完了。謝謝大家??！