本文主要介紹如何使用R語言實現(xiàn)K最近鄰算法。
本文目錄：
算法簡介
準備數(shù)據(jù)
class包
建立模型
模型評價
超參數(shù)調(diào)優(yōu)
kknn包
建立模型
模型評價
超參數(shù)調(diào)優(yōu)
算法簡介
K最近鄰（K-Nearest-Neighbor，KNN）是一種非線性的分類算法，KNN處理分類問題的方法是：找K個距離待預測樣本最近的點，然后根據(jù)這幾個點的類別來確定新樣本的類別。
近朱者赤，近墨者黑。
下面以一個二分類問題為例說明KNN的思想。
下圖有兩個特征可以用來預測腫瘤是”良性”還是”惡性”。圖中的X表示我們要預測的新樣本。如果算法設定k=3，那么圓圈中包含的3個觀測就是樣本X的最近鄰。因為其中占多數(shù)比例的類別是”惡性”，所以樣本X被分類為”惡性”。
思想是不是很簡單？K的選擇對于KNN的預測結(jié)果是非常重要的。
KNN中另一個需要指出的重要問題是距離的計算方法，或者說是特征空間中數(shù)據(jù)點的臨近度的計算。默認的距離是歐氏距離，也就是從點A到點B的簡單直線距離。
Note
兩點間的距離強烈依賴于測量特征時使用的單位，所以必須對其進行標準化，而且要求數(shù)據(jù)不能有缺失值。
準備數(shù)據(jù)
演示數(shù)據(jù)為印第安人糖尿病數(shù)據(jù)集，這個數(shù)據(jù)一共有768行，9列，其中diabetes是結(jié)果變量，為二分類，其余列是預測變量。
該數(shù)據(jù)集的原始版本是有缺失值的，我這里使用的是插補過的版本，詳細過程請參考數(shù)據(jù)準備這一章。
rm(list = ls())
load(file = "../000機器學習/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 ...
各個變量的含義：
insulin：2小時血清胰島素（mu U/ml）pedigree：糖尿病譜系功能，是一種用于預測糖尿病發(fā)病風險的指標，該指標是基于家族史的糖尿病遺傳風險因素的計算得出的。它計算了患者的家族成員是否患有糖尿病以及他們與患者的親緣關系，從而得出一個綜合評分，用于預測患糖尿病的概率。
先對數(shù)據(jù)進行標準化：
# 對數(shù)值型變量進行標準化
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ù)劃分為訓練集和測試集，劃分比例為7：3。
但是R語言里class包在使用時需要把真實結(jié)果去掉，所以我們把真實結(jié)果去掉，只保留預測變量。
# 劃分是隨機的，設置種子數(shù)可以讓結(jié)果復現(xiàn)
set.seed(123)
ind <- sample(1:nrow(pimadiabetes), size = 0.7*nrow(pimadiabetes))
# 去掉真實結(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 ...
# 把真實結(jié)果列單獨拿出來，后面用
truth_train <- pimadiabetes[ind,9]
truth_test <- pimadiabetes[-ind,9]
建立模型
在訓練集建立模型，1行代碼搞定：
library(class)
f <- knn(train = train, # 訓練集，只有預測變量，沒有結(jié)果變量
test = test, # 測試集，沒有結(jié)果變量
cl = truth_train, # 訓練集的真實結(jié)果
k = 8, # 使用的近鄰個數(shù)
prob = TRUE # 需要計算概率
)
# 查看測試集的預測結(jié)果，只看前6個
head(f)
## [1] neg neg pos neg neg pos
## Levels: pos neg
# 查看測試集的預測概率，只看前6個
prob <- attr(f,"prob")
head(prob)
## [1] 0.500 0.750 1.000 0.750 0.625 0.625
此時得到的f這個結(jié)果是一個因子型的向量，而且是有名字和屬性的，大多數(shù)模型擬合結(jié)果的格式都是不一樣的，使用時需要注意！
還要注意這里的概率，并不是陽性結(jié)果的概率，而是預測結(jié)果的概率！比如第一個概率0.500是neg的概率，第二個概率0.750是neg的概率，第三個概率1.000是pos的概率！
所以你如果想要陽性結(jié)果(pos)的概率，需要自己計算一下：
prob <- ifelse(f == "pos", prob, 1-prob)
head(prob)
## [1] 0.500 0.250 1.000 0.250 0.375 0.625
模型評價
模型擬合好之后，下一步就是查看模型的各種指標，來看看這個模型在訓練集中的表現(xiàn)如何，比如混淆矩陣、AUC值、準確率等。
不管是什么類型的模型，如果我們想要評價它的模型表現(xiàn)，都是需要用到模型的預測結(jié)果和真實結(jié)果的。
對于回歸任務來說，預測結(jié)果也是數(shù)值型的，對于分類任務來說，模型的預測結(jié)果可以是某一種類別的概率，也可以是預測出的具體類別。大多數(shù)模型都是既支持計算類別概率又支持計算具體類別的，但是有些模型可能只支持一種類型。
首先查看混淆矩陣，我們借助caret包展示，這個包目前仍然是查看混線矩陣最全面的，，沒有之一，非常好用！
# 借助caret包，f是預測的類別，truth_test是真實的結(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é)果非常全面，最上面是混淆矩陣，然后給出了：
- 準確率（Accuracy）和準確率的可信區(qū)間（95% CI）
- 無信息率（No Information Rate）和P值（P-Value Acc>NIR]）、
- Mcnemar檢驗的P值（Mcnemar’s Test P-Value）
最后告訴你參考類別是pos。
當然這些值你也可以單獨計算：
caret::precision(f,truth_test) # 精準率
## [1] 0.803681
caret::recall(f,truth_test) # 召回率，靈敏度
## [1] 0.8733333
caret::F_meas(f,truth_test) # F1分數(shù)
## [1] 0.8370607
然后再畫個ROC曲線。首先使用ROCR進行演示，不管是什么包，都遵循前面說過的規(guī)律，畫ROC曲線是需要真實結(jié)果和預測概率的！
library(ROCR)
# ROCR畫ROC曲線就是2步，先prediction，再performance
pred <- prediction(prob,truth_test) # 預測概率，真實類別
perf <- performance(pred, "tpr","fpr") # ROC曲線的橫縱坐標，不要寫錯了
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) # 添加對角線
# 添加圖例
legend("bottomright", legend=paste("AUC: ",auc), col="tomato", lwd=2,bty = "n")
另一種方法，使用pROC進行演示，還是那句話，不管是哪種方法，畫ROC曲線都是需要提供真實結(jié)果和預測概率！
library(pROC)
rocc <- roc(truth_test, prob) # 預測概率，真實結(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)
關于ROC曲線繪制的合集，共13篇文章，鏈接：ROC曲線繪制合集
順手再展示下PR曲線，也是用ROCR實現(xiàn)：
library(ROCR)
# ROCR畫ROC曲線就是2步，先prediction，再performance
pred <- prediction(prob,truth_test) # 預測概率，真實類別
perf <- performance(pred, "rec","prec") # ROC曲線的橫縱坐標，不要寫錯了
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")
image-20240502163045511是不是非常easy？
順手再畫個校準曲線，公眾號后臺回復校準曲線即可獲取合集，也是非常簡單：
我這里給大家介紹最新的方法（其實之前也介紹過了），用probably這個包繪制：
library(probably)
cali_data <- data.frame(.pred_pos = prob, diabetes=truth_test)
cal_plot_breaks(cali_data,diabetes, .pred_pos,conf_level = 0.95)
但是目前這個版本（1.0.3）有個bug，第3個參數(shù)estimate，必須是.pred_xxx，其中的xxx必須是真實結(jié)果中的某一個類別，比如我這個數(shù)據(jù)diabetes中的類別就是pos和neg，那么這個名字就必須是.pred_pos或者.pred_neg，其他都會報錯（下標出界）?。?/p>
順手再畫個決策曲線，這個決策曲線是臨床預測模型中才有的內(nèi)容，其他內(nèi)容基本上都是機器學習的基礎知識。后臺回復決策曲線即可獲取合集：
source("datasets/dca.r")
# 把概率加到測試集中
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, # 測試集
outcome = "diabetes",
predictors = "prob",
probability = T
)
前幾年stdca.r和dca.r這兩個腳本是可以在文中給出的網(wǎng)址中免費下載的，但是從2022年底左右這個網(wǎng)站就不提供這兩段代碼的下載了。因為我很早就下載好了，所以我把這兩段代碼放在粉絲qq群文件里，大家有需要的加群下載即可。
超參數(shù)調(diào)優(yōu)
KNN算法只有一個超參數(shù)，就是近鄰的數(shù)量，所以KNN的超參數(shù)調(diào)優(yōu)其實就是如何確定最佳的K值，到底用幾個近鄰是能得到最好的結(jié)果呢？
現(xiàn)在有很多好用的工具可以實現(xiàn)調(diào)優(yōu)過程了，比如caret、tidymodels、mlr3等，但是這里我給大家演示下for循環(huán)的做法，因為它只有1個超參數(shù)，很適合這種方法，還可以繪制學習曲線。
模型評價指標選擇AUC。
aucs <- list()
for (i in 1:50) { # K的值選擇1~50
f <- knn(train = train, # 訓練集
test = test, # 測試集
cl = truth_train, # 訓練集的真實類別
k = i, # 使用的近鄰個數(shù)
prob = TRUE # 需要計算概率
)
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個AUC值，分別對應50個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
畫個圖看下不同的K值對應的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é)果顯示當K=9的時候，AUC值是最大的，此時是0.8586。
這個圖其實是一個學習曲線圖，是一種經(jīng)典的進行超參數(shù)調(diào)優(yōu)時使用的圖，我在介紹決策樹的超參數(shù)調(diào)優(yōu)時介紹過了，不知道大家有沒有印象？
所以此時你可以用K=9再重新跑一遍模型，作為你最終的結(jié)果。
final_f <- knn(train = train, # 訓練集，只有預測變量，沒有結(jié)果變量
test = test, # 測試集，沒有結(jié)果變量
cl = truth_train, # 訓練集的真實結(jié)果
k = 9, # 這里的K值選擇9哦！??！
prob = TRUE # 需要計算概率
)
后續(xù)模型評價、畫ROC曲線就是一樣的代碼了，就不再重復了。
kknn包
數(shù)據(jù)劃分為訓練集和測試集，劃分比例為7：3。
kknn包不需要把結(jié)果變量去掉。
# 劃分是隨機的，設置種子數(shù)可以讓結(jié)果復現(xiàn)
set.seed(123)
ind <- sample(1:nrow(pimadiabetes), size = 0.7*nrow(pimadiabetes))
# 訓練集、測試集
train <- pimadiabetes[ind,]
test <- pimadiabetes[-ind, ]
# 把真實結(jié)果列單獨拿出來，后面用
truth_train <- pimadiabetes[ind,9]
truth_test <- pimadiabetes[-ind,9]
建立模型
在訓練集擬合模型，支持R語言經(jīng)典的formula形式：
library(kknn)
fit <- kknn::kknn(diabetes ~ ., train, test,
scale = F, # w我們已經(jīng)對數(shù)據(jù)進行過標準化了，這里就不用了
)
# 直接summary可以查看預測類別和預測概率，太長不展示
#summary(fit)
我們最關心的東西其實只有預測類別和預測概率而已，所以可以單獨查看它們：
# 預測類別
pred_class <- fit[["fitted.values"]]
head(pred_class)
## [1] neg neg pos neg pos pos
## Levels: pos neg
# 預測概率
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
而且這個包的結(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) # 預測概率，真實類別
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的畫法就不再重復了，大家想要學習的就自己寫一下即可。
超參數(shù)調(diào)優(yōu)
借助for循環(huán)也可以，這里再給大家演示下如何使用e1071包實現(xiàn)輕量化的超參數(shù)調(diào)優(yōu)。
library(e1071)
set.seed(123)
tune.knn(x=train[,-9], # 預測變量
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
1行代碼出結(jié)果，默認是使用10折交叉驗證，比我們的手動for循環(huán)更加穩(wěn)健，結(jié)果最佳的k值是25。使用的評價指標不同，具體計算的步驟也不一樣，得出的結(jié)果不一樣是很正常的。
