随机森林中的变量重要性

随机森林中的变量重要性研究' (Study on Variable Importance in Random Forest)

传统方法与新发展

Distributional(Random) Forest的特点。本文中：能够产生变量重要性。来源：作者。

随机森林以及其推广（尤其是广义随机森林（GRF）和分布随机森林（DRF））是强大且易于使用的机器学习方法，任何数据科学家的工具箱都不能缺少。它们不仅在大范围的数据集上表现出稳健的性能，而且无需调整就可以轻松处理缺失值，甚至可以提供置信区间。在本文中，我们着重讨论它们能够提供的另一个特性：特征重要性的概念。具体而言，我们关注以下内容：

传统随机森林（RF），用于预测给定p个预测变量X的变量Y的条件期望。
分布随机森林，用于预测给定p个预测变量X的d维变量Y的全部条件分布。

不幸的是，像许多现代机器学习方法一样，这两种方法都缺乏解释性。也就是说，涉及了很多操作，似乎不可能确定预测变量和Y之间的实际功能关系。解决这个问题的一种常见方法是定义变量重要性度量（VIMP），至少帮助确定哪些预测变量是重要的。通常，这有两个不同的目标：

(1) 找到少量具有最大准确性的变量。

(2) 检测和排名所有具有影响力的变量，以便进一步研究。

(1) 和 (2) 之间的差异在于，只要X的元素之间存在依赖关系（几乎总是存在），就会产生影响。例如，如果两个变量彼此高度相关，并且与Y高度相关，则可以在不损害目标 (1) 准确性的情况下去除其中一个输入，因为两个变量传达了相同的信息。然而，对于目标 (2)，两个变量都应该包括在内，因为这两个变量在实践中可能具有不同的含义，对于领域专家来说。

今天我们关注的是目标 (1)，并试图找到一个较小数量的预测变量，这些变量显示出更或多或少相同的预测准确性。例如，在下面的工资示例中，我们能够将预测变量的数量从79减少到约20个，而准确性只有轻微减少。这些最重要的预测变量包含了一些公认会影响工资的变量，例如年龄和教育。关于目标 (2)，VoAGI上也有许多优秀的文章，使用了Shapley值，比如这篇文章或者这篇文章。关于用随机森林高效计算Shapley值的最新而令人兴奋的学术文献也有。但这是另一篇文章的内容。

今天我们要看的这两个度量实际上是更通用的变量重要性度量，可以用于任何方法，基于我们将在下面介绍的扔掉和重新学习的原则。在这里，我们仅关注基于树的方法。此外，我们不会详细解释这些方法，而是尝试关注它们的应用以及为什么较新版本比传统版本更可取。

随机森林的变量重要性度量概览。平均减少不纯度（MDI）和平均减少准确性（MDA）均由Breiman提出。但由于它们的经验性质，存在一些问题，这些问题最近由Sobol-MDA得到解决。来源：作者

起源

随机森林（RF）的变量重要性度量实际上就像RF本身一样古老。第一个准确度度量平均减少准确率（MDA）是由Breiman在他的开创性随机森林论文[1]中提出的。这个想法很简单：对于每个维度j=1,…,p，比较完整预测的准确率与将X_j随机置换后的预测准确率。这样做的目的是破坏X_j与Y之间的关系，并将X_j不帮助预测Y与潜在有用性的准确度进行比较。

在R和Python中，有不同版本的MDA实现：

不幸的是，以这种方式置换变量X_j不仅会破坏它与Y的关系，也会破坏它与其他变量X的关系。如果X_j与所有其他变量都是独立的，这不是一个问题，但一旦存在依赖关系，就会成为一个问题。因此，[3]能够证明一旦X中存在依赖性，MDA就会收敛到一个荒谬的结果。特别是，MDA可能会赋予与预测Y不相关但与另一个变量X_l高度相关的变量X_j高重要性（如下面的示例所示）。同时，它可能无法检测到实际相关的变量，正如[3，第2.1节]中的一系列论文所示。直观地说，我们想要测量的是如果不包括X_j，则模型的性能，而我们测量的是具有置换过的X_j变量的模型的性能。

第二个传统的准确度度量是平均减少纯度（MDI），它将在森林中的所有树分裂在给定协变量上的减少纯度加权求和，再求平均。不幸的是，MDI从一开始就没有明确定义（不清楚它应该衡量什么），一些论文强调了这种方法的实际问题（例如[5]）。因此，我们不会详细介绍MDI，因为MDA通常是首选。

现代发展 I：Sobol-MDA

很长一段时间里，我一直以为这些有些非正式的度量是我们能做到的最好的了。最近，一篇论文改变了这种看法。在这篇论文中，作者从理论上证明了上述流行的度量实际上相当有缺陷，没有衡量我们想衡量的东西。那么第一个问题可能是：我们到底想衡量什么？一个潜在的答案是：Sobol指数（最初在计算机科学文献中提出）：

让我们来解释一下。首先，tau(X)=E[ Y | X]是我们想要估计的条件期望函数。这是一个随机变量，因为它是随机变量X的函数。现在X^{(-j)}是p-1维的向量，去掉了协变量j。因此，ST^{(j)}是如果去掉第j个输出变量后输出解释方差的减少量。

上面是更传统的度量方式。然而，对我来说，以下方式更直观：

这里d是两个随机向量之间的距离，对于上面的ST^{(j)}，这个距离就是通常的欧几里得距离。因此，ST^{(j)}的上部简单地衡量了我们想要的（tau(X)）与没有变量j的情况之间的平均平方距离。后者是

问题变成如何高效地估计这个问题。事实证明，直观的放弃和重新学习原则就足够了：只需使用RF估计tau(X)，然后删除X_j并重新拟合RF以获得tau(X^{(-j)})的估计，就可以获得一致的估计：

其中tau_n(X_i)是使用所有p个预测变量对测试点X_i进行的RF估计，类似地，tau_n(X_i^{(-j)})是只使用p-1个预测变量重新拟合的森林。

然而，这意味着森林需要重新拟合p次，当p很大时效率不高！因此，[3]中的作者开发了他们称之为Sobol-MDA的方法。与每次重新拟合森林不同，森林只需拟合一次。然后，测试点可以通过同一森林进行投影，从而得到预测结果“投影”形成(1)中的度量。也就是说，忽略X_j上的分割（记住目标是获得不包含X_j的估计结果）。作者能够证明使用这种投影方法计算上述(1)也会得到一致的估计！这确实是一个很好的想法，并使算法即使在高维情况下也是适用的。

投影方法的示意图。左边是RF对二维空间的划分。右边的投影方法忽略X^(2)上的分割，在进行预测时将其删除。可以看到根据这个原则，点X在右边被投影到X^{(-j)}上。来源：[3]中的图1

这种方法在R中实现了soboldMDA包，该包基于非常快速的ranger包。

现代发展II：基于MMD的灵敏度指数

从使用距离d的公式来看，一个自然的问题是是否可以使用不同的距离来获得更复杂问题的变量重要性度量。最近的一个例子是使用MMD距离作为d：

MMD距离是一个很好的工具，可以使用核函数k（如高斯核）相对容易地构建分布之间的距离：

暂时不详细说明，重要的一点是I^{(j)}考虑了比条件期望更一般的目标。它认识到如果变量X_j以任何方式影响Y的分布，那么它就是重要的。可能是X_j仅仅改变方差或分位数，而不影响Y的条件均值（见下面的示例）。在这种情况下，Sobol-MDA将不会将X_j识别为重要的，但是MMD方法会。这并不一定使它更好，它只是一种不同的工具：如果想预测条件期望，ST^{(j)}是正确的度量。然而，如果想预测分布的其他方面，特别是分位数，I^{(j)}更适用。同样，I^{(j)}可以使用放弃和重新学习的原则（对于每个$j$，使用移除变量$j$时重新拟合DRF），或者可以使用与Sobol-MDA相同的投影方法。本文末尾附有基于放弃和重新学习的实现。我们在这里将这个方法称为MMD-MDA。

模拟数据

我们接下来通过一个简单的模拟例子来演示这两种现代测量方法：我们首先从Gitlab下载并安装 Sobol-MDA 软件包，然后加载所有必要的软件包以供本例使用：

library(kernlab)library(drf)library(Matrix)library(DescTools)library(mice)library(sobolMDA)source("compute_drf_vimp.R") ##该文件内容可以在下面找到source("evaluation.R") ##该文件内容可以在下面找到

然后我们从这个简单例子中进行模拟：我们取独立于 (-1,1) 之间的 X_1、X_2、X_4、…、X_10 的均匀分布，并通过取 X_3=X_1 + 均匀误差来创建 X_1 和 X_3 之间的相关性。然后我们按照以下方式模拟 Y：

##模拟既有均值移动又有标准差移动的数据# 从 X_1 <- runif(n,-1,1) 运算的结果中模拟数据x2 <- runif(n,-1,1)X0 <- matrix(runif(7*n,-1,1), nrow=n, ncol=7)x3 <- x1+ runif(n,-1,1)X <- cbind(x1,x2, x3, X0)# 模拟相关变量 YY <- as.matrix(rnorm(n,mean = 0.8*(x1 > 0), sd = 1 + 1*(x2 > 0)))colnames(X)<-paste0("X", 1:10)head(cbind(Y,X))

然后我们使用 Sobol-MDA 方法对给定的 X 估计条件期望值进行分析：

##用于条件期望值估计的变量重要性XY <- as.data.frame(cbind(Xfull, Y))colnames(XY) <- c(paste('X', 1:(ncol(XY)-1), sep=''), 'Y')num.trees <- 500forest <- sobolMDA::ranger(Y ~., data = XY, num.trees = num.trees, importance = 'sobolMDA')sobolMDA <- forest$variable.importancenames(sobolMDA) <- colnames(X)sort(sobolMDA, decreasing = T)          X1           X8           X7           X6           X5           X9  0.062220958  0.021946135  0.016818860  0.016777223 -0.001290326 -0.001540919           X3          X10           X4           X2 -0.001578540 -0.007400854 -0.008299478 -0.020334150

可以看到，它正确地确定了 X_1 是最重要的变量，而其他变量的重要性排名相同(不)重要。这是因为给定 X_1，Y 的条件期望值才会发生变化。尽管 X_1 和 X_3 之间存在相关性，该测量方法仍然能够实现这一点。因此，在这个例子中，我们成功地实现了上述目标 (1)。另一方面，我们也可以观察传统的 MDA 方法：

forest <- sobolMDA::ranger(Y ~., data = XY, num.trees = num.trees, importance = 'permutation')MDA <- forest$variable.importancenames(MDA) <- colnames(X)sort(MDA, decreasing = T)          X1           X3           X6           X7           X8           X2  0.464516976  0.118147061  0.063969310  0.032741521  0.029004312 -0.004494380           X4           X9          X10           X5 -0.009977733 -0.011030996 -0.014281844 -0.018062544

在这种情况下，虽然它正确地确定了 X_1 是最重要的变量，但它也将 X_3 置于第二位，并且其值似乎比其他变量高得多。尽管如此，实际上 X_3 和 X_2、X_4、…、X_10 一样不重要！

但是，如果我们对Y的分布更感兴趣，例如用于估计分位数，那该怎么办呢？在这种情况下，我们需要一个能够识别X_2对Y的条件方差影响的度量。在这种情况下，MMD变量重要性度量就发挥作用了：

MMDVimp <- compute_drf_vimp(X=X,Y=Y)sort(MMDVimp, decreasing = T)         X2          X1         X10          X6          X8          X3 0.683315006 0.318517259 0.014066410 0.009904518 0.006859128 0.005529749          X7          X9          X4          X5 0.003476256 0.003290550 0.002417677 0.002036174

再次，这个度量能够正确识别出重要的因素：X_1和X_2是最重要的两个变量。同样，在X_1和X_3之间存在依赖关系的情况下，它也能够做到这一点。有趣的是，它也比X_1的期望变化更重视X_2的方差变化。

真实数据

最后，我呈现一个真实数据应用程序，以演示变量重要性度量。请注意，使用DRF，我们甚至可以查看多变量Y，但为了简化问题，我们选择了单变量设置，并考虑了来自美国普查局的2018年美国社区调查的美国工资数据。在第一篇DRF论文中，我们从2018年美国社区调查的大约100万全职员工中获取了数据，并从中提取了薪资信息和可能与薪资相关的所有协变量。这样丰富的数据对于尝试像DRF这样的方法非常理想（实际上，我们只会使用其中的一个小子集进行分析）。加载的数据可以在这里找到。

# 加载数据（https://github.com/lorismichel/drf/blob/master/applications/wage_data/data/datasets/wage_benchmark.Rdata）load("wage_benchmark.Rdata")##定义训练数据n<-1000Xtrain<-X[1:n,] Ytrain<-Y[1:n,]Xtrain<-cbind(Xtrain,Ytrain[,"male"])colnames(Xtrain)[ncol(Xtrain)]<-"male"Ytrain<-Ytrain[,1, drop=F]##定义测试数据ntest<-2000Xtest<-X[(n+1):(n+ntest),]  Ytest<-Y[(n+1):(n+ntest),]Xtest<-cbind(Xtest,Ytest[,"male"])colnames(Xtest)[ncol(Xtest)]<-"male"Ytest<-Ytest[,1, drop=F]

现在我们计算两种变量重要性度量（由于DRF只实现了drop-and-relearn方法，所以这将需要一些时间）：

# 计算两种度量的变量重要性# 1. Sobol-MDAXY <- as.data.frame(cbind(Xtrain, Ytrain))colnames(XY) <- c(paste('X', 1:(ncol(XY)-1), sep=''), 'Y')num.trees <- 500forest <- sobolMDA::ranger(Y ~., data = XY, num.trees = num.trees, importance = 'sobolMDA')SobolMDA <- forest$variable.importancenames(SobolMDA) <- colnames(Xtrain)# 2. MMD-MDAMMDVimp <- compute_drf_vimp(X=Xtrain,Y=Ytrain,silent=T)print("用于条件期望估计的前10个最重要变量")sort(SobolMDA, decreasing = T)[1:10]print("用于条件分布估计的前5个最重要变量")sort(MMDVimp, decreasing = T)[1:10]

Sobol-MDA:education_level                   age                  male           0.073506769           0.027079349           0.013722756         occupation_11         occupation_43           industry_54           0.013550320           0.010025332           0.007744589           industry_44         occupation_23         occupation_15           0.006657918           0.005772662           0.004610835 marital_never married           0.004545964

MMD-MDA:education_level                   age                  male           0.420316085           0.109212519           0.027356393         occupation_43         occupation_11 marital_never married           0.016861954           0.014122583           0.003449910         occupation_29       marital_married           industry_81           0.002272629           0.002085207           0.001152210           industry_72           0.000984725

在这种情况下，两个变量重要性度量在哪些变量是重要的方面达成了相当大的一致。虽然这不是一种因果分析，但两个度量同时认为对于预测工资来说重要的变量，特别是“age”、“education_level”和“gender”，确实非常重要。

为了得到一小组预测变量，现在可以做如下步骤：

(I) 去除最不重要的变量

(II) 在测试集上计算损失（例如均方误差）

(III) 重新计算剩余变量的重要性

(IV) 重复以上步骤，直到满足某个停止准则

例如，如果损失增加超过5%，可以停止。为了方便在本文中，我只使用保存在“SobolMDA”和“MMDVimp”中的相同的变量重要性值。也就是说，我忽略了步骤（III），只考虑步骤（I）、（II）和（IV）。当估计的目标是完全条件分布时，步骤（II）也不是完全清楚的。我们使用我们在论文中所称的MMD损失，更详细地描述了错误我们在分布预测中所犯的。对于条件均值，我们简单地使用均方误差。这在下面的函数“evalall”中完成：

# 按照存储在SobolMDA和MMDVimp中的重要性值逐个删除变量
evallistSobol <- evalall(SobolMDA, X=Xtrain ,Y=Ytrain ,Xtest, Ytest, metrics=c("MSE"), num.trees )
evallistMMD <- evalall(MMDVimp, X=Xtrain ,Y=Ytrain ,Xtest, Ytest, metrics=c("MMD"), num.trees )
plot(evallistSobol$evalMSE, type="l", lwd=2, cex=0.8, col="darkgreen", main="MSE损失", xlab="删除的变量数", ylab="值")
plot(evallistMMD$evalMMD, type="l", lwd=2, cex=0.8, col="darkgreen", main="MMD损失", xlab="删除的变量数", ylab="值")

结果如下两张图片所示：

请注意，两者都有些起伏的线条，首先是因为我没有重新计算重要性度量，即没有进行步骤（III），第二是由于森林的随机性。除此之外，这些图表很好地展示了随着每个被移除的变量的增加，错误逐渐增加。最不重要的变量的增加速度较慢，而最重要的变量的增加速度则更快，正如人们所预期的那样。特别是，如果删除50个最不重要的变量，两种情况下的损失几乎没有变化！实际上，在两种情况下，可以删除约70个变量，而损失增加不超过6%。然而，必须注意的是，许多预测变量是独热编码的分类变量的一部分，因此在删除预测变量时需要略微谨慎，因为它们对应于一个分类变量的级别。然而，在实际应用中，这可能仍然是可取的。

结论

在本文中，我们探讨了随机森林中变量重要性的现代方法，目标是获得一小组与条件期望和一般情况下的条件分布相关的预测变量或协变量。通过工资数据的例子，我们看到这可以大幅减少使用的预测变量，同时几乎不影响准确性。

如上所述，所提出的度量方法并不严格限于随机森林，而在原则上可以更广泛地应用。然而，随机森林允许优雅的投影方法，可以在不每次重新拟合森林的情况下计算所有变量j的重要性度量！这在[3]和[4]中都有描述。

文献

[1] Breiman, L. (2001). 随机森林. 机器学习, 45(1):5–32.

[2] Breiman, L. (2003a). 随机森林 v3.1 的设置、使用和理解. 技术报告, 加州大学伯克利分校统计系

[3] Bénard, C., Da Veiga, S., and Scornet, E. (2022). 随机森林的平均减少精确度: 不一致性及通过Sobol-MDA的实用解决方案. 生物统计学, 109(4):881–900.

[4] Clément Bénard, Jeffrey Näf, and Julie Josse. 基于MMD的分布随机森林变量重要性, 2023.

[5] Strobl, C., Boulesteix, A.-L., Zeileis, A., and Hothorn, T. (2007). 随机森林变量重要性测度中的偏差: 实例、来源和解决方案. BMC生物信息学, 8:25.

附录：代码

#### compute_drf_vimp.R 的内容 #######' Distributional Random Forests 的变量重要性#'#' @param X 用于训练的输入矩阵.#' @param Y 用于训练的输出矩阵.#' @param X_test 用于测试的输入矩阵。如果为空，则使用包外估计值。#' @param num.trees 需要拟合DRF的树的数量。默认值是500棵树。#' @param silent 如果为FALSE，则打印变量迭代数，否则不打印任何内容。默认为FALSE。#' @return 所有输入变量的重要性值的列表。#' @export#'#' @examplescompute_drf_vimp <- function(X, Y, X_test = NULL, num.trees = 500, silent = FALSE){    # 拟合初始DRF  bandwidth_Y <- drf:::medianHeuristic(Y)  k_Y <- rbfdot(sigma = bandwidth_Y)  K <- kernelMatrix(k_Y, Y, Y)  DRF <- drf(X, Y, num.trees = num.trees)  wall <- predict(DRF, X_test)$weights    # 计算归一化常数  wbar <- colMeans(wall)  wall_wbar <- sweep(wall, 2, wbar, "-")  I0 <- as.numeric(sum(diag(wall_wbar %*% K %*% t(wall_wbar))))    # 逐个剔除变量计算drf重要性  I <- sapply(1:ncol(X), function(j) {    if (!silent){print(paste0('正在计算变量 X', j, ' 的重要性...'))}    DRFj <- drf(X = X[, -j, drop=F], Y = Y, num.trees = num.trees)     DRFpredj <- predict(DRFj, X_test[, -j])    wj <- DRFpredj$weights    Ij <- sum(diag((wj - wall) %*% K %*% t(wj - wall)))/I0    return(Ij)  })    # 计算重新训练的偏差  DRF0 <- drf(X = X, Y = Y, num.trees = num.trees)  DRFpred0 = predict(DRF0, X_test)  w0 <- DRFpred0$weights  vimp0 <- sum(diag((w0 - wall) %*% K %*% t(w0 - wall)))/I0    # 计算最终的重要性（删除偏差并截断负值）  vimp <- sapply(I - vimp0, function(x){max(0,x)})    names(vimp)<-colnames(X)    return(vimp)  }

#### evaluation.R 的内容 ######compute_mmd_loss <- function(Y_train, Y_test, weights){  # Y_train <- scale(Y_train)  # Y_test <- scale(Y_test)  bandwidth_Y <- (1/drf:::medianHeuristic(Y_train))^2  k_Y <- rbfdot(sigma = bandwidth_Y)  K_train <- matrix(kernelMatrix(k_Y, Y_train, Y_train), ncol = nrow(Y_train))  K_cross <- matrix(kernelMatrix(k_Y, Y_test, Y_train), ncol = nrow(Y_train))  weights <- matrix(weights, ncol = ncol(weights))  t1 <- diag(weights%*%K_train%*%t(weights))  t2 <- diag(K_cross%*%t(weights))  mmd_loss <- mean(t1) - 2*mean(t2)  mmd_loss}evalall <- function(Vimp, X ,Y ,Xtest, Ytest, metrics=c("MMD","MSE"), num.trees ){    if (ncol(Ytest) > 1 & "MSE" %in% metrics){    metrics <- metrics[!( metrics %in% "MSE") ]  }    # 根据重要性的增加顺序排序，使最不重要的变量先被移除  Vimp<-sort(Vimp)    if ( is.null(names(Vimp)) ){    stop("需要后续的命名")    }      evalMMD<-matrix(0, nrow=ncol(X))  evalMSE<-matrix(0, nrow=ncol(X))    ###思路：创建一个函数，接受变量重要性测度并进行此循环！！    for (j in 1:ncol(X)){               if (j==1){            if ("MMD" %in% metrics){                DRFred<- drf(X=X,Y=Y)        weights<- predict(DRFred, newdata=Xtest)$weights        evalMMD[j]<-compute_mmd_loss(Y_train=Y, Y_test=Ytest, weights)          }            if ("MSE" %in% metrics){                XY <- as.data.frame(cbind(X, Y))        colnames(XY) <- c(paste('X', 1:(ncol(XY)-1), sep=''), 'Y')        RFfull <- sobolMDA::ranger(Y ~., data = XY, num.trees = num.trees)        XtestRF<-Xtest        colnames(XtestRF) <- paste('X', 1:ncol(XtestRF), sep='')        predRF<-predict(RFfull, data=XtestRF)        evalMSE[j] <- mean((Ytest - predRF$predictions)^2)            }    }else{                  if ("MMD" %in% metrics){                DRFred<- drf(X=X[,!(colnames(X) %in% names(Vimp[1:(j-1)])), drop=F],Y=Y)        weights<- predict(DRFred, newdata=Xtest[,!(colnames(Xtest) %in% names(Vimp[1:(j-1)])), drop=F])$weights        evalMMD[j]<-compute_mmd_loss(Y_train=Y, Y_test=Ytest, weights)              }                        if ("MSE" %in% metrics){                XY <- as.data.frame(cbind(X[,!(colnames(X) %in% names(Vimp[1:(j-1)])), drop=F], Y))        colnames(XY) <- c(paste('X', 1:(ncol(XY)-1), sep=''), 'Y')        RFfull <- sobolMDA::ranger(Y ~., data = XY, num.trees = num.trees)        XtestRF<-Xtest[,!(colnames(Xtest) %in% names(Vimp[1:(j-1)])), drop=F]        colnames(XtestRF) <- paste('X', 1:ncol(XtestRF), sep='')        predRF<-predict(RFfull, data=XtestRF)        evalMSE[j] <- mean((Ytest - predRF$predictions)^2)                # DRFall <- drf(X=X[,!(colnames(X) %in% names(Vimp[1:(j-1)])), drop=F], Y=Y, num.trees=num.trees)        # quantpredictall<-predict(DRFall, newdata=Xtest[,!(colnames(Xtest) %in% names(Vimp[1:(j-1)])), drop=F], functional="

Deep Dives,Interpretable Ml,Machine learning,Math,Random Forest

随机森林中的变量重要性