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 lowdimensional graphical renderings of the prediction functionthat 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 stateofthe 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 5–6). 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 stateoftheart 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 Modelbased approaches to variable importance
Decision trees probably offer the most natural modelbased 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 improvementbased 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 outofbag (OOB) data to construct validationset 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 crossvalidation (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 Filterbased approaches to variable importance
Filterbased 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 locallyweighted 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 onevsall 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.(1)  
However, classical regression and model building is rather illsuited 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 welltuned 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 (
1).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 13) 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
(2) 
where is the marginal probability density of : . Equation (2) can be estimated from a set of training data by
(3) 
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: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 multicore 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 singlevariable 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 10fold crossvalidation. The model fit is reasonable, with a crossvalidated (pseudo) of 91.54%. Like other treebased 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 1—Street 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.)
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.
3 A partial dependencebased 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(4) 
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 dependencebased algorithm matches closely with the results from the GBM. In particular, Figure LABEL:fig:amesgbmviboth 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 dependencebased 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.
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 partialdependencebased 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 dependencebased 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 5fold crossvalidation) to 500 observations simulated from the above model with . The crossvalidated 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.
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.
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
and
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 ).
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 5fold crossvalidation. We used the R package gbm (Ridgeway, 2017) which has builtin 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 importancebased interaction statistic, on the other hand, clearly suggests the pair as having the strongest interaction effect.
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 crossvalidated 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 10fold crossvalidation. The crossvalidated 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 “levelone” data. Next, we trained a metalearning algorithm—in this case a GLM—on the levelone 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 builtin 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.
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 autosklearn
(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, metalearning, 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 selfnoise data set (Lichman, 2013) which are available from the University of California Machine Learning Repository: https://archive.ics.uci.edu/ml/datasets/airfoil+selfnoise. 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), freestream velocity (m/s) (free_stream_velocity), and suction side displacement thickness (m) (suction_side_displacement_thickness). H2O’s current AutoML implementation trains and crossvalidates an RF, an extremelyrandomized 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 10fold crossvalidation 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 10fold crossvalidated RMSE and Rsquared 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.
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 
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) modelbased 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 modelbased 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).
SUPPLEMENTARY MATERIAL
 R package:

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

The R script ‘vip2018.R‘ contains the R code to reproduce all of the results and figures in this paper.
References
 (1)

Breiman (1996)
Breiman, L. (1996), ‘Bagging predictors’,
Machine Learning 8(2), 209–218.
https://doi.org/10.1023/A:1018054314350  Breiman et al. (1984) Breiman, L., Friedman, J. & Charles J. Stone, R. A. O. (1984), Classification and Regression Trees, The Wadsworth and BrooksCole statisticsprobability 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.
http://ww2.amstat.org/publications/jse/v19n3/decock.pdf 
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.
http://papers.nips.cc/paper/5872efficientandrobustautomatedmachinelearning.pdf 
Friedman (1991a)
Friedman, J. H. (1991a), ‘Multivariate
adaptive regression splines’, The Annals of Statistics 19(1), 1–67.
https://doi.org/10.1214/aos/1176347963 
Friedman (1991b)
Friedman, J. H. (1991b), ‘Multivariate
adaptive regression splines’, The Annals of Statistics 19(1), 1–67.
https://doi.org/10.1214/aos/1176347963 
Friedman (2001)
Friedman, J. H. (2001), ‘Greedy function
approximation: A gradient boosting machine’, The Annals of Statistics
29, 1189–1232.
https://doi.org/10.1214/aos/1013203451 
Friedman & Popescu (2008)
Friedman, J. H. & Popescu, B. E. (2008), ‘Predictive learning via rule ensembles’, Annals
of Applied Statistics 2(3), 916–954.
https://doi.org/10.1214/07aoas148  Garson (1991) Garson, D. G. (1991), ‘Interpreting neuralnetwork 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.
https://doi.org/10.1007/s1099400662261 
Geurts et al. (2006)
Geurts, P., Ernst, D. & Wehenkel, L. (2006), ‘Extremely randomized trees’, Machine Learning
63(1), 3–42.
https://doi.org/10.1007/s1099400662261  Goh (1995) Goh, A. (1995), ‘Backpropagation 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.
https://doi.org/10.1080/10618600.2014.907095 
Greenwell (2017)
Greenwell, B. M. (2017), ‘pdp: An r package
for constructing partial dependence plots’, The R Journal 9(1), 421–436.
https://journal.rproject.org/archive/2017/RJ2017016/index.html 
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.
https://doi.org/10.1016/00950696(78)900062  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, SpringerVerlag.

Kuhn (2017)
Kuhn, M. (2017), AmesHousing: The Ames
Iowa Housing Data.
R package version 0.0.3.
https://CRAN.Rproject.org/package=AmesHousing  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’.
http://archive.ics.uci.edu/ml 
Mayo (2017)
Mayo, M. (2017), ‘The current state of
automated machine learning’.
https://www.kdnuggets.com/2017/01/currentstateautomatedmachinelearning.html  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.
https://www.Rproject.org/ 
Ridgeway (2017)
Ridgeway, G. (2017), gbm: Generalized
Boosted Regression Models.
R package version 2.1.3.
https://CRAN.Rproject.org/package=gbm 
Sarkar (2008)
Sarkar, D. (2008),
Lattice: Multivariate Data Visualization with R
, Springer, New York. ISBN 9780387759685.
http://lmdvr.rforge.rproject.org 
The H2O.ai team (2017)
The H2O.ai team (2017), h2o: R
Interface for H2O.
R package version 3.14.0.3.
https://CRAN.Rproject.org/package=h2o  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,
SpringerVerlag, New York.
http://www.stats.ox.ac.uk/pub/MASS4 
Wickham (2009)
Wickham, H. (2009), ggplot2: Elegant
Graphics for Data Analysis, SpringerVerlag New York.
http://ggplot2.org  Wolpert (1992) Wolpert, D. H. (1992), ‘Stacked generalization’, Neural Networks 5, 241–259.
Comments
There are no comments yet.