R語言機器學習：caret包使用及其黑箱模型解釋（離散變數預測）

作者：黃天元，復旦大學博士在讀，熱愛數據科學與開源工具（R），致力於利用數據科學迅速積累行業經驗優勢和學術知識發現。知乎專欄：R語言數據挖掘 郵箱：[email protected].歡迎合作交流。

上次在

HopeR：R語言機器學習：caret包使用及其黑箱模型解釋（連續變數預測)?

zhuanlan.zhihu.com

中介紹瞭如何使用caret包建模並採用DALEX包進行模型的解釋。當時是針對連續型變數進行探索，這次我們針對響應變數為離散變數（分類變數）的模型進行黑箱解釋。

1 包的載入與數據導入

安裝4個包。

library(pacman) p_load(DALEX,caret,tidyverse,breakDown)

觀察我們要使用的目標數據：

library(breakDown) data(wine) ? wine %>% as_tibble

# A tibble: 4,898 x 12
fixed.acidity volatile.acidity citric.acid residual.sugar chlorides free.sulfur.dio~ total.sulfur.di~ density
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 7 0.27 0.36 20.7 0.045 45 170 1.00
2 6.3 0.3 0.34 1.6 0.049 14 132 0.994
3 8.1 0.28 0.4 6.9 0.05 30 97 0.995
4 7.2 0.23 0.32 8.5 0.058 47 186 0.996
5 7.2 0.23 0.32 8.5 0.058 47 186 0.996
6 8.1 0.28 0.4 6.9 0.05 30 97 0.995
7 6.2 0.32 0.16 7 0.045 30 136 0.995
8 7 0.27 0.36 20.7 0.045 45 170 1.00
9 6.3 0.3 0.34 1.6 0.049 14 132 0.994
10 8.1 0.22 0.43 1.5 0.044 28 129 0.994
# ... with 4,888 more rows, and 4 more variables: pH <dbl>, sulphates <dbl>, alcohol <dbl>, quality <fct>

下面，我們對酒的質量進行因子花，如果質量數值大於5，則標誌為1，否則標誌為0.

wine$quality <- factor(ifelse(wine$quality > 5, 1, 0))

2 使用caret包迅速建模

本次建模過程中，我們會手動劃分訓練集和測試集。其中，訓練集的比例為60%。

trainIndex <- createDataPartition(wine$quality, p = 0.6, list = FALSE, times = 1) wineTrain <- wine[ trainIndex,] wineTest <- wine[-trainIndex,]

隨後，我們開始建立模型。我們分別創建隨機森林、邏輯回歸和支持向量機模型。除了酒的質量以外，其他變數都設為解釋變數。

classif_rf <- train(quality~., data = wineTrain, method="rf", ntree = 100, tuneLength = 1) ? classif_glm <- train(quality~., data = wineTrain, method="glm", family="binomial") ? classif_svm <- train(quality~., data = wineTrain, method="svmRadial", prob.model = TRUE, tuneLength = 1)

3 對模型進行解釋

這裡直接利用DALEX包的explain函數對三個模型進行解釋性分析。需要注意的是，做這個分析需要包含4個信息：1.模型信息；2.標籤信息（如果沒有，會自動從模型抽取）；3.驗證數據集；4.驗證數據集中哪個是響應變數。除此以外，在這個問題中特別需要注意的是：這些建模函數默認響應變數是因子變數的時候，才會進行分類。但是就得到的結果而言，必須是似然函數得到的概率值，纔能夠對其進行更加精細的解釋。所以，必須把預測函數進行修飾，讓它能夠返回極大似然估計，而不是0和1。此外，對於驗證數據集而言，必須要把它的響應變數轉化為數值，纔能夠與預測值進行比較。代碼如下：

#定義求得似然估計的函數 p_fun <- function(object, newdata){predict(object, newdata=newdata, type="prob")[,2]} ? #把測試集的響應變數轉化為數值 yTest <- as.numeric(as.character(wineTest$quality)) ? explainer_classif_rf <- DALEX::explain(classif_rf, label = "rf", data = wineTest, y = yTest, predict_function = p_fun) ? explainer_classif_glm <- DALEX::explain(classif_glm, label = "glm", data = wineTest, y = yTest, predict_function = p_fun) ? explainer_classif_svm <- DALEX::explain(classif_svm, label = "svm", data = wineTest, y = yTest, predict_function = p_fun)

建模可能很久，但是解釋性驗證是非常快的，直接是黑箱的映射關係。

4 模型表現

後面的步驟，基本與之前連續變數的相似，不再進行贅述，感興趣請觀看https://zhuanlan.zhihu.com/p/59329696。

mp_classif_rf <- model_performance(explainer_classif_rf) mp_classif_glm <- model_performance(explainer_classif_glm) mp_classif_svm <- model_performance(explainer_classif_svm) ? plot(mp_classif_rf, mp_classif_glm, mp_classif_svm)

plot(mp_classif_rf, mp_classif_glm, mp_classif_svm, geom = "boxplot")

就這個結果看來，隨機森林的效果最好。

5 變數重要性分析

需要看每個模型中，不同變數對於模型預測的相對重要性，可以用如下方法。

vi_classif_rf <- variable_importance(explainer_classif_rf, loss_function = loss_root_mean_square) vi_classif_glm <- variable_importance(explainer_classif_glm, loss_function = loss_root_mean_square) vi_classif_svm <- variable_importance(explainer_classif_svm, loss_function = loss_root_mean_square) ? plot(vi_classif_rf, vi_classif_glm, vi_classif_svm)

損失函數使用的是RMSE，這裡解釋為：如果模型少了這個變數，將會給響應變數的預測值帶來多大影響（我們已經把因子編程數值，因此最小化的還是RMSE）？自信觀察可以發現，rf和glm兩個模型大相徑庭，也就是說這些模型對特徵的利用方法不同，導致它們對特徵看重的程度不同。我自己把它稱之為模型的特徵利用率。

6 變數解析

下面，我們分別做PDP圖和ALE圖。我們先對pH值這個變數進行分析，PDP圖如下：

pdp_classif_rf <- variable_response(explainer_classif_rf, variable = "pH", type = "pdp") pdp_classif_glm <- variable_response(explainer_classif_glm, variable = "pH", type = "pdp") pdp_classif_svm <- variable_response(explainer_classif_svm, variable = "pH", type = "pdp") ? plot(pdp_classif_rf, pdp_classif_glm, pdp_classif_svm)

我們看到，在前面認為殘差較小的隨機森林模型中，隨著pH值的變化，整個響應變數的變化是非線性的。SVM也稍微能夠捕捉到這個信息，但是邏輯回歸則完全不行。下面來看ALE圖：

ale_classif_rf <- variable_response(explainer_classif_rf, variable = "alcohol", type = "ale") ale_classif_glm <- variable_response(explainer_classif_glm, variable = "alcohol", type = "ale") ale_classif_svm <- variable_response(explainer_classif_svm, variable = "alcohol", type = "ale") ? plot(ale_classif_rf, ale_classif_glm, ale_classif_svm)