【原】R語言KNN（K最近鄰、ROC、PR曲線、校準(zhǔn)曲線、決策曲線）

阿越就是我 2024-05-13 發(fā)布于上海

展開全文

本文主要介紹如何使用R語言實(shí)現(xiàn)K最近鄰算法。

本文目錄：

算法簡(jiǎn)介
準(zhǔn)備數(shù)據(jù)
class包

建立模型
模型評(píng)價(jià)
超參數(shù)調(diào)優(yōu)

kknn包

建立模型
模型評(píng)價(jià)
超參數(shù)調(diào)優(yōu)

算法簡(jiǎn)介

K最近鄰（K-Nearest-Neighbor，KNN）是一種非線性的分類算法，KNN處理分類問題的方法是：找K個(gè)距離待預(yù)測(cè)樣本最近的點(diǎn)，然后根據(jù)這幾個(gè)點(diǎn)的類別來確定新樣本的類別。

近朱者赤，近墨者黑。

下面以一個(gè)二分類問題為例說明KNN的思想。

下圖有兩個(gè)特征可以用來預(yù)測(cè)腫瘤是”良性”還是”惡性”。圖中的X表示我們要預(yù)測(cè)的新樣本。如果算法設(shè)定k=3，那么圓圈中包含的3個(gè)觀測(cè)就是樣本X的最近鄰。因?yàn)槠渲姓级鄶?shù)比例的類別是”惡性”，所以樣本X被分類為”惡性”。

思想是不是很簡(jiǎn)單？K的選擇對(duì)于KNN的預(yù)測(cè)結(jié)果是非常重要的。

KNN中另一個(gè)需要指出的重要問題是距離的計(jì)算方法，或者說是特征空間中數(shù)據(jù)點(diǎn)的臨近度的計(jì)算。默認(rèn)的距離是歐氏距離，也就是從點(diǎn)A到點(diǎn)B的簡(jiǎn)單直線距離。

Note
兩點(diǎn)間的距離強(qiáng)烈依賴于測(cè)量特征時(shí)使用的單位，所以必須對(duì)其進(jìn)行標(biāo)準(zhǔn)化，而且要求數(shù)據(jù)不能有缺失值。

準(zhǔn)備數(shù)據(jù)

演示數(shù)據(jù)為印第安人糖尿病數(shù)據(jù)集，這個(gè)數(shù)據(jù)一共有768行，9列，其中diabetes是結(jié)果變量，為二分類，其余列是預(yù)測(cè)變量。

該數(shù)據(jù)集的原始版本是有缺失值的，我這里使用的是插補(bǔ)過的版本，詳細(xì)過程請(qǐng)參考數(shù)據(jù)準(zhǔn)備這一章。

rm(list = ls())

load(file = "../000機(jī)器學(xué)習(xí)/pimadiabetes.rdata")

dim(pimadiabetes)
## [1] 768   9
str(pimadiabetes)
## 'data.frame':    768 obs. of  9 variables:
##  $ pregnant: num  6 1 8 1 0 5 3 10 2 8 ...
##  $ glucose : num  148 85 183 89 137 116 78 115 197 125 ...
##  $ pressure: num  72 66 64 66 40 ...
##  $ triceps : num  35 29 22.9 23 35 ...
##  $ insulin : num  202.2 64.6 217.1 94 168 ...
##  $ mass    : num  33.6 26.6 23.3 28.1 43.1 ...
##  $ pedigree: num  0.627 0.351 0.672 0.167 2.288 ...
##  $ age     : num  50 31 32 21 33 30 26 29 53 54 ...
##  $ diabetes: Factor w/ 2 levels "pos","neg": 2 1 2 1 2 1 2 1 2 2 ...

各個(gè)變量的含義：

pregnant：懷孕次數(shù)
glucose：血漿葡萄糖濃度（葡萄糖耐量試驗(yàn)）
pressure：舒張壓（毫米汞柱）
triceps：三頭肌皮褶厚度（mm）
insulin：2小時(shí)血清胰島素（mu U/ml）
mass：BMI
pedigree：糖尿病譜系功能，是一種用于預(yù)測(cè)糖尿病發(fā)病風(fēng)險(xiǎn)的指標(biāo)，該指標(biāo)是基于家族史的糖尿病遺傳風(fēng)險(xiǎn)因素的計(jì)算得出的。它計(jì)算了患者的家族成員是否患有糖尿病以及他們與患者的親緣關(guān)系，從而得出一個(gè)綜合評(píng)分，用于預(yù)測(cè)患糖尿病的概率。
age：年齡
diabetes：是否有糖尿病

先對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化：

# 對(duì)數(shù)值型變量進(jìn)行標(biāo)準(zhǔn)化
pimadiabetes[,-9] <- scale(pimadiabetes[,-9])
str(pimadiabetes)
## 'data.frame':    768 obs. of  9 variables:
##  $ pregnant: num  0.64 -0.844 1.233 -0.844 -1.141 ...
##  $ glucose : num  0.863 -1.203 2.011 -1.072 0.503 ...
##  $ pressure: num  -0.0314 -0.5244 -0.6887 -0.5244 -2.6607 ...
##  $ triceps : num  0.63124 -0.00231 -0.64853 -0.63586 0.63124 ...
##  $ insulin : num  0.478 -0.933 0.63 -0.631 0.127 ...
##  $ mass    : num  0.172 -0.844 -1.323 -0.626 1.551 ...
##  $ pedigree: num  0.468 -0.365 0.604 -0.92 5.481 ...
##  $ age     : num  1.4251 -0.1905 -0.1055 -1.0409 -0.0205 ...
##  $ diabetes: Factor w/ 2 levels "pos","neg": 2 1 2 1 2 1 2 1 2 2 ...

class包

數(shù)據(jù)劃分為訓(xùn)練集和測(cè)試集，劃分比例為7：3。

但是R語言里class包在使用時(shí)需要把真實(shí)結(jié)果去掉，所以我們把真實(shí)結(jié)果去掉，只保留預(yù)測(cè)變量。

# 劃分是隨機(jī)的，設(shè)置種子數(shù)可以讓結(jié)果復(fù)現(xiàn)
set.seed(123)
ind <- sample(1:nrow(pimadiabetes), size = 0.7*nrow(pimadiabetes))

# 去掉真實(shí)結(jié)果列
train <- pimadiabetes[ind,-9]
test <- pimadiabetes[-ind, -9]

dim(train)
## [1] 537   8
dim(test)
## [1] 231   8

str(train)
## 'data.frame':    537 obs. of  8 variables:
##  $ pregnant: num  -1.141 1.233 0.343 -0.251 1.233 ...
##  $ glucose : num  0.535 -1.564 0.699 -1.137 -1.203 ...
##  $ pressure: num  -1.017 -0.196 0.462 -1.017 -1.428 ...
##  $ triceps : num  0.631 1.159 0.86 -1.164 -0.953 ...
##  $ insulin : num  0.117 -1.092 1.093 -0.881 -0.767 ...
##  $ mass    : num  0.317 0.419 1.827 -1.541 -1.164 ...
##  $ pedigree: num  0.1875 0.7036 -0.8507 -0.0841 -1.0137 ...
##  $ age     : num  -1.041 0.49 1.17 -1.041 0.745 ...

# 把真實(shí)結(jié)果列單獨(dú)拿出來，后面用
truth_train <- pimadiabetes[ind,9]
truth_test <- pimadiabetes[-ind,9]

建立模型

在訓(xùn)練集建立模型，1行代碼搞定：

library(class)

f <- knn(train = train, # 訓(xùn)練集，只有預(yù)測(cè)變量，沒有結(jié)果變量
         test = test, # 測(cè)試集，沒有結(jié)果變量
         cl = truth_train, # 訓(xùn)練集的真實(shí)結(jié)果
         k = 8, # 使用的近鄰個(gè)數(shù)
         prob = TRUE # 需要計(jì)算概率
         )

# 查看測(cè)試集的預(yù)測(cè)結(jié)果，只看前6個(gè)
head(f)
## [1] neg neg pos neg neg pos
## Levels: pos neg

# 查看測(cè)試集的預(yù)測(cè)概率，只看前6個(gè)
prob <- attr(f,"prob")
head(prob)
## [1] 0.500 0.750 1.000 0.750 0.625 0.625

此時(shí)得到的f這個(gè)結(jié)果是一個(gè)因子型的向量，而且是有名字和屬性的，大多數(shù)模型擬合結(jié)果的格式都是不一樣的，使用時(shí)需要注意！

還要注意這里的概率，并不是陽性結(jié)果的概率，而是預(yù)測(cè)結(jié)果的概率！比如第一個(gè)概率0.500是neg的概率，第二個(gè)概率0.750是neg的概率，第三個(gè)概率1.000是pos的概率！

所以你如果想要陽性結(jié)果(pos)的概率，需要自己計(jì)算一下：

prob <- ifelse(f == "pos", prob, 1-prob)
head(prob)
## [1] 0.500 0.250 1.000 0.250 0.375 0.625

模型評(píng)價(jià)

模型擬合好之后，下一步就是查看模型的各種指標(biāo)，來看看這個(gè)模型在訓(xùn)練集中的表現(xiàn)如何，比如混淆矩陣、AUC值、準(zhǔn)確率等。

不管是什么類型的模型，如果我們想要評(píng)價(jià)它的模型表現(xiàn)，都是需要用到模型的預(yù)測(cè)結(jié)果和真實(shí)結(jié)果的。

對(duì)于回歸任務(wù)來說，預(yù)測(cè)結(jié)果也是數(shù)值型的，對(duì)于分類任務(wù)來說，模型的預(yù)測(cè)結(jié)果可以是某一種類別的概率，也可以是預(yù)測(cè)出的具體類別。大多數(shù)模型都是既支持計(jì)算類別概率又支持計(jì)算具體類別的，但是有些模型可能只支持一種類型。

首先查看混淆矩陣，我們借助caret包展示，這個(gè)包目前仍然是查看混線矩陣最全面的，，沒有之一，非常好用！

# 借助caret包，f是預(yù)測(cè)的類別，truth_test是真實(shí)的結(jié)果
caret::confusionMatrix(f,truth_test)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction pos neg
##        pos 131  32
##        neg  19  49
##                                          
##                Accuracy : 0.7792         
##                  95% CI : (0.7201, 0.831)
##     No Information Rate : 0.6494         
##     P-Value [Acc > NIR] : 1.269e-05      
##                                          
##                   Kappa : 0.4966         
##                                          
##  Mcnemar's Test P-Value : 0.09289        
##                                          
##             Sensitivity : 0.8733         
##             Specificity : 0.6049         
##          Pos Pred Value : 0.8037         
##          Neg Pred Value : 0.7206         
##              Prevalence : 0.6494         
##          Detection Rate : 0.5671         
##    Detection Prevalence : 0.7056         
##       Balanced Accuracy : 0.7391         
##                                          
##        'Positive' Class : pos            
##

結(jié)果非常全面，最上面是混淆矩陣，然后給出了：

準(zhǔn)確率（Accuracy）和準(zhǔn)確率的可信區(qū)間（95% CI）
無信息率（No Information Rate）和P值（P-Value Acc>NIR]）、
Kappa值（Kappa一致性指數(shù)）
Mcnemar檢驗(yàn)的P值（Mcnemar’s Test P-Value）
敏感性、特異性
陽性預(yù)測(cè)值、陰性預(yù)測(cè)值
流行率（Prevalence）
檢出率（Detection Rate）
（Detection Prevalence）
均衡準(zhǔn)確率（Balanced Accuracy）

最后告訴你參考類別是pos。

當(dāng)然這些值你也可以單獨(dú)計(jì)算：

caret::precision(f,truth_test) # 精準(zhǔn)率
## [1] 0.803681
caret::recall(f,truth_test) # 召回率，靈敏度
## [1] 0.8733333
caret::F_meas(f,truth_test) # F1分?jǐn)?shù)
## [1] 0.8370607

然后再畫個(gè)ROC曲線。首先使用ROCR進(jìn)行演示，不管是什么包，都遵循前面說過的規(guī)律，畫ROC曲線是需要真實(shí)結(jié)果和預(yù)測(cè)概率的！

library(ROCR)

# ROCR畫ROC曲線就是2步，先prediction，再performance
pred <- prediction(prob,truth_test) # 預(yù)測(cè)概率，真實(shí)類別
perf <- performance(pred, "tpr","fpr") # ROC曲線的橫縱坐標(biāo)，不要寫錯(cuò)了
auc <- round(performance(pred, "auc")@y.values[[1]],digits = 4) # 提取AUC值
auc
## [1] 0.8477

# 畫圖
plot(perf,lwd=2,col="tomato")
abline(0,1,lty=2) # 添加對(duì)角線
# 添加圖例
legend("bottomright", legend=paste("AUC: ",auc), col="tomato", lwd=2,bty = "n")

另一種方法，使用pROC進(jìn)行演示，還是那句話，不管是哪種方法，畫ROC曲線都是需要提供真實(shí)結(jié)果和預(yù)測(cè)概率！

library(pROC)

rocc <- roc(truth_test, prob) # 預(yù)測(cè)概率，真實(shí)結(jié)果
rocc # 看下結(jié)果
## 
## Call:
## roc.default(response = truth_test, predictor = prob)
## 
## Data: prob in 150 controls (truth_test pos) > 81 cases (truth_test neg).
## Area under the curve: 0.8477

# 畫圖
plot(rocc, 
     print.auc=TRUE, 
     auc.polygon=TRUE, 
     max.auc.polygon=TRUE, 
     auc.polygon.col="skyblue", 
     grid=c(0.1, 0.2), 
     grid.col=c("green", "red"), 
     print.thres=TRUE)

關(guān)于ROC曲線繪制的合集，共13篇文章，鏈接：ROC曲線繪制合集

順手再展示下PR曲線，也是用ROCR實(shí)現(xiàn)：

library(ROCR)

# ROCR畫ROC曲線就是2步，先prediction，再performance
pred <- prediction(prob,truth_test) # 預(yù)測(cè)概率，真實(shí)類別
perf <- performance(pred, "rec","prec") # ROC曲線的橫縱坐標(biāo)，不要寫錯(cuò)了
auc <- round(performance(pred, "auc")@y.values[[1]],digits = 4) # 提取AUC值
auc

# 畫圖
plot(perf,lwd=2,col="tomato")
# 添加圖例
legend("bottomright", legend=paste("AUC: ",auc), col="tomato", lwd=2,bty = "n")

是不是非常easy？

順手再畫個(gè)校準(zhǔn)曲線，公眾號(hào)后臺(tái)回復(fù)校準(zhǔn)曲線即可獲取合集，也是非常簡(jiǎn)單：

我這里給大家介紹最新的方法（其實(shí)之前也介紹過了），用probably這個(gè)包繪制：

library(probably)

cali_data <- data.frame(.pred_pos = prob, diabetes=truth_test)

cal_plot_breaks(cali_data,diabetes, .pred_pos,conf_level = 0.95)

但是目前這個(gè)版本（1.0.3）有個(gè)bug，第3個(gè)參數(shù)estimate，必須是.pred_xxx，其中的xxx必須是真實(shí)結(jié)果中的某一個(gè)類別，比如我這個(gè)數(shù)據(jù)diabetes中的類別就是pos和neg，那么這個(gè)名字就必須是.pred_pos或者.pred_neg，其他都會(huì)報(bào)錯(cuò)（下標(biāo)出界）??！

順手再畫個(gè)決策曲線，這個(gè)決策曲線是臨床預(yù)測(cè)模型中才有的內(nèi)容，其他內(nèi)容基本上都是機(jī)器學(xué)習(xí)的基礎(chǔ)知識(shí)。后臺(tái)回復(fù)決策曲線即可獲取合集：

source("datasets/dca.r")

# 把概率加到測(cè)試集中
dca_data <- pimadiabetes[-ind,]
dca_data$prob <- prob

# 結(jié)果變量變成0，1
dca_data$diabetes <- ifelse(dca_data$diabetes=="pos",1,0)

dc <- dca(data = dca_data, # 測(cè)試集
          outcome = "diabetes",
          predictors = "prob",
          probability = T
          )

前幾年stdca.r和dca.r這兩個(gè)腳本是可以在文中給出的網(wǎng)址中免費(fèi)下載的，但是從2022年底左右這個(gè)網(wǎng)站就不提供這兩段代碼的下載了。因?yàn)槲液茉缇拖螺d好了，所以我把這兩段代碼放在粉絲qq群文件里，大家有需要的加群下載即可。

超參數(shù)調(diào)優(yōu)

KNN算法只有一個(gè)超參數(shù)，就是近鄰的數(shù)量，所以KNN的超參數(shù)調(diào)優(yōu)其實(shí)就是如何確定最佳的K值，到底用幾個(gè)近鄰是能得到最好的結(jié)果呢？

現(xiàn)在有很多好用的工具可以實(shí)現(xiàn)調(diào)優(yōu)過程了，比如caret、tidymodels、mlr3等，但是這里我給大家演示下for循環(huán)的做法，因?yàn)樗挥?個(gè)超參數(shù)，很適合這種方法，還可以繪制學(xué)習(xí)曲線。

模型評(píng)價(jià)指標(biāo)選擇AUC。

aucs <- list()
for (i in 1:50) { # K的值選擇1~50
  f <- knn(train = train, # 訓(xùn)練集
         test = test, # 測(cè)試集
         cl = truth_train, # 訓(xùn)練集的真實(shí)類別
         k = i, # 使用的近鄰個(gè)數(shù)
         prob = TRUE # 需要計(jì)算概率
         )
  prob <- attr(f,"prob")
  prob <- ifelse(f == "pos", prob, 1-prob)
  pred <- prediction(prob,truth_test)
  perf <- performance(pred, "tpr","fpr")
  auc <- round(performance(pred, "auc")@y.values[[1]],digits = 4)
  aucs[[i]] <- auc
}
aucs <- do.call(rbind,aucs)[,1]
aucs # 50個(gè)AUC值，分別對(duì)應(yīng)50個(gè)K值
##  [1] 0.7068 0.7870 0.8180 0.8264 0.8505 0.8510 0.8481 0.8477 0.8586 0.8540
## [11] 0.8543 0.8580 0.8537 0.8513 0.8500 0.8451 0.8449 0.8453 0.8488 0.8467
## [21] 0.8450 0.8473 0.8464 0.8462 0.8474 0.8492 0.8483 0.8508 0.8519 0.8500
## [31] 0.8499 0.8504 0.8497 0.8491 0.8460 0.8455 0.8460 0.8427 0.8439 0.8417
## [41] 0.8422 0.8402 0.8393 0.8384 0.8393 0.8377 0.8377 0.8374 0.8356 0.8342

畫個(gè)圖看下不同的K值對(duì)應(yīng)的AUC變化的情況，看圖更加直觀：

plot_df <- data.frame(k=1:50,auc=aucs)
library(ggplot2)

ggplot(plot_df, aes(k,auc))+
  geom_line(linewidth=1)+
  geom_point(size=2)+
  geom_hline(yintercept = 0.85,linetype = 2)+
  geom_vline(xintercept = 9,linetype = 2,color="red")

結(jié)果顯示當(dāng)K=9的時(shí)候，AUC值是最大的，此時(shí)是0.8586。

這個(gè)圖其實(shí)是一個(gè)學(xué)習(xí)曲線圖，是一種經(jīng)典的進(jìn)行超參數(shù)調(diào)優(yōu)時(shí)使用的圖，我在介紹決策樹的超參數(shù)調(diào)優(yōu)時(shí)介紹過了，不知道大家有沒有印象？

所以此時(shí)你可以用K=9再重新跑一遍模型，作為你最終的結(jié)果。

final_f <- knn(train = train, # 訓(xùn)練集，只有預(yù)測(cè)變量，沒有結(jié)果變量
         test = test, # 測(cè)試集，沒有結(jié)果變量
         cl = truth_train, # 訓(xùn)練集的真實(shí)結(jié)果
         k = 9,            # 這里的K值選擇9哦?。?！
         prob = TRUE # 需要計(jì)算概率
         )

后續(xù)模型評(píng)價(jià)、畫ROC曲線就是一樣的代碼了，就不再重復(fù)了。

kknn包

數(shù)據(jù)劃分為訓(xùn)練集和測(cè)試集，劃分比例為7：3。

kknn包不需要把結(jié)果變量去掉。

# 劃分是隨機(jī)的，設(shè)置種子數(shù)可以讓結(jié)果復(fù)現(xiàn)
set.seed(123)
ind <- sample(1:nrow(pimadiabetes), size = 0.7*nrow(pimadiabetes))

# 訓(xùn)練集、測(cè)試集
train <- pimadiabetes[ind,]
test <- pimadiabetes[-ind, ]

# 把真實(shí)結(jié)果列單獨(dú)拿出來，后面用
truth_train <- pimadiabetes[ind,9]
truth_test <- pimadiabetes[-ind,9]

建立模型

在訓(xùn)練集擬合模型，支持R語言經(jīng)典的formula形式：

library(kknn)
fit <- kknn::kknn(diabetes ~ ., train, test,
                  scale = F, # w我們已經(jīng)對(duì)數(shù)據(jù)進(jìn)行過標(biāo)準(zhǔn)化了，這里就不用了
            )
# 直接summary可以查看預(yù)測(cè)類別和預(yù)測(cè)概率，太長(zhǎng)不展示
#summary(fit)

我們最關(guān)心的東西其實(shí)只有預(yù)測(cè)類別和預(yù)測(cè)概率而已，所以可以單獨(dú)查看它們：

# 預(yù)測(cè)類別
pred_class <- fit[["fitted.values"]]
head(pred_class)
## [1] neg neg pos neg pos pos
## Levels: pos neg

# 預(yù)測(cè)概率
pred_prob <- fit[["prob"]]
head(pred_prob)
##            pos       neg
## [1,] 0.3860350 0.6139650
## [2,] 0.3440688 0.6559312
## [3,] 1.0000000 0.0000000
## [4,] 0.2768308 0.7231692
## [5,] 0.6679614 0.3320386
## [6,] 0.5939196 0.4060804

而且這個(gè)包的結(jié)果給出了兩種類別的概率，不用再自己計(jì)算了。

模型評(píng)價(jià)

首先還是借助caret包查看混淆矩陣等各種信息：

caret::confusionMatrix(pred_class,truth_test)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction pos neg
##        pos 131  34
##        neg  19  47
##                                           
##                Accuracy : 0.7706          
##                  95% CI : (0.7109, 0.8232)
##     No Information Rate : 0.6494          
##     P-Value [Acc > NIR] : 4.558e-05       
##                                           
##                   Kappa : 0.4738          
##                                           
##  Mcnemar's Test P-Value : 0.05447         
##                                           
##             Sensitivity : 0.8733          
##             Specificity : 0.5802          
##          Pos Pred Value : 0.7939          
##          Neg Pred Value : 0.7121          
##              Prevalence : 0.6494          
##          Detection Rate : 0.5671          
##    Detection Prevalence : 0.7143          
##       Balanced Accuracy : 0.7268          
##                                           
##        'Positive' Class : pos             
##

然后是繪制ROC曲線，完全一樣的代碼：

library(ROCR)

pred <- prediction(pred_prob[,1],truth_test) # 預(yù)測(cè)概率，真實(shí)類別
perf <- performance(pred, "tpr","fpr")
auc <- round(performance(pred, "auc")@y.values[[1]],digits = 4)
auc
## [1] 0.8491

plot(perf,lwd=2,col="tomato")
abline(0,1,lty=2)
legend("bottomright", legend=paste("AUC: ",auc), col="tomato", lwd=2,bty = "n")

easy！一樣的用法，基本沒啥變化，所以pROC的畫法就不再重復(fù)了，大家想要學(xué)習(xí)的就自己寫一下即可。

超參數(shù)調(diào)優(yōu)

借助for循環(huán)也可以，這里再給大家演示下如何使用e1071包實(shí)現(xiàn)輕量化的超參數(shù)調(diào)優(yōu)。

library(e1071)

set.seed(123)
tune.knn(x=train[,-9], # 預(yù)測(cè)變量
         y=truth_train,# 結(jié)果變量
         k=1:50        # k的值
         )
## 
## Parameter tuning of 'knn.wrapper':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##   k
##  25
## 
## - best performance: 0.2198113