Interpretability: Cracking open the black box, Part 2
The second part in a series on leveraging techniques to take a look inside the black box of AI, this guide considers posthoc interpretation that is useful when the model is not transparent.
By Manu Joseph, Problem Solver, Practitioner, Researcher at Thoucentric Analytics.
In the last post in the series, we defined what interpretability is and looked at a few interpretable models and the quirks and ‘gotchas’ in it. Now let’s dig deeper into the posthoc interpretation techniques which is useful when your model itself is not transparent. This resonates with most realworld use cases because whether we like it or not, we get better performance with a black box model.
Data Set
For this exercise, I have chosen the Adult dataset a.k.a Census Income dataset. Census Income is a pretty popular dataset which has demographic information like age, occupation, along with a column which tells us if the income of the particular person >50k or not. We are using this column to run a binary classification using Random Forest. The reasons for choosing Random Forest are twofold:
 Random Forest is one of the most popularly used algorithm, along with Gradient Boosted Trees. Both of them are from the family of ensemble algorithms with Decision Trees.
 There are a few techniques that are specific to treebased models, which I want to discuss.
Overview of Census Income Dataset.
Sample Data from Census Income Dataset.
Posthoc Interpretation
Now let’s look at techniques to do posthoc interpretation to understand our black box models. All through the rest of the blog, the discussion will be based on machine learning models(and not deep learning) and will be based on structured data. While many of the methods here are model agnostic, since there are a lot of specific ways to interpret deep learning models, especially on unstructured data, we leave that out of our current scope.(Maybe another blog, another day.)
DATA PREPROCESSING
 Encoded the target variable into numerical variables
 Dealt with missing values
 Transformed marital_statusinto a binary variable by combining a few values
 Dropped educationbecause education_num gives the same info, but in numerical format
 Dropped capital_gainand capital_loss because they do not have any information. More than 90% of them are zeroes
 Dropped native_countrybecause of high cardinality and skew towards US
 Dropped relationshipbecause of a lot of overlap with marital_status
A Random Forest algorithm was tuned and trained on the data with 83.58% performance. It is a decent score considering the best scores vary from 7886% based on the way you model and test set. But for our purposes, the model performance is more than enough.
1) Mean Decrease in Impurity
This is by far the most popular way of explaining a treebased model and it’s ensembles. A lot of it is because of SciKit Learn, and its easy to use implementation of the same. Fitting a Random Forest or a Gradient Boosting Model and plotting the “feature importance” has become the most used and abused technique among Data Scientist.
The mean decrease in impurity importance of a feature is computed by measuring how effective the feature is at reducing uncertainty (classifiers) or variance (regressors) when creating decision trees within any ensemble Decision Tree method(Random Forest, Gradient Boosting, etc.).
The advantages of the technique are:
 A fast and easy way of getting feature importance
 Readily available in Scikit Learn and Decision Tree implementation in R
 It is pretty intuitive to explain to a layman
ALGORITHM
 During tree construction, whenever a split is made, we keep track of which feature made the split, what was the Gini impurity before and after, and how many samples did it affect
 At the end of the tree building process, you calculate the total gain in Gini Index attributed to each feature
 And in case of a Random Forest or a Gradient Boosted Trees, we average this score over all the trees in the ensemble
IMPLEMENTATION
Scikit Learn implements this by default in the “feature importance” in treebased models. So retrieving them and plotting the top 25 features is very simple.
feat_imp = pd.DataFrame({'features': X_train.columns.tolist(), "mean_decrease_impurity": rf.feature_importances_}).sort_values('mean_decrease_impurity', ascending=False) feat_imp = feat_imp.head(25) feat_imp.iplot(kind='bar', y='mean_decrease_impurity', x='features', yTitle='Mean Decrease Impurity', xTitle='Features', title='Mean Decrease Impurity', )
Click here for full interactive plot
We can also retrieve and plot the mean decrease in the impurity of each tree as a box plot.
# get the feature importances from each tree and then visualize the # distributions as boxplots all_feat_imp_df = pd.DataFrame(data=[tree.feature_importances_ for tree in rf], columns=X_train.columns) order_column = all_feat_imp_df.mean(axis=0).sort_values(ascending=False).index.tolist() all_feat_imp_df[order_column[:25]].iplot(kind='box', xTitle = 'Features', yTitle='Mean Decease Impurity')
Click for full interactive plot
INTERPRETATION
 The top 4 features are marital status, education_num, age,and hours_worked. This makes perfect sense, as they have a lot to do with how much you earn
 Notice the two features fnlwgt and random in there? Are they more important than the occupation of a person?
 One other caveat here is that we are looking at the onehot features as separate features and it may have some bearing on why the occupation features are ranked lower than random. Dealing with Onehot features when looking at feature importance is a whole other topic
Let’s take a look at what fnlwgt and random are.
 The description of the dataset for fnlwgt is a long and convoluted description of how the census agency uses sampling to create “weighted tallies” of any specified socioeconomic characteristics of the population. In short, it is a sampling weight and nothing to do with how much a person earns
 And randomis just what the name says. Before fitting the model, I made a column with random numbers and called it random
Now, surely, these features cannot be more important than other features like occupation, work_class, sex etc. If that is the case, then something is wrong.
THE JOKER IN THE PACK A.K.A. THE ‘GOTCHA’
Of course… there is. The mean decrease in impurity measure is a biased measure of feature importance. It favours continuous features and features with high cardinality. In 2007, Strobl et al [1] also pointed out in Bias in random forest variable importance measures: Illustrations, sources and a solution that “the variable importance measures of Breiman’s original Random Forest method … are not reliable in situations where potential predictor variables vary in their scale of measurement or their number of categories.”
Let’s try to understand why it is biased. Remember how the mean decrease in impurity is calculated? Each time a node is split on a feature, the decrease in gini index is recorded. And when a feature is continuous or has high cardinality, the feature may be split many more times than other features. This inflates the contribution of that particular feature. And what do our two culprit features have in common they are both continuous variables.
2) Drop Column Importance a.k.a Leave One Covariate Out (LOOC)
Drop Column feature importance is another intuitive way of looking at the feature importance. As the name suggests, it’s a way of iteratively removing a feature and calculating the difference in performance.
The advantages of the technique are:
 Gives a pretty accurate picture of the predictive power of each feature
 One of the most intuitive way to look at feature importance
 Model agnostic. Can be applied to any model
 The way it is calculated, it automatically takes into account all the interactions in the model. If the information in a feature is destroyed, all its interactions are also destroyed
ALGORITHM
 Use your trained model parameters and calculate the metric of your choice on an OOB sample. You can use cross validation to get the score. This is your baseline.
 Now, drop one column at a time from your training set, and retrain the model (with the same parametersand random state) and calculate the OOB score.
 Importance = OOB score – Baseline
IMPLEMENTATION
def dropcol_importances(rf, X_train, y_train, cv = 3): rf_ = clone(rf) rf_.random_state = 42 baseline = cross_val_score(rf_, X_train, y_train, scoring='accuracy', cv=cv) imp = [] for col in X_train.columns: X = X_train.drop(col, axis=1) rf_ = clone(rf) rf_.random_state = 42 oob = cross_val_score(rf_, X, y_train, scoring='accuracy', cv=cv) imp.append(baseline  oob) imp = np.array(imp) importance = pd.DataFrame( imp, index=X_train.columns) importance.columns = ["cv_{}".format(i) for i in range(cv)] return importance
Let’s do a 50 fold cross validation to estimate our OOB score. (I know it’s excessive, but let’s keep it to increase the samples for our boxplot) Like before, we are plotting the mean decrease in accuracy as well as the boxplot to understand the distribution across cross validation trials.
drop_col_imp = dropcol_importances(rf, X_train, y_train, cv=50) drop_col_importance = pd.DataFrame({'features': X_train.columns.tolist(), "drop_col_importance": drop_col_imp.mean(axis=1).values}).sort_values('drop_col_importance', ascending=False) drop_col_importance = drop_col_importance.head(25) drop_col_importance.iplot(kind='bar', y='drop_col_importance', x='features', yTitle='Drop Column Importance', xTitle='Features', title='Drop Column Importances', ) all_feat_imp_df = drop_col_imp.T order_column = all_feat_imp_df.mean(axis=0).sort_values(ascending=False).index.tolist() all_feat_imp_df[order_column[:25]].iplot(kind='box', xTitle = 'Features', yTitle='Drop Column Importance')
Click for full interactive plot
Click for full interactive plot
INTERPRETATION
 The top 4 features are still marital status, education_num, age,and hours_worked.
 fnlwgtis pushed down the list and now features after some of the onehot encoded occupations.
 randomstill occupies a high rank, positioning itself right after the hours_worked
As expected, the fnlwgt was much less important that we were led to believe from the Mean Decrease in Impurity importance. The high position of the random perplexed me a little bit and I reran the importance calculation considering all onehot features as one. i.e., dropping all the occupation columns and checking the predictive power of the occupation. When I do that, I can see random and fnlwgt rank lower than occupation, and workclass. At the risk of making the post bigger than it already is, let’s leave that investigation for another day.
So, have we got the perfect solution? The results are aligned with the Mean Decrease in Impurity, they make coherent sense, and they can be applied to any model.
JOKER IN THE PACK
The kicker here is the computation involved. To carry out this kind of importance calculation, you have to train a model multiple times, one for each feature you have and repeat that for the number of cross validation loops you want to do. Even if you have a model that trains under a minute, the time required to calculate this explodes as you have more features. To give you an idea, it took 2 hr 44 mins for me to calculate the feature importance with 36 features and 50 cross validation loops (which, of course, can be improved with parallel processing, but you get the point). And if you have a large model which is takes two days to train, then you can forget about this technique.
Another concern I have with this method is that since we are retraining the model every time with a new set of features, we are not doing a fair comparison. We remove one column and train the model again, it will find another way to derive the same information if it can, and this gets amplifies when there are collinear features. So we are mixing two things when we investigate – the predictive power of the feature and the way the model configures itself.
3) Permutation Importance
The permutation feature importance is defined to be the decrease in a model score when a single feature value is randomly shuffled [2]. This technique measures the difference in performance if you permute or shuffle a feature vector. The key idea is that a feature is important if the model performance drops if that feature is shuffled.
The advantages of this technique are:
 It is very intuitive. What is the drop in performance if the information in a feature is destroyed by shuffling it?
 Model agnostic. Even though the method was initially developed for Random Forest by Breiman, it was soon adapted to a model agnostic framework
 The way it is calculated, it automatically takes into account all the interactions in the model. If the information in a feature is destroyed, all its interactions are also destroyed
 The model need not be retrained, and hence we save on computation
ALGORITHM
 Calculate a baseline score using the metric, trained model, the feature matrix and the target vector
 For each feature in the feature matrix, make a copy of the feature matrix.
 Shuffle the feature column, pass it through the trained model to get a prediction and use the metric to calculate the performance.
 Importance = Baseline – Score
 Repeat for N times for statistical stability and take an average importance across trials
IMPLEMENTATION
The permutation importance is implemented in at least three libraries in python – ELI5, mlxtend, and in a development branch of Scikit Learn. I’ve picked the mlxtend version for purely no other reason other than convenience. According to Strobl et al. [3], “the raw [permutation] importance… has better statistical properties.” as opposed to normalizing the importance values by dividing by the standard deviation. I have checked the source code for mlxtend and Scikit Learn, and they do not normalize them.
from mlxtend.evaluate import feature_importance_permutation #This takes sometime. You can reduce this number to make the process faster num_rounds = 50 imp_vals, all_trials = feature_importance_permutation( predict_method=rf.predict, X=X_test.values, y=y_test.values, metric='accuracy', num_rounds=num_rounds, seed=1) permutation_importance = pd.DataFrame({'features': X_train.columns.tolist(), "permutation_importance": imp_vals}).sort_values('permutation_importance', ascending=False) permutation_importance = permutation_importance.head(25) permutation_importance.iplot(kind='bar', y='permutation_importance', x='features', yTitle='Permutation Importance', xTitle='Features', title='Permutation Importances', )
We also plot a box plot of all trials to get a sense of the deviation.
all_feat_imp_df = pd.DataFrame(data=np.transpose(all_trials), columns=X_train.columns, index = range(0,num_rounds)) order_column = all_feat_imp_df.mean(axis=0).sort_values(ascending=False).index.tolist() all_feat_imp_df[order_column[:25]].iplot(kind='box', xTitle = 'Features', yTitle='Permutation Importance')
Click for full interactive plot
Click for full interactive plot
INTERPRETATION
 The top 4 remains the same, but the first three(marital_status, education, age) are much more pronounced in the permutation importance
 fnlwgtand random does not even make it to the top 25
 Being an Exec Manager, or Profspeciality has a lot to do with whether you are earning >50k or not
 All in all, it resonates with our mental model of the process
Everything is hunkydory in feature importance land? Have we got the best way of explaining what features the model is using for predictions?
THE JOKER IN THE PACK
We know from life that nothing is perfect and neither is this technique. It’s Achilles’ Heel is the correlation between features. Just like drop column importance, this technique is also affected by the effect of correlation between features. Strobl et al. in Conditional variable importance for random forests [3] showed that “permutation importance overestimates the importance of correlated predictor variables.” Especially in ensemble of trees, if you have two correlated variables, some of the trees might have picker feature A and some others might have picked feature B. And while doing this analysis, in the absence of feature A, the trees which picked feature B would work well and keep the performance high and vice versa. What this will result in is that both the correlated features A and B will have inflated importance.
Another drawback of the technique is that the core idea in the technique is about permuting a feature. But that is essentially randomness that is not in our control. And because of this, the results may vary greatly. We don’t see it here, but if the box plot shows a large variation in importance for a feature across trials, I’ll be wary in my interpretation.
Correlation Coefficient [7] (inbuilt in pandas profiling which considers categorical variables as well).
There is yet another drawback to this technique, which in my opinion, is the most concerning. Giles Hooker et al. [6] says, “When features in the training set exhibit statistical dependence, permutation methods can be highly misleading when applied to the original model.”
Let’s consider occupation and education. We can understand this from two perspectives:
 Logical: If you think about it, occupationand education have a definite dependence. You can only get a few jobs if you have sufficient education and statistically, you can draw parallels between them. So if we are permuting any one of those columns, it would create feature combinations which do not make sense. A person with education as 10th and occupation as Profspeciality doesn’t make a lot of sense, does it? So, when we are evaluating the model, we are evaluating nonsensical cases like these, which muddles up the metric which we use to assess the feature importance.
 Mathematical: occupationand education have strong statistical dependence(we can see that from the correlation plot above). So, when we are permuting any one of these features, we are forcing the model to explore unseen subspaces in the highdimensional feature space. And this forces the model to extrapolate and this extrapolation is a significant source of error.
Giles Hooker et al. [6] suggests alternative methodologies that combine LOOC and Permutation methods, but all the alternatives are computationally more intensive and do not have a strong theoretical guarantee of having better statistical properties.
DEALING WITH CORRELATED FEATURES
After identifying the highly correlated features, there are two ways of dealing with correlated features.
 Group the highly correlated variables together and evaluate only one feature from the group as a representative of the group
 When you permute the columns, permute the whole group of features in one trial.
N.B. The second method is the same method that I would suggest to deal with onehot variables.
SIDENOTE (TRAIN OR VALIDATION)
During the discussion of both Drop Column importance and Permutation importance, one question should have come to your mind. We passed the test/validation set to the methods to calculate the importance. Why not train set?
This is a grey area in the application of some of these methods. There is no right or wrong here because there are arguments for and against both. In Interpretable Machine Learning, Christoph Molnar argues a case for both train and validation sets.
The case for test/validation data is a nobrainer. For the same reason why we do not judge a model by the error in the training set, we cannot judge the feature importance on the performance on the training set (especially since the importance is intrinsically linked to the error).
The case for train data is counterintuitive. But if you think about it, you’ll see that what we want to measure is how the model is using the features. And what better data to judge that than the training set on which the model was trained? Another trivial issue is also that we would ideally train the model on all available data, and in such an ideal scenario, there will not be a test or validation data to check performances on. In Interpretable Machine Learning [5], section 5.5.2 discusses this issue at length and even with a synthetic example of an overfitting SVM.
It all comes down to whether you want to know what features the model relies on to make predictions or the predictive power of each feature on unseen data. For e.g., if you are evaluating feature importance in the context of feature selection, do not use test data in any circumstances (there you are overfitting your feature selection to fit the test data).
4) Partial Dependence Plots (PDP) and Individual Conditional Expectation (ICE) plots
All the techniques we reviewed until now looked at the relative importance of different features. Now let’s move slightly in a different direction and look at a few techniques which explore how a particular feature interacts with the target variable.
Partial Dependence Plots and Individual Conditional Expectation plots help us understand the functional relationship between the features and the target. They are graphical visualizations of the marginal effect of a given variable(or multiple variables) on an outcome. Friedman(2001) introduced this technique in his seminal paper Greedy function approximation: A gradient boosting machine[8].
Partial Dependence Plots show an average effect, whereas ICE plots show the functional relationship for individual observations. PD plots show the average effect whereas ICE plots show the dispersion or heterogeneity of the effect.
The advantages of this technique are:
 The calculation is very intuitive and is easy to explain in layman terms
 We can understand the relationship between a feature or a combination of features with the target variable, i.e. if it is linear, monotonic, exponential etc.
 They are easy to compute and implement
 They give a causal interpretation, as opposed to a feature importance style interpretation. But what we have to keep in mind is that the causal interpretation of how the model sees the world and now the real world.
ALGORITHM
Let’s consider a simple situation where we plot the PD plot for a single feature x, with unique values . The PD plot can be constructed as follows:
 For
 Copy the training data and replace the original values of x with
 Use the trained model to generate predictions for the modified copy of the entire training data
 Store all the predictions against in a maplike data structure
 For PD plot:
 Calculate the average predictions for each
 Plot the pairs
 For ICE plot:
 Plot all the pairs . N is the total number of observations in the training set.
 In practice, instead of taking all the possible values of a feature, we define a grid of intervals for the continuous variable to save on computation.
 For a categorical variable also, this definition holds, But we won’t be defining a grid there. Instead we take all the unique values in the category(or all the onehot encoded variables pertaining to a categorical feature) and calculate the ICE and PD plots using the same methodology.
 If the process is still unclear to you, I suggest looking at this medium post (by the author of PDPbox, a python library for plotting PD plots.
IMPLEMENTATION
I have found the PD plots implemented in PDPbox, skater and Scikit Learn. And the ICE plots in PDPbox, pyCEbox, and skater. Out of all of these, I found PDPbox to be the most polished. And they also support 2 variable PDP plots as well.
from pdpbox import pdp, info_plots pdp_age = pdp.pdp_isolate( model=rf, dataset=X_train, model_features=X_train.columns, feature='age' ) #PDP Plot fig, axes = pdp.pdp_plot(pdp_age, 'Age', plot_lines=False, center=False, frac_to_plot=0.5, plot_pts_dist=True,x_quantile=True, show_percentile=True) #ICE Plot fig, axes = pdp.pdp_plot(pdp_age, 'Age', plot_lines=True, center=False, frac_to_plot=0.5, plot_pts_dist=True,x_quantile=True, show_percentile=True)
ICE plot for Age.
Let me take some time to explain the plot. On the xaxis, you can find the values of the feature you are trying to understand, i.e., age. On the yaxis, you find the prediction. In the case of classification, it is the prediction probability and in the case of regression, it is the realvalued prediction. The bar on the bottom represents the distribution of training data points in different quantiles. It helps us gauge the goodness of the inference. The parts where the number of points is very less, the model has seen fewer examples and the interpretation can be tricky. The single line in the PD plot shows the average functional relationship between the feature and target. All the lines in the ICE plot show the heterogeneity in the training data, i.e., how all the observations in the training data vary with the different values of the age.
INTERPRETATION
 agehas a largely monotonic relationship with the earning capacity of a person. The older a person is, the more likely he is to earn above 50k
 The ICE plots show a lot of dispersion. But all of it shows the same kind of behaviour that we see in the PD plot
 The training observations are considerable well balanced across the different quantiles.
Now, let’s also take an example with a categorical feature, like education. PDPbox has a very nice feature where it lets you pass a list of features as an input and it will calculate the PDP for them considering them as categorical features.
# All the onehot variables for the occupation feature occupation_features = ['occupation_ ?', 'occupation_ Admclerical', 'occupation_ ArmedForces', 'occupation_ Craftrepair', 'occupation_ Execmanagerial', 'occupation_ Farmingfishing', 'occupation_ Handlerscleaners', 'occupation_ Machineopinspct', 'occupation_ Otherservice', 'occupation_ Privhouseserv', 'occupation_ Profspecialty', 'occupation_ Protectiveserv', 'occupation_ Sales', 'occupation_ Techsupport', 'occupation_ Transportmoving'] #Notice we are passing the list of features as a list with the feature parameter pdp_occupation = pdp.pdp_isolate( model=rf, dataset=X_train, model_features=X_train.columns, feature=occupation_features ) #PDP fig, axes = pdp.pdp_plot(pdp_occupation, 'Occupation', center = False, plot_pts_dist=True) #Processing the plot for aesthetics _ = axes['pdp_ax']['_pdp_ax'].set_xticklabels([col.replace("occupation_","") for col in occupation_features]) axes['pdp_ax']['_pdp_ax'].tick_params(axis='x', rotation=45) bounds = axes['pdp_ax']['_count_ax'].get_position().bounds axes['pdp_ax']['_count_ax'].set_position([bounds[0], 0, bounds[2], bounds[3]]) _ = axes['pdp_ax']['_count_ax'].set_xticklabels([])
INTERPRETATION
 Most of the occupations have very minimal effect on how much you earn.
 The ones that stand out from the rest are, Execmanagerial, Profspeciality, and Tech Support
 But, from the distribution, we know that there were very little training examples for Techsupport, and hence we take that with a grain of salt.
INTERACTION BETWEEN MULTIPLE FEATURES
PD plots can be theoretically drawn for any number of features to show their interaction effect, as well. But practically, we can only do it for two, at the max three. Let’s take a look at an interaction plot between two continuous features age and education(education and age are not truly continuous, but for lack of better example we are choosing them).
There are two ways you can plot a PD plot between two features. There are three dimensions here, feature value 1, feature value 2, and the target prediction. Either, we can plot a 3D plot or a 2D plot with the 3rd dimension depicted as color. I prefer the 2D plot because I think it conveys the information in a much more crisp manner than a 3D plot where you have to look at the 3D shape to infer the relationship. PDPbox implements the 2D interaction plots, both as a contour plot and grid. Contour works best for continuous features and grid for categorical features.
# Age and Education inter1 = pdp.pdp_interact( model=rf, dataset=X_train, model_features=X_train.columns, features=['age', 'education_num'] ) fig, axes = pdp.pdp_interact_plot( pdp_interact_out=inter1, feature_names=['age', 'education_num'], plot_type='contour', x_quantile=False, plot_pdp=False ) axes['pdp_inter_ax'].set_yticklabels([edu_map.get(col) for col in axes['pdp_inter_ax'].get_yticks()])
INTERPRETATION
 Even though we observed a monotonic relationship with age when looked at isolation, now we know that it is not universal. For e.g., look at the contour line to the right of the 12theducation level. It’s pretty flat as compared to the lines for somecollege and above. What it really shows that your probability of getting more than 50k doesn’t only increase with age, but it also has a bearing on your education. A college degree ensures you increase your earning potential as you get older.
This is also a very useful technique to investigate bias(the ethical kind) in your algorithms. Suppose if we want to look at the algorithmic bias in the sex dimension.
#PDP Sex pdp_sex = pdp.pdp_isolate( model=rf, dataset=X_train, model_features=X_train.columns, feature='sex' ) fig, axes = pdp.pdp_plot(pdp_sex, 'Sex', center=False, plot_pts_dist=True) _ = axes['pdp_ax']['_pdp_ax'].set_xticklabels(sex_le.inverse_transform(axes['pdp_ax']['_pdp_ax'].get_xticks())) # marital_status and sex inter1 = pdp.pdp_interact( model=rf, dataset=X_train, model_features=X_train.columns, features=['marital_status', 'sex'] ) fig, axes = pdp.pdp_interact_plot( pdp_interact_out=inter1, feature_names=['marital_status', 'sex'], plot_type='grid', x_quantile=False, plot_pdp=False ) axes['pdp_inter_ax'].set_xticklabels(marital_le.inverse_transform(axes['pdp_inter_ax'].get_xticks())) axes['pdp_inter_ax'].set_yticklabels(sex_le.inverse_transform(axes['pdp_inter_ax'].get_yticks()))
PD Plot of sex on the left and PD interaction plot of sex and marital_status on the right.
 If we look at just the PD plot of sex, we would conclude that there is no real discrimination based on the sex of the person.
 But, just take a look at the interaction plot with marital_status. On the lefthand side(married), both the squares have the same color and value, but on the righthand side(single) there is a difference between Female and Male
 We can conclude that being a single male gives you a much better chance of getting more than 50k than being a single female. (Although I wouldn’t start an allout war against sexual discrimination based on this, it would definitely be a starting point in the investigation.
JOKER IN THE PACK
The assumption of independence between the features is the biggest flaw in this approach. The same flaw that is present in LOOC importance and Permutation Importance is applicable to PDP and ICE plots. Accumulated Local Effects plots are a solution to this problem. ALE plots solve this problem by calculating – also based on the conditional distribution of the features – differences in predictions instead of averages.
To summarize how each type of plot (PDP, ALE) calculates the effect of a feature at a certain grid value v:
Partial Dependence Plots: “Let me show you what the model predicts on average when each data instance has the value v for that feature. I ignore whether the value v makes sense for all data instances.”
ALE plots: “Let me show you how the model predictions change in a small ‘window’ of the feature around v for data instances in that window.”
In the python environment, there is no good and stable library for ALE. I’ve only found one ALEpython, which is still very much in development. I managed to get an ALE plot of age, which you can find below. But got an error when I tried an interaction plot. It’s also not developed for categorical features.
ALE plot for age.
This is where we break off again and push the rest of the stuff to the next blog post. In the next part, we take a look at LIME, SHAP, Anchors, and more.
Full Code is available in my Github.
Original. Reposted with permission.
Bio: Manu Joseph (@manujosephv) is an inherently curious and selftaught Data Scientist with about 8+ years of professional experience working with Fortune 500 companies, including a researcher at Thoucentric Analytics.
Related:
Top Stories Past 30 Days

