Relative variable importance for Boosting

I'm looking for an explanation of how relative variable importance is computed in Gradient Boosted Trees that is not overly general/simplistic like:

The measures are based on the number of times a variable is selected for splitting, weighted by the squared improvement to the model as a result of each split, and averaged over all trees. [Elith et al. 2008, A working guide to boosted regression trees]

And that is less abstract than:

$\hat{I_{j}^2}(T)=\sum\limits_{t=1}^{J-1} \hat{i_{t}^2} 1(v_{t}=j)$

Where the summation is over the nonterminal nodes $t$ of the $J$-terminal node tree $T$, $v_{t}$ is the splitting variable associated with node $t$, and $\hat{i_{t}^2}$ is the corresponding empirical improvement in squared error as a result of the split, defined as $i^2(R_{l},R_{r})=\frac{w_{l}w_{r}}{w_{l}+w_{r}}(\bar{y_{l}}-\bar{y_{r}})^2$, where $\bar{y_{l}}, \bar{y_{r}}$ are the left and right daughter response means respectively, and $w_{l}, w_{r}$ are the corresponding sums of the weights. [Friedman 2001, Greedy function approximation: a gradient boosting machine]

Finally, I did not find the Elements of Statistical Learning (Hastie et al. 2008) to be a very helpful read here, as the relevant section (10.13.1 page 367) tastes very similar to the second reference above (which might be explained by the fact that Friedman is a co-author of the book).

PS: I know relative variable importance measures are given by the summary.gbm in the gbm R package. I tried to explore the source code, but I can't seem to find where the actual computation takes place.

Brownie points: I'm wondering how to get these plots in R.

• I just added a new answer to the linked question about how to plot variable importance by class which may be helpful stackoverflow.com/a/51952918/3277050 – see24 Aug 21 '18 at 17:03

1 Answer

I'll use the sklearn code, as it is generally much cleaner than the R code.

Here's the implementation of the feature_importances property of the GradientBoostingClassifier (I removed some lines of code that get in the way of the conceptual stuff)

def feature_importances_(self):
total_sum = np.zeros((self.n_features, ), dtype=np.float64)
for stage in self.estimators_:
stage_sum = sum(tree.feature_importances_
for tree in stage) / len(stage)
total_sum += stage_sum

importances = total_sum / len(self.estimators_)
return importances


This is pretty easy to understand. self.estimators_ is an array containing the individual trees in the booster, so the for loop is iterating over the individual trees. There's one hickup with the

stage_sum = sum(tree.feature_importances_
for tree in stage) / len(stage)


this is taking care of the non-binary response case. Here we fit multiple trees in each stage in a one-vs-all way. Its simplest conceptually to focus on the binary case, where the sum has one summand, and this is just tree.feature_importances_. So in the binary case, we can rewrite this all as

def feature_importances_(self):
total_sum = np.zeros((self.n_features, ), dtype=np.float64)
for tree in self.estimators_:
total_sum += tree.feature_importances_
importances = total_sum / len(self.estimators_)
return importances


So, in words, sum up the feature importances of the individual trees, then divide by the total number of trees. It remains to see how to calculate the feature importances for a single tree.

The importance calculation of a tree is implemented at the cython level, but it's still followable. Here's a cleaned up version of the code

cpdef compute_feature_importances(self, normalize=True):
"""Computes the importance of each feature (aka variable)."""

while node != end_node:
if node.left_child != _TREE_LEAF:
# ... and node.right_child != _TREE_LEAF:
left = &nodes[node.left_child]
right = &nodes[node.right_child]

importance_data[node.feature] += (
node.weighted_n_node_samples * node.impurity -
left.weighted_n_node_samples * left.impurity -
right.weighted_n_node_samples * right.impurity)
node += 1

importances /= nodes[0].weighted_n_node_samples

return importances


This is pretty simple. Iterate through the nodes of the tree. As long as you are not at a leaf node, calculate the weighted reduction in node purity from the split at this node, and attribute it to the feature that was split on

importance_data[node.feature] += (
node.weighted_n_node_samples * node.impurity -
left.weighted_n_node_samples * left.impurity -
right.weighted_n_node_samples * right.impurity)


Then, when done, divide it all by the total weight of the data (in most cases, the number of observations)

importances /= nodes[0].weighted_n_node_samples


It's worth recalling that the impurity is a common name for the metric to use when determining what split to make when growing a tree. In that light, we are simply summing up how much splitting on each feature allowed us to reduce the impurity across all the splits in the tree.

In the context of gradient boosting, these trees are always regression trees (minimize squared error greedily) fit to the gradient of the loss function.

• Thanks a lot for this very detailed answer. Let me some time to carefully go through it before I accept it. – Antoine Jul 30 '15 at 10:38
• While it seems that various impurity criteria can be used, the Gini index was not the criterion used by Friedman. As mentioned in my question and line 878 of your third link, Friedman used the mean squared error impurity criterion with improvement score. If you could update this section of your answer, that would be great. And yes, your're right, it seems that the weights are indeed the number of observations. – Antoine Jul 31 '15 at 6:40
• or maybe it would make your answer even better to keep both the parts about the Gini index and Friedman's original criterion, stressing that the first is used for classification and the second for regression? – Antoine Jul 31 '15 at 6:53
• Antoine, thanks for this update. It is really helpful to know that the mean squared error is the improvement criteria used for the regression trees. It was not obvious how that would be used for classification. However, even in gradient boosting for classification, I think that regression trees are still being used, as opposed to classification trees. At least in python, the regression analysis is being done on the current error at each boosting stage. – Fairly Nerdy Dec 3 '16 at 19:49
• You guys are correct about regression trees. – Matthew Drury Nov 30 '17 at 19:28