A Simple and Effective Model-Based Variable Importance Measure

05/12/2018 ∙ by Brandon M. Greenwell, et al. ∙ 0

In the era of "big data", it is becoming more of a challenge to not only build state-of-the-art predictive models, but also gain an understanding of what's really going on in the data. For example, it is often of interest to know which, if any, of the predictors in a fitted model are relatively influential on the predicted outcome. Some modern algorithms---like random forests and gradient boosted decision trees---have a natural way of quantifying the importance or relative influence of each feature. Other algorithms---like naive Bayes classifiers and support vector machines---are not capable of doing so and model-free approaches are generally used to measure each predictor's importance. In this paper, we propose a standardized, model-based approach to measuring predictor importance across the growing spectrum of supervised learning algorithms. Our proposed method is illustrated through both simulated and real data examples. The R code to reproduce all of the figures in this paper is available in the supplementary materials.



There are no comments yet.


page 15

page 17

page 19

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Complex supervised learning algorithms, such as neural networks (NNs) and support vector machines (SVMs), are more common than ever in predictive analytics, especially when dealing with large observational databases that don’t adhere to the strict assumptions imposed by traditional statistical techniques (e.g., multiple linear regression which typically assumes linearity, homoscedasticity, and normality). However, it can be challenging to understand the results of such complex models and explain them to management. Graphical displays such as variable importance plots (when available) and partial dependence plots (PDPs)

(Friedman, 2001) offer a simple solution (see, for example, Hastie et al., 2009, pp. 367–380). PDPs are low-dimensional graphical renderings of the prediction function

that allow analysts to more easily understand the estimated relationship between the outcome and predictors of interest. These plots are especially useful in interpreting the output from “black box” models. While PDPs can be constructed for any predictor in a fitted model, variable importance scores are more difficult to define, and when available, their interpretation often depends on the model fitting algorithm used.

In this paper, we consider a standardized method to computing variable importance scores using PDPs. There are a number of advantages to using our approach. First, it offers a standardized procedure to quantifying variable importance across the growing spectrum of supervised learning algorithms. For example, while popular statistical learning algorithms like random forests (RFs) and gradient boosted decision trees (GBMs) have there own natural way of measuring variable importance, each is interpreted differently (these are briefly described in Section 2.1

). Secondly, our method is suitable for use with any trained supervised learning algorithm, provided predictions on new data can be obtained. For example, it is often beneficial (from an accuracy standpoint) to train and tune multiple state-of-the art predictive models (e.g., multiple RFs, GBMs, and deep learning NNs (DNNs)) and then combine them into an ensemble called a

super learner through a process called model stacking. Even if the base learners can provide there own measures of variable importance, there is no logical way to combine them to form an overall score for the super learner. However, since new predictions can be obtained from the super learner, our proposed variable importance measure is still applicable (examples are given in Sections 56). Thirdly, as shown in Section 3.2, our proposed method can be modified to quantify the strength of potential interaction effects. Finally, since our approach is based on constructing PDPs for all the main effects, the analyst is forced to also look at the estimated functional relationship between each feature and the target—which should be done in tandem with studying the importance of each feature.

2 Background

We are often confronted with the task of extracting knowledge from large databases. For this task we turn to various statistical learning algorithms which, when tuned correctly, can have state-of-the-art predictive performance. However, having a model that predicts well is only solving part of the problem. It is also desirable to extract information about the relationships uncovered by the learning algorithm. For instance, we often want to know which predictors, if any, are important by assigning some type of variable importance score to each feature. Once a set of influential features has been identified, the next step is summarizing the functional relationship between each feature, or subset thereof, and the outcome of interest. However, since most statistical learning algorithms are “black box” models, extracting this information is not always straightforward. Luckily, some learning algorithms have a natural way of defining variable importance.

2.1 Model-based approaches to variable importance

Decision trees probably offer the most natural model-based approach to quantifying the importance of each feature. In a binary decision tree, at each node

, a single predictor is used to partition the data into two homogeneous groups. The chosen predictor is the one that maximizes some measure of improvement . The relative importance of predictor is the sum of the squared improvements over all internal nodes of the tree for which was chosen as the partitioning variable; see Breiman et al. (1984) for details. This idea also extends to ensembles of decision trees, such as RFs and GBMs. In ensembles, the improvement score for each predictor is averaged across all the trees in the ensemble. Fortunately, due to the stabilizing effect of averaging, the improvement-based variable importance metric is often more reliable in large ensembles (see Hastie et al., 2009, pg. 368). RFs offer an additional method for computing variable importance scores. The idea is to use the leftover out-of-bag (OOB) data to construct validation-set errors for each tree. Then, each predictor is randomly shuffled in the OOB data and the error is computed again. The idea is that if variable is important, then the validation error will go up when is perturbed in the OOB data. The difference in the two errors is recorded for the OOB data then averaged across all trees in the forest.

In multiple linear regression, the absolute value of the

-statistic is commonly used as a measure of variable importance. The same idea also extends to generalized linear models (GLMs). Multivariate adaptive regression splines (MARS), which were introduced in

Friedman (1991a), is an automatic regression technique which can be seen as a generalization of multiple linear regression and generalized linear models. In the MARS algorithm, the contribution (or variable importance score) for each predictor is determined using a generalized cross-validation (GCV) statistic.

For NNs, two popular methods for constructing variable importance scores are the Garson algorithm (Garson, 1991), later modified by Goh (1995), and the Olden algorithm (Olden et al., 2004)

. For both algorithms, the basis of these importance scores is the network’s connection weights. The Garson algorithm determines variable importance by identifying all weighted connections between the nodes of interest. Olden’s algorithm, on the other hand, uses the product of the raw connection weights between each input and output neuron and sums the product across all hidden neurons. This has been shown to outperform the Garson method in various simulations. For DNNs, a similar method due to

Gedeon (1997) considers the weights connecting the input features to the first two hidden layers (for simplicity and speed); but this method can be slow for large networks.

2.2 Filter-based approaches to variable importance

Filter-based approaches, which are described in Kuhn & Johnson (2013), do not make use of the fitted model to measure variable importance. They also do not take into account the other predictors in the model.

For regression problems, a popular approach to measuring the variable importance of a numeric predictor is to first fit a flexible nonparametric model between and the target ; for example, the locally-weighted polynomial regression (LOWESS) method developed by Cleveland (1979). From this fit, a pseudo- measure can be obtained from the resulting residuals and used as a measure of variable importance. For categorical predictors, a different method based on standard statistical tests (e.g., -tests and ANOVAs) is employed; see Kuhn & Johnson (2013) for details.

For classification problems, an area under the ROC curve (AUC) statistic can be used to quantify predictor importance. The AUC statistic is computed by using the predictor as input to the ROC curve. If can reasonably separate the classes of , that is a clear indicator that is an important predictor (in terms of class separation) and this is captured in the corresponding AUC statistic. For problems with more than two classes, extensions of the ROC curve or a one-vs-all approach can be used.

2.3 Partial dependence plots

Harrison & Rubinfeld (1978) analyzed a data set containing suburban Boston housing data from the 1970 census. They sought a housing value equation using an assortment of features; see Table IV of Harrison & Rubinfeld (1978)

for a description of each variable. The usual regression assumptions, such as normality, linearity, and constant variance, were clearly violated, but through an exhausting series of transformations, significance testing, and grid searches, they were able to build a model which fit the data reasonably well (

). Their prediction equation is given in Equation (1). This equation makes interpreting the model easier. For example, the average number of rooms per dwelling () is included in the model as a quadratic term with a positive coefficient. This means that there is a monotonic increasing relationship between and the predicted median home value, but larger values of have a greater impact.


However, classical regression and model building is rather ill-suited for more contemporary big data sets, like the Ames housing data described in Cock (2011)

which has a total of 79 predictors (and many more that can be created through feature engineering). Fortunately, using modern computing power, many supervised learning algorithms can fit such data sets in seconds, producing powerful, highly accurate models. The downfall of many of these machine learning algorithms, however, is decreased interpretability. For example, fitting a well-tuned RF to the Boston housing data will likely produce a model with more accurate predictions, but no interpretable prediction formula such as the one in Equation (


To help understand the estimated functional relationship between each predictor and the outcome of interest in a fitted model, we can construct PDPs. PDPs are particularly effective at helping to explain the output from “black box” models, such as RFs and SVMs. Not only do PDPs visually convey the relationship between low cardinality subsets of the feature set (usually 1-3) and the response (while accounting for the average effect of the other predictors in the model), they can also be used to rank and score the predictors in terms of their relative influence on the predicted outcome, as will be demonstrated in this paper.

Let represent the predictors in a model whose prediction function is . If we partition into an interest set, , and its complement, , then the “partial dependence” of the response on is defined as


where is the marginal probability density of : . Equation (2) can be estimated from a set of training data by


where are the values of that occur in the training sample; that is, we average out the effects of all the other predictors in the model.

Constructing a PDP (3) in practice is rather straightforward. To simplify, let

be the predictor variable of interest with unique values

. The partial dependence of the response on can be constructed as follows:

Input: the unique predictor values ; Output: the estimated partial dependence values . for  do
       (1) copy the training data and replace the original values of with the constant ; (2) compute the vector of predicted values from the modified copy of the training data; (3) compute the average prediction to obtain .
end for
The PDP for is obtained by plotting the pairs for .
Algorithm 1 A simple algorithm for constructing the partial dependence of the response on a single predictor .

Algorithm 1 can be computationally expensive since it involves

passes over the training records. Fortunately, it is embarrassingly parallel and computing partial dependence functions for each predictor can be done rather quickly on a machine with a multi-core processor. For large data sets, it may be worthwhile to reduce the grid size by using specific quantiles for each predictor, rather than all the unique values. For example, the partial dependence function can be approximated very quickly by using the deciles of the unique predictor values. The exception is classification and regression trees based on single-variable splits which can make use of the efficient weighted tree traversal method described in

Friedman (2001).

While PDPs are an invaluable tool in understanding the relationships uncovered by complex nonparametric models, they can be misleading in the presence of substantial interaction effects (Goldstein et al., 2015). To overcome this issue, Goldstein et al. introduced the concept of individual conditional expectation (ICE) curves. ICE curves display the estimated relationship between the response and a predictor of interest for each observation; in other words, skipping step 1 (c) in Algorithm 1. Consequently, the PDP for a predictor of interest can be obtained by averaging the corresponding ICE curves across all observations. Although ICE curves provide a refinement over traditional PDPs in the presence of substantial interaction effects, in Section 3.2, we show how to use partial dependence functions to evaluate the strength of potential interaction effects.

2.4 The Ames housing data set

For illustration, we will use the Ames housing data set—a modernized and expanded version of the often cited Boston Housing data set. These data are available in the AmesHousing package (Kuhn, 2017). Using the R package h2o (The H2O.ai team, 2017), we trained and tuned a GBM using 10-fold cross-validation. The model fit is reasonable, with a cross-validated (pseudo) of 91.54%. Like other tree-based ensembles, GBMs have a natural way of defining variable importance which was described in Section 2.1. The variable importance scores for these data are displayed in the left side of Figure 1. This plot indicates that the overall quality of the material and finish of the house (Overall_Qual), physical location within the Ames city limits (Neighborhood), and the above grade (ground) living area (Gr_Liv_Area) are highly associated with the logarithm of the sales price (Log_Sale_Price). The variable importance scores also indicate that pool quality (Pool_QC), the number of kitchens above grade (Kitchen_AbvGr), and the low quality finished square feet for all floors (Low_Qual_Fin_SF) have little association with Log_Sale_Price. (Note that the bottom two features in Figure 1Street and Utilities—have a variable importance score of exactly zero; in other words, they were never used to partion the data at any point in the GBM ensemble.)

Figure 1: Variable importance scores for the Ames housing data. Left: Traditional approach. Right: Partial dependence-based approach.

The PDPs for these six variables are displayed in Figure 2. These plots indicate that Overall_Qual, Neighborhood, and Gr_Liv_Area have a strong nonlinear relationship with the predicted outcome. For instance, it seems that GrLivArea has a monotonically increasing relationship with Log_Sale_Price until about 12 sq ft, after which the relationship flattens out. Notice how the PDPs for Pool_QC, Kitchen_AbvGr, and Low_Qual_Fin_SF are relatively flat in comparison. It is this notion of “flatness” which we will use as a basis to define our variable importance measure.

Figure 2: Partial dependence of log selling price on the three highest (top row) and lowest (bottom row) ranked predictors. The dashed red line in each plot represents the mean of the log selling prices

3 A partial dependence-based variable importance measure

The PDPs for Pool_QC, Kitchen_AbvGr, and Low_Qual_Fin_SF in Figure 2 are relatively flat, indicating that they do not have much influence on the predicted value of Log_Sale_Price. In other words, the partial dependence values display little variability. One might conclude that any variable for which the PDP is “flat” is likely to be less important than those predictors whose PDP varies across a wider range of the response.

Our notion of variable importance is based on any measure of the “flatness” of the partial dependence function. In general, we define

where is any measure of the “flatness” of

. A simple and effective measure to use is the sample standard deviation for continuous predictors and the range statistic divided by four for factors with

levels; the range divided by four provides an estimate of the standard deviation for small to moderate sample sizes. Based on Algorithm 1, our importance measure for predictor is simply


Note that our variable importance metric relies on the fitted model; hence, it is crucial to properly tune and train the model to attain the best performance possible.

To illustrate, we applied Algorithm 1 to all of the predictors in the Ames GBM model and computed (4). The results are displayed in right side of Figure 1. In this case, our partial dependence-based algorithm matches closely with the results from the GBM. In particular, Figure LABEL:fig:ames-gbm-vi-both indicates that Overall_Qual, Neighborhood, and Gr_Liv_Area are still the most important variables in predicting Log_Sale_Price; though, Neighborhood and Gr_Liv_Area have swapped places.

3.1 Linear models

As mentioned earlier, a natural choice for measuring the importance of each term in a linear model is to use the absolute value of the corresponding coefficient divided by its estimated standard error (i.e., the absolute value of the

-statistic). This turns out to be equivalent to the partial dependence-based metric (4

) when the predictors are independently and uniformly distributed over the same range.

For example, suppose we have a linear model of the form

where () is a constant, and are both independent random variables, and . Since we know the distribution of and , we can easily find and . For instance,

where . Simple calculus then leads to

Because is additive, the true partial dependence functions are just simple linear regressions in each predictor with their original coefficient and an adjusted intercept. Taking the variance of each gives

Hence, the standard deviations are just the absolute values of the original coefficients (scaled by the same constant).

To illustrate, we simulated observations from the following linear model

where and are both independent random variables, and . For this example, we have

These are plotted as red lines in Figure 3. Additionally, the black lines in Figure 3 correspond to the estimated partial dependence functions using Algorithm 1.

Figure 3: Estimated (black) and true (red) partial dependence functions from a linear model with two predictors.

Based on these plots, is more influential than . Taking the absolute value of the ratio of the slopes in and gives . In other words, is roughly 1.67 times more influential on than . Using the partial-dependence-based variable importance metric, we obtain and which gives the ratio . In fact, we can compute the ratio of the true variances of and :

Taking the square root gives .

Using the absolute value of the -statistic becomes less useful in linear models when, for example, a predictor appears in multiple terms (e.g., interaction effects and polynomial terms). The partial dependence approach, on the other hand, does not suffer from such drawbacks.

3.2 Detecting interaction effects

As it turns out, our partial dependence-based variable importance measure (4) can also be used to quantify the strength of potential interaction effects. Let be the standard deviation of the joint partial dependence values for and . Essentially, a weak interaction effect of and on would suggest that has little variation when either or is held constant while the other varies.

Let , , be any two predictors in the feature space . Construct the partial dependence function and compute for each unique value of , denoted , and take the standard deviation of the resulting importance scores. The same can be done for and the results are averaged together. Large values (relative to each other) would be indicative of possible interaction effects.

4 Friedman’s regression problem

To further illustrate, we will use one of the regression problems described in Friedman (1991b) and Breiman (1996). The feature space consists of ten independent random variables; however, only five out of these ten actually appear in the true model. The response is related to the features according to the formula

where . Using the R package nnet (Venables & Ripley, 2002), we fit a NN with one hidden layer containing eight units and a weight decay of 0.01 (these parameters were chosen using 5-fold cross-validation) to 500 observations simulated from the above model with . The cross-validated value was 0.94.

Variable importance plots are displayed in Figure 4. Notice how the Garson and Olden algorithms incorrectly label some of the features not in the true model as “important”. For example, the Garson algorithm incorrectly labels (which is not included in the true model) as more important than (which is in the true model). Similarly, Oden’s method incorrectly labels as being more important than . Our method, on the other hand, clearly labels all five of the predictors in the true model as the most important features in the fitted NN.

Figure 4: Variable importance plots for the NN fit to the Friedman regression data. Left: partial dependence-based method. Middle: Garson’s method. Right: Olden’s method.

We also constructed the partial dependence functions for all pairwise interactions and computed the interaction statistic discussed in Section 3.2. The top ten interaction statistics are displayed in Figure 5. There is a clear indication of an interaction effect between features and , the only interaction effect present in the true model.

Figure 5: Variable importance-based interaction statistics from the NN fit to the Friedman regression data set.

In fact, since we know the distributions of the predictors in the true model, we can work out the true partial dependence functions. For example, for the pairs and , we have


Next, we simulated the standard deviation of for a wide range of fixed values of ; this is what is trying to estimate. The results from doing this for both predictors in each model are displayed in Figure~6. The top row of Figure~6 illustrates that the importance of (i.e., the strength of its relationship to the predicted outcome) heavily depends on the value of and vice versa (i.e., an interaction effect between and ). In the bottom row, on the other hand, we see that the importance of does not depend on the value of and vice versa (i.e., no interaction effect between and ).

Figure 6: Results from a small Monte Carlo simulation on the interaction effects between and (top row), and and (bottom row).

4.1 Friedman’s -statistic

An alternative measure for the strength of interaction effects is known as Friedman’s -statistic (Friedman & Popescu, 2008). Coincidentally, this method is also based on the estimated partial dependence functions of the corresponding predictors, but uses a different approach.

For comparison, we fit a GBM to the Friedman regression data from the previous section. The parameters were chosen using 5-fold cross-validation. We used the R package gbm (Ridgeway, 2017) which has built-in support for computing Friedman’s -statistic for any combination of predictors. The results are displayed in Figure 7. To our surprise, the -statistic did not seem to catch the true interaction between and . Instead, the -statistic ranked the pairs and as having the strongest interaction effects, even though these predictors do not appear in the true model. Our variable importance-based interaction statistic, on the other hand, clearly suggests the pair as having the strongest interaction effect.

Figure 7: Interaction statistics for the GBM model fit to the Friedman regression data. Left: Friedman’s -statistic. Right: Our variable importance-based interaction statistic.

5 Application to model stacking

In the Ames housing example, we used a GBM to illustrate the idea of using the partial dependence function to quantify the importance of each predictor on the predicted outcome (in this case, Log_Sale_Price). While the GBM model achieved a decent cross-validated of 91.54%, better predictive performance can often be attained by creating an ensemble of learning algorithms. While GBMs are themselves ensembles, we can use a method called stacking to combine the GBM with other cutting edge learning algorithms to form a super learner (Wolpert, 1992). Such stacked ensembles tend to outperform any of the individual base learners (e.g., a single RF or GBM) and have been shown to represent an asymptotically optimal system for learning (van der Laan et al., 2003).

For the Ames data set, in addition to the GBM, we trained and tuned a RF using 10-fold cross-validation. The cross-validated predicted values from each of these models were combined to form a matrix, where is the number of training records. This matrix, together with the observed values of Log_Sale_Price, formed the “level-one” data. Next, we trained a metalearning algorithm—in this case a GLM—on the level-one data to create a stacked ensemble. To generate new predictions, we first generate predictions from the RF and GBM learners, then feed those into the GLM metalearner to generate the ensemble prediction.

Even though the base learners in this example have the built-in capability to compute variable importance scores, there is no way of constructing them from the super learner. However, since we can generate predictions from the super learner, we can easily construct partial dependence functions for each predictor and use Equation (4). The left and middle plots in Figure 8 display the top 15 variable importance scores from the individual base learners. The right plot displays the top 15 variable importance scores constructed using Equation (4) for the combined super learner. All models seem to agree on the top three predictors of Log_Sale_Price; namely, Overall_Qual, Neighborhood, and Gr_Liv_Area.

Figure 8: Variable importance plots for the Ames data set. Left: RF. Middle: GBM. Right: Stacked ensemble (i.e., super learner).

6 Application to automatic machine learning

The data science field has seen an explosion in interest over the last decade. However, the supply of data scientists and machine learning experts has not caught up with the demand. Consequently, many organizations are turning to automated machine learning (AutoML) approaches to predictive modelling. AutoML has been a topic of increasing interest over the last couple of years and open source implementations—like those in H2O and auto-sklearn

(Feurer et al., 2015)—have made it a simple and viable options for real supervised learning problems.

Current AutoML algorithms leverage recent advantages in Bayesian optimization, meta-learning, and ensemble construction (i.e., model stacking). The benefit, as compared to the stacked ensemble formed in Section 5

, is that AutoML frees the analyst from having to do algorithm selection (e.g., “Do I fit an RF or GBM to my data?”) and hyperparameter tuning (e.g., how many hidden layers to use in a DNN). While AutoML has great practical potential, it does not automatically provide any useful interpretations—AutoML is

not automated data science (Mayo, 2017). For instance, which variables are the most important in making accurate predictions? How do these variables functionally related to the outcome of interest? Fortunately, our approach can still be used to answer these two questions simultaneously.

To illustrate, we used the R implementation of H2O’s AutoML algorithm to model the airfoil self-noise data set (Lichman, 2013) which are available from the University of California Machine Learning Repository: https://archive.ics.uci.edu/ml/datasets/airfoil+self-noise. These data are from a NASA experiment investigating different size NACA 0012 airfoils at various wind tunnel speeds and angles of attack. The objective is to accurately predict the scaled sound pressure level (dB) (scaled_sound_pressure_level) using frequency (Hz) (frequency), angle of attack (degrees) (angle_of_attack), chord length (m) (chord_length), free-stream velocity (m/s) (free_stream_velocity), and suction side displacement thickness (m) (suction_side_displacement_thickness). H2O’s current AutoML implementation trains and cross-validates an RF, an extremely-randomized forest (XRF) (Geurts et al., 2006), a random grid of GBMs, a random grid of DNNs, and then trains a stacked ensemble using the approach outlined in Section 5. We used 10-fold cross-validation and RMSE for the validation metric. Due to hardware constraints, we used a max run time of 15 minutes (the default is one hour). The final model consisted of one DNN, a random grid of 37 GBMs, one GLM, one RF, one XRT, and two stacked ensembles. The final stacked ensemble achieved a 10-fold cross-validated RMSE and R-squared of 1.43 and 95.68, respectively.

Since predictions can be obtained from the automated stacked ensemble, we can easily apply Algorithm 1 and construct PDPs for all the features; these are displayed in Figure 9. It seems frequency (Hz) and suction side displacement thickness (m) have a strong monotonically decreasing relationship with scaled sound pressure level (dB). They also appear to be more influential than the other three. Using Equation (4), we computed the variable importance of each of the five predictors which are displayed in Table 1.

Figure 9: PDPs for the five features in the airfoil self-noise data set based on an AutoML algorithm consisting of DNNs, GBMs, GLMs, RFs, XRTs, and stacked ensembles. The dashed red line in each plot represents the mean scaled sound pressure level (dB) in the training data.
Variable Importance
Angle of attack (degrees) 0.4888482
Free stream velocity (m/s) 1.1158569
Chord length (m) 1.4168613
Suction side displacement thickness (m) 3.3389105
Frequency (Hz) 5.4821362
Table 1: Partial dependence-based variable importance scores for the five predictors in the airfoil self-noise data set based on an AutoML algorithm consisting of DNNs, GBMs, GLMs, RFs, XRTs, and stacked ensembles.

7 Discussion

We have discussed a new variable importance measure that is (i) suitable for use with any supervised learning algorithm, provided new predictions can be obtained, (ii) model-based and takes into account the effect of all the features in the model, (iii) consistent and has the same interpretation regardless of the learning algorithm employed, and (iv) has the potential to help identify possible interaction effects. Since our algorithm is model-based it requires that the model be properly trained and tuned to achieve optimum performance. While this new approach appears to have high utility, more research is needed to determine where its deficiencies may lie. For example, outliers in the feature space can cause abnormally large fluctuations in the partial dependence values

. Therefore, it may be advantageous to use more robust measures of spread to describe the variability in the estimated partial dependence values; a reasonable choice would be the median absolute deviation which has a finite sample breakdown point of . It is also possible to replace the mean in step (3) of Algorithm 1 with a more robust estimate such as the median or trimmed mean. Another drawback is the computational burden imposed by Algorithm 1 on large data sets, but this can be mitigated using the methods discussed in Greenwell (2017).

All the examples in this article were produced using R version 3.4.0 (R Core Team, 2017); a software environment for statistical computing. With the exception of Figure 6, all graphics were produced using the R package ggplot2 (Wickham, 2009); Figure 6 was produced using the R package lattice (Sarkar, 2008).


R package:

The R package ‘vip‘, hosted on GitHub at https://github.com/AFIT-R/vip, contains functions for computing variable importance scores and constructing variable importance plots for various types of fitted models in R using the partial dependence-based approach discussed in this paper.

R code:

The R script ‘vip-2018.R‘ contains the R code to reproduce all of the results and figures in this paper.


  • (1)
  • Breiman (1996) Breiman, L. (1996), ‘Bagging predictors’, Machine Learning 8(2), 209–218.
  • Breiman et al. (1984) Breiman, L., Friedman, J. & Charles J. Stone, R. A. O. (1984), Classification and Regression Trees, The Wadsworth and Brooks-Cole statistics-probability series, Taylor & Francis.
  • Cleveland (1979) Cleveland, W. S. (1979), ‘Robust locally weighted regression and smoothing scatterplots’, Journal of the American Statistical Association 74(368), 829–836.
  • Cock (2011) Cock, D. D. (2011), ‘Ames, iowa: Alternative to the boston housing data as an end of semester regression project’, Journal of Statistics Education 19(3), 1–15.
  • Feurer et al. (2015) Feurer, M., Klein, A., Eggensperger, K., Springenberg, J., Blum, M. & Hutter, F. (2015), Efficient and robust automated machine learning, in C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama & R. Garnett, eds, ‘Advances in Neural Information Processing Systems 28’, Curran Associates, Inc., pp. 2962–2970.
  • Friedman (1991a) Friedman, J. H. (1991a), ‘Multivariate adaptive regression splines’, The Annals of Statistics 19(1), 1–67.
  • Friedman (1991b) Friedman, J. H. (1991b), ‘Multivariate adaptive regression splines’, The Annals of Statistics 19(1), 1–67.
  • Friedman (2001) Friedman, J. H. (2001), ‘Greedy function approximation: A gradient boosting machine’, The Annals of Statistics 29, 1189–1232.
  • Friedman & Popescu (2008) Friedman, J. H. & Popescu, B. E. (2008), ‘Predictive learning via rule ensembles’, Annals of Applied Statistics 2(3), 916–954.
  • Garson (1991) Garson, D. G. (1991), ‘Interpreting neural-network connection weights’, Artificial Intelligence Expert 6(4), 46–51.
  • Gedeon (1997) Gedeon, T. (1997), ‘Data mining of inputs: Analysing magnitude and functional measures’, International Journal of Neural Systems 24(2), 123–140.
  • Geurts et al. (2006) Geurts, P., Ernst, D. & Wehenkel, L. (2006), ‘Extremely randomized trees’, Machine Learning 63(1), 3–42.
  • Goh (1995) Goh, A. (1995), ‘Back-propagation neural networks for modeling complex systems’, Artificial Intelligence in Engineering 9(3), 143–151.
  • Goldstein et al. (2015) Goldstein, A., Kapelner, A., Bleich, J. & Pitkin, E. (2015), ‘Peeking inside the black box: Visualizing statistical learning with plots of individual conditional expectation’, Journal of Computational and Graphical Statistics 24(1), 44–65.
  • Greenwell (2017) Greenwell, B. M. (2017), ‘pdp: An r package for constructing partial dependence plots’, The R Journal 9(1), 421–436.
  • Harrison & Rubinfeld (1978) Harrison, D. & Rubinfeld, D. L. (1978), ‘Hedonic housing prices and the demand for clean air’, Journal of Environmental Economics and Management 5(1), 81–102.
  • Hastie et al. (2009) Hastie, T., Tibshirani, R. & Friedman, J. (2009), The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition, Springer Series in Statistics, Springer-Verlag.
  • Kuhn (2017) Kuhn, M. (2017), AmesHousing: The Ames Iowa Housing Data. R package version 0.0.3.
  • Kuhn & Johnson (2013) Kuhn, M. & Johnson, K. (2013), Applied Predictive Modeling, SpringerLink : Bücher, Springer New York.
  • Lichman (2013) Lichman, M. (2013), ‘UCI machine learning repository’.
  • Mayo (2017) Mayo, M. (2017), ‘The current state of automated machine learning’.
  • Olden et al. (2004) Olden, J. D., Joy, M. K. & Death, R. G. (2004), ‘An accurate comparison of methods for quantifying variable importance in artificial neural networks using simulated data’, Ecological Modelling 178(3), 389–397.
  • R Core Team (2017) R Core Team (2017), R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
  • Ridgeway (2017) Ridgeway, G. (2017), gbm: Generalized Boosted Regression Models. R package version 2.1.3.
  • Sarkar (2008) Sarkar, D. (2008),

    Lattice: Multivariate Data Visualization with R

    , Springer, New York.
    ISBN 978-0-387-75968-5.
  • The H2O.ai team (2017) The H2O.ai team (2017), h2o: R Interface for H2O. R package version
  • van der Laan et al. (2003) van der Laan, M. J., Polley, E. C. & Hubbard, A. E. (2003), ‘Super learner’, Statistical Applications in Genetics and Molecular Biology 6(1).
  • Venables & Ripley (2002) Venables, W. N. & Ripley, B. D. (2002), Modern Applied Statistics with S, 4th edn, Springer-Verlag, New York.
  • Wickham (2009) Wickham, H. (2009), ggplot2: Elegant Graphics for Data Analysis, Springer-Verlag New York.
  • Wolpert (1992) Wolpert, D. H. (1992), ‘Stacked generalization’, Neural Networks 5, 241–259.