Relative variable importance for Boosting

Question

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.

Answer 1

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.