如何将pdp::partial映射到嵌套randomForest模型并解决变量找不到报错

错误原因 #

• 拟合分组模型时调用randomForest()指定的训练数据是全局全量的boston数据集，而非每个分组对应的子集（即map迭代的.x），不仅会导致两个分组的模型实际是同一个全量模型，也会引发后续偏依赖计算的训练数据匹配问题。
• 调用map2执行pdp::partial时没有正确引用迭代临时变量：map2中第一个迭代对象fit对应临时变量.x，第二个迭代对象分组data对应临时变量.y，直接写fit和data会指向全局环境下的同名变量，而非嵌套tibble每行对应的模型和分组数据。

修正后完整代码 #

library(tidyverse)
library(pdp)
library(randomForest)
data(boston)
glimpse(boston)

# 分组嵌套+拟合分组随机森林
boston %>%
  mutate(grp = ifelse(age < 80, "young", "old")) %>%
  nest(data = -grp) %>%
  # 修正点：使用迭代的分组子集作为训练数据
  mutate(fit = map(data, ~ randomForest(cmedv ~ ., data = .x, importance = TRUE))) %>%
  {. ->> GrpModels}

# 计算分组偏依赖
GrpModels %>%
  # 修正点：引用map2迭代的临时变量
  mutate(p = map2(fit, data, ~ pdp::partial(.x, train = .y, pred.var = c("lstat", "rm")))) %>%
  unnest(p) %>%
  {. ->> checkpdp}

分面3D偏依赖图绘制参考 #

可以搭配plotly实现分面3D展示：

library(plotly)

checkpdp %>%
  group_by(grp) %>%
  do(
    p = plot_ly(., x = ~lstat, y = ~rm, z = ~yhat, type = "mesh3d") %>%
      layout(scene = list(
        xaxis = list(title = "lstat"),
        yaxis = list(title = "rm"),
        zaxis = list(title = "cmedv偏依赖值")
      ), title = ~paste0("分组：", .$grp[1]))
  ) %>%
  subplot(nrows = 1, shareX = TRUE, shareY = TRUE)

内容的提问来源于stack exchange，提问作者Kevin