R语言tidymodels那些事（十五） 模型的可解释性

缘由

艺多不压身

碎碎念

我们继续飞行在云层和雾气之间，飞行在一个无光无影的空洞世界中。 ——《风沙星辰》

Ref

18 Explaining Models and Predictions | Tidy Modeling with R

解释

本文使用 Ames 房价数据对比了两个模型（线性回归 vs 随机森林），并用局部分解（breakdown）、SHAP、置换重要性（permutation）与部分依赖图（PDP）来解释模型。局部解释能告诉你“这套房为什么贵/便宜”，而全局解释告诉你“哪些变量对整体预测最关键”**。

---

1. 加载包与拆分数据

library(DALEXtra)    # 模型解释（predict_parts, model_parts, model_profile 等）
library(embed)       #
library(tidymodels)  # 建模生态（parsnip, recipes 等）
library(modeldata)   # ames 数据集
library(forcats)     # 因子重排（绘图）
library(ggplot2)     # 绘图
library(dplyr)       # 数据处理
library(stringr)     # 字符串操作
 
tidymodels_prefer()  # 确保 tidymodels 的函数优先

data(ames)
ames <- mutate(ames, Sale_Price = log10(Sale_Price))  # 目标取对数，便于建模与解释
 
set.seed(502)
ames_split <- initial_split(ames, prop = 0.80, strata = Sale_Price)
ames_train <- training(ames_split)
ames_test  <- testing(ames_split)

2.拟合模型：线性回归与随机森林

 
lm_model <- linear_reg() %>% set_engine("lm")
lm_fit <- fit(
  lm_model,
  Sale_Price ~ Neighborhood + Gr_Liv_Area + Year_Built + Bldg_Type + Latitude + Longitude,
  data = ames_train
)
 
rf_model <- rand_forest(trees = 1000, min_n = 5) %>%
  set_engine("ranger", verbose = TRUE) %>%
  set_mode("regression")
 
rf_fit <- fit(
  rf_model,
  Sale_Price ~ Neighborhood + Gr_Liv_Area + Year_Built + Bldg_Type + Latitude + Longitude,
  data = ames_train
)

3.构建结果的解释模型 explainer

#选择需要解释带特征 特征需要和建模的选择的特征一致
vip_features <- c("Neighborhood", "Gr_Liv_Area", "Year_Built", 
                  "Bldg_Type", "Latitude", "Longitude")
 
vip_train <- ames_train %>% select(all_of(vip_features))
 
 
explainer_lm <- explain_tidymodels(
  lm_fit,
  data = vip_train,
  y = ames_train$Sale_Price,
  label = "lm + interactions",
  verbose = FALSE
)
 
explainer_rf <- explain_tidymodels(
  rf_fit,
  data = vip_train,
  y = ames_train$Sale_Price,
  label = "random forest",
  verbose = FALSE
)

4.局部解释（单个观测的分解 / SHAP）

单个观测

把单个观测（如 duplex）的预测拆成「截距 + 各特征贡献 = 预测值」。告诉你“对于这套房子，哪些特征把预测往上拉，哪些往下拉”，直观。各 feature 的 contribution（正为推高预测，负为拉低）
prediction（最后的预测值 = 截距 + 各贡献之和）

# 选择单个观测（示例）：duplex
duplex <- vip_train[120,]
print(duplex)  # 查看该房子的特征
 
# 1) 单次分解（additive breakdown）
lm_breakdown <- predict_parts(explainer = explainer_lm, new_observation = duplex)
rf_breakdown <- predict_parts(explainer = explainer_rf, new_observation = duplex)
 
# 查看结果（表格或打印）
print(lm_breakdown)
print(rf_breakdown)
 
 

SHAP

通过近似 Shapley 值为单个观测计算每个特征的贡献分布（重复抽样得到贡献的分布/不确定性）。比单次 breakdown 更稳健（考虑不同特征加入的次序），并给出贡献的分布信息（箱线/波动），便于判断贡献稳定性。 要点：B 控制采样次数（越大越稳定，但计算越慢）。

set.seed(1801)
shap_duplex <- predict_parts(explainer = explainer_rf, new_observation = duplex, type = "shap", B = 20)
 
# 绘制 SHAP 的可视化（DALEX 自带 plot 方法）
plot(shap_duplex)
 
# 或者自定义箱线图来展示 SHAP 分布（更美观）
shap_duplex %>%
  group_by(variable) %>%
  mutate(mean_val = mean(contribution)) %>%
  ungroup() %>%
  mutate(variable = fct_reorder(variable, abs(mean_val))) %>%
  ggplot(aes(contribution, variable, fill = mean_val > 0)) +
  geom_col(data = ~distinct(., variable, mean_val), aes(mean_val, variable), alpha = 0.5) +
  geom_boxplot(width = 0.5) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(
    title = "SHAP contributions for a single observation (box = distribution over permutations)",
    x = "Contribution to prediction (log10 price)",
    y = NULL
  )

SHAP 可以用来比较不同样本的贡献分布差异（比如小房 vs 大房）。

# 再选一个样本（大房子）做对比
big_house <- vip_train[1269,]
print(big_house)
 
set.seed(1802)
shap_house <- predict_parts(explainer = explainer_rf, new_observation = big_house, type = "shap", B = 20)
plot(shap_house)

5.全局解释（Permutation importance）

把某一特征的值在数据集中随机置换，观察模型的性能下降（这里用 RMSE）有多少，用下降幅度衡量该特征对模型整体性能的重要性。给出“全局层面”的特征重要性：如果打乱这个变量，模型预测会变多糟（越糟表示越重要）。

注意：对强相关特征（例如经纬度分别与地段信息相关）置换重要性可能会低估某些变量，因为信息可被其它相关变量代替。permutation 是衡量“对模型性能的敏感度”，不是衡量“真实因果性”。

set.seed(1803)
vip_lm <- model_parts(explainer_lm, loss_function = loss_root_mean_square)
 
set.seed(1804)
vip_rf <- model_parts(explainer_rf, loss_function = loss_root_mean_square)
 
# 自定义绘图函数，比较两个模型的置换重要性分布
library(forcats)
ggplot_imp <- function(...) {
  obj <- list(...)
  metric_name <- attr(obj[[1]], "loss_name")
  metric_lab <- paste(metric_name, "after permutations\n(higher indicates more important)")
  
  full_vip <- bind_rows(obj) %>% filter(variable != "_baseline_")
  
  perm_vals <- full_vip %>% 
    filter(variable == "_full_model_") %>% 
    group_by(label) %>% 
    summarise(dropout_loss = mean(dropout_loss))
  
  p <- full_vip %>%
    filter(variable != "_full_model_") %>% 
    mutate(variable = fct_reorder(variable, dropout_loss)) %>%
    ggplot(aes(dropout_loss, variable))
  
  if (length(obj) > 1) {
    p <- p + 
      facet_wrap(vars(label)) +
      geom_vline(data = perm_vals, aes(xintercept = dropout_loss, color = label),
                 linewidth = 1.4, lty = 2, alpha = 0.7) +
      geom_boxplot(aes(color = label, fill = label), alpha = 0.2)
  } else {
    p <- p + 
      geom_vline(data = perm_vals, aes(xintercept = dropout_loss),
                 linewidth = 1.4, lty = 2, alpha = 0.7) +
      geom_boxplot(fill = "#91CBD765", alpha = 0.4)
  }
  
  p + theme_minimal() + theme(legend.position = "none") +
    labs(x = metric_lab, y = NULL, fill = NULL, color = NULL)
}
 
# 绘图
ggplot_imp(vip_lm, vip_rf)

RMSE越大，说明特征越重要

6. 部分依赖图（PDP / ICE），包括分组 PDP（查看交互）

展示当我们把某个变量（例如 Year_Built 或 Gr_Liv_Area）从小到大改变时，模型预测的平均响应（并可按组展示交互差异）。可视化变量的平均效应和非线性趋势；当按 Bldg_Type 分组时，还能看到交互效应（不同建筑类型的面积效应不同）。
注意：PDP 展示的是平均效果，可能掩盖不同个体之间的差异（这时要看 ICE/分组曲线）。

# Year_Built 的 PDP（平均效应）
set.seed(1805)
pdp_age <- model_profile(explainer_rf, N = 500, variables = "Year_Built")
 
# 绘制 PDP 的函数（显示 ICE 灰线 + 平均曲线）
ggplot_pdp <- function(obj, x) {
  p <- as_tibble(obj$agr_profiles) %>%
    mutate(`_label_` = str_remove(`_label_`, "^[^_]*_")) %>%
    ggplot(aes(`_x_`, `_yhat_`)) +
    geom_line(data = as_tibble(obj$cp_profiles),
              aes(x = !!rlang::sym(x), group = `_ids_`),
              linewidth = 0.5, alpha = 0.05, color = "gray50")
  
  num_colors <- n_distinct(obj$agr_profiles$`_label_`)
  
  if (num_colors > 1) {
    p <- p + geom_line(aes(color = `_label_`), linewidth = 1.2, alpha = 0.8)
  } else {
    p <- p + geom_line(color = "midnightblue", linewidth = 1.2, alpha = 0.8)
  }
  
  p + theme_minimal()
}
 
# 绘图并加注
ggplot_pdp(pdp_age, "Year_Built") +
  labs(x = "Year built", y = "Sale Price (log10)", color = NULL, 
       title = "Partial dependence of Year_Built (random forest)")
 
# Gr_Liv_Area 的 PDP，按 Bldg_Type 分组（查看交互）
set.seed(1806)
pdp_liv <- model_profile(explainer_rf, N = 1000, variables = "Gr_Liv_Area", groups = "Bldg_Type")
 
# 用原始 agr_profiles + cp_profiles 画分组 PDP（并用对数 x 轴）
as_tibble(pdp_liv$agr_profiles) %>%
  mutate(Bldg_Type = str_remove(`_label_`, "random forest_")) %>%
  ggplot(aes(`_x_`, `_yhat_`, color = Bldg_Type)) +
  geom_line(data = as_tibble(pdp_liv$cp_profiles),
            aes(x = Gr_Liv_Area, group = `_ids_`),
            linewidth = 0.5, alpha = 0.1, color = "gray50") +
  geom_line(linewidth = 1.2, alpha = 0.8, show.legend = FALSE) +
  scale_x_log10() +
  facet_wrap(~Bldg_Type) +
  theme_minimal() +
  labs(x = "Gross living area (log scale)", y = "Sale Price (log10)", 
       title = "Partial dependence of area by building type")

常见问题解答

1. 为什么要“局部→全局→再局部聚合”？三者各自的独特价值

• **局部（单样本）解释：回答“这套具体房子为什么被预测成这个价？”

  • 优点：直观、可操作（便于单条反馈或案例研究）。
  • 局限：只是个别样本的解释，可能不代表总体。

• **全局（Permutation / PDP）：回答“总体上哪些特征最重要？这些特征与目标之间的平均关系是什么？”

  • 优点：把模型放大到整体看，帮助把握宏观规律（例如面积和地理位置总体重要）。
  • 局限：平均化会掩盖个体差异与交互结构。

• 局部聚合（把很多个体的 SHAP 取平均/绝对值等）：把大量局部解释汇总成稳定的全局视角（例如用平均绝对 SHAP 值做全局重要性），同时保留个体层面的分布信息（不只单个数值）。

  • 优点：结合了局部的真实分布（不再是“单一平均”），能揭示“总体重要性同时是否稳定/有分歧”。