Lab 8 - Model Explainability

This lab covers model agnostic approaches to global and local explainability.

We'll use the lime and shap libraries in addition to our standard set. You'll likely need to install each using pip:

# pip install lime
import lime

# pip install shap
import shap

## Standard Libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn
import warnings   # Turns off the future warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
import xgboost as xgb

The lab will use the Ames Housing data set, which is becoming slightly dated (it contains homes sold between 2006 and 2010) but remains one of the most comprehensive publicly available housing data sets.

## Read the data
ames = pd.read_csv("https://remiller1450.github.io/data/AmesHousing.csv")

## Separate y and X, and keep only a few of the more interesting features to make things manageable
ames_y = ames['SalePrice']
ames_X = ames.drop('SalePrice', axis=1)[['Overall.Qual', 'Overall.Cond', 'Lot.Area', 'Year.Built', 'Total.Bsmt.SF',
                                         'Full.Bath', 'Half.Bath', 'Bedroom.AbvGr', 'Fireplaces', 'Gr.Liv.Area',
                                        'Garage.Area', 'Wood.Deck.SF']]

## Train-test split
from sklearn.model_selection import train_test_split
train_X, test_X, train_y, test_y = train_test_split(ames_X, ames_y, test_size=0.2, random_state=7)

## Create a separate validation set to save the test set for later
train_X, val_X, train_y, val_y = train_test_split(train_X, train_y, test_size=0.25, random_state=7)

Part 1 - Global Approaches

This section will introduce key functions and a common workflow for assessing and understanding the global importance of variables in a machine learning model. The basic steps in this process are:

1. Develop a strong predictive model (using the methods described in our previous labs)
2. Calculate global importance metrics for each feature (permutation importance).
3. Use partial dependence plots (PDPs) to understand the roles of features identified as important.
4. Use individual conditional expectation plots to assess whether any features exhibit heterogeneous effects.

To demonstrate these steps, we'll first train a random forest model:

## Fit random forest
from sklearn.ensemble import RandomForestRegressor
rf_mod = RandomForestRegressor(n_estimators=500, max_depth=6, random_state=7).fit(train_X, train_y)

Next, the permutation_importance() function from sklearn is used to find permutation importance scores.

We'll perform 20 repetitions for each variable, though in practice you might consider doing more to promote stability (depending on the size of your data and any computational constraints you might have):

## Permutation importance
from sklearn.inspection import permutation_importance
rf_imp = permutation_importance(rf_mod, val_X, val_y, n_repeats=20, scoring="neg_root_mean_squared_error", random_state=7)

## Assemble and print results
imp_results = pd.DataFrame({"Feature": val_X.columns,
                       "Importance": rf_imp.importances_mean,
                       "Std": rf_imp.importances_std})
imp_results.sort_values("Importance", ascending=False)

	Feature	Importance	Std
0	Overall.Qual	32612.282287	1824.206044
9	Gr.Liv.Area	12894.538614	1045.266957
4	Total.Bsmt.SF	6605.847408	649.986663
3	Year.Built	1695.798056	312.379495
2	Lot.Area	1424.600545	157.522896
10	Garage.Area	1020.410589	221.871921
5	Full.Bath	661.940552	123.076584
8	Fireplaces	642.781767	130.901568
1	Overall.Cond	302.105425	131.889216
7	Bedroom.AbvGr	175.250275	39.985661
6	Half.Bath	39.889757	53.468468
11	Wood.Deck.SF	39.693638	26.591396

Each importance score measures the deterioration in the evaluation metric (RMSE in this example) when the predictive information in a feature is destroyed by permuting its values. For example, RMSE increases by $32612.28 when the most important feature, Overall.Qual, is permuted. Overall, our results suggest features like Overall.Qual and Gr.Liv.Area are much bigger contributors to the model's performance than features like Half.Bath and Wood.Deck.SF.

As a side note: There is some debate about whether feature importance should be calculated using held out data (validation set) or the training data. In our example we use the validation data, which is generally a better practice, as it avoids overestimating the importance of variables that the model may have "memorized" during training.

As our next step, we might look at the PDPs of some of the more influential predictors:

from sklearn.inspection import PartialDependenceDisplay
import matplotlib.pyplot as plt
PartialDependenceDisplay.from_estimator(rf_mod, val_X, features=["Overall.Qual", "Gr.Liv.Area", "Year.Built"])
plt.show()

When supplied together as a list, PartialDependenceDisplay() will use a common Y-axis for all features, which allows us to draw a nice connection with the permutation importance (greater changes for Overall.Qual than for Year.Built), but it can also hide the form of less influential predictors.

For example, our initial plot seemingly showed a relatively linear effect of Year.Built, but upon closer inspection there is actually more of a piecewise relationship where newer build years aren't associated with any change in the average predicted home price until around 1920.

PartialDependenceDisplay.from_estimator(rf_mod, val_X, features=["Year.Built"])
plt.show()

We can construct ICE plots using the PartialDependenceDisplay() with the argument kind = "both":

PartialDependenceDisplay.from_estimator(rf_mod, val_X,features=["Overall.Qual", "Gr.Liv.Area"], kind="both")
plt.show()

In this example we can judge that there is relatively little effect heterogeneity in how Overall.Qual contributes to sale price, most of the lines are roughly parallel to the average trend. But there is some evidence of effect heterogeneity in Gr.Liv.Area, as we see the lines for some observations cross (particularly noticeable around 1500, 2000, and 2300 sqft).

Question #1:

• Part A: Follow the workflow from Lab 8 to build an effective xgboost model that predicts sale price using the Ames Housing data. Report the best estimator and its cross-validated performance.
• Part B: An important starting point is to ask whether a model is accurate enough to actually warrant explaining. With this in mind, find the cross-validated performance of the random forest model used in this sections examples and compare it with the results you found in Part A.
• Part C: Compute the permutation importance of each feature in the validation set using the xgboost model from Part A. Do both the random forest and xgboost models agree on the most important features? Do any features appear substantially more important in one of these models?
• Part D: Create PDPs for Overall.Qual and Gr.Liv.Area. Briefly comment on how the shapes of these plots compare with those for the random forest model.
• Part E: Now create ICE plots for Overall.Qual and Gr.Liv.Area. Which feature seems to show more heterogeneity? How does the amount of heterogeneity compare with the random forest model?

Part 2 - Local Approaches

While PDPs and ICE plots help us understand the overall behavior of certain features within the fitted model, it can often be useful to understand how the model arrived at a specific prediction. For example, a prospective home seller might want to understand why a model predicts a certain sale price for the specific home they are trying to sell.

The example below shows how to use Local Interpretable Model-Agnostic Explanations (LIME) to explain the predicted sale price for the first observation in the validation set:

from lime.lime_tabular import LimeTabularExplainer
warnings.filterwarnings("ignore", message="X does not have valid feature names")

## First, fit the explainer
lime_explainer = LimeTabularExplainer(training_data=train_X.values,feature_names=train_X.columns,mode="regression",random_state=15)

## Next, provide the observation whose prediction is to be explained
lime_exp = lime_explainer.explain_instance(val_X.iloc[0].values, rf_mod.predict,num_features=10)

## Plot the local explanation
lime_exp.as_pyplot_figure()
plt.show()

Noting that the predicted sale price for this home is $$$162,822.97, a simple interpretation of the LIME results is that this home is a low quality, older home that is being rewarded for its larger living area, garage, and deck.

While LIME provides an approximate explanation using local perturbations, it does not directly attempt to allocate responsibility for the prediction among the input features. Additionally, LIME results can vary depending on the random_state used (compare random_state=7 and random_state=15 in the example). SHAP takes a different approach by attributing the prediction directly to the input features, providing a decomposition of the prediction into a baseline value plus feature-specific contributions.

For these reasons, SHAP is now much more widely used in practice than LIME. Additionally, the SHAP library is such that the typical workflow will apply the explainer to all of the validation set, rather than one observation at a time. We can then use these values for both local and global explanations:

## Fit the explainer
shap_explainer = shap.TreeExplainer(rf_mod)

## Apply the explainer
shap_values = shap_explainer.shap_values(val_X)

## Waterfall Plot - a local explanation of first home in the validation set
home1_shap = shap.Explanation(values=shap_values[0], base_values= shap_explainer.expected_value, data=val_X.iloc[0], feature_names=val_X.columns)
shap.plots.waterfall(home1_shap)

## Beeswarm Plot - a global explanation of each feature's role in the model
shap.summary_plot(shap_values, val_X)

The waterfall plot explains a single prediction by showing how each feature moved its prediction from the baseline of $\$$181,341.00 to reach its model prediction of $$$162,822.97.

The second plot, known as a beeswarm plot, shows the SHAP values for each variable for each individual observation in the validation set. The color of the point reflects where an observation's feature value falls in the observed data range, and the X-axis position shows the SHAP value. This plot simultaneously tells us:

1. Importance - features with larger average absolute SHAP values across all observations appear near the top
2. Direction - because positive SHAP values increase predictions, seeing red colored observations (high feature values) tending to appear on the right side of the plot indicates a positive relationship. Similarly, seeing blue colored observations (low feature values) on the right side of the plot would indicate a negative relationship.
3. Heterogeneity - features whose SHAP values exhibit greater spread indicate the effect of that feature tends to vary more across observations, generally because of interactions with other predictors.

Additionally, the beeswarm plot can help us judge:

1. Monotonicity - A clear blue-to-red gradient from low to high SHAP values suggests that larger feature values tend to consistently increase (or decrease) predictions. An approximate linear relationship would show SHAP values that change proportionally over a feature's range, but this is difficult to judge from the beeswarm plot alone.
2. Interactions - If observations with similar feature values (similar colors) receive very different SHAP values, the model is likely using that feature differently depending upon the values of other predictors.

Question #2: For this question you should use the xgboost model you fit as part of Question #1.

• Part A: Create a LIME explanation plot for the 5th observation in the validation set. Provide a short, 1-2 sentence description of how you would summarize why this home received the predicted value it did in the xgboost model.
• Part B: Create a SHAP waterfall plot for the same home (the 5th observation in the validation set). Does this plot generally tell the same story as the LIME plot? Briefly explain any similarities or differences.
• Part C: Create a SHAP beeswarm plot and identify:
  • i - the most important and least important feature, according to mean absolute SHAP value.
  • ii - a feature that might not have a monotonic relationship with the outcome in the model.
• Part D: Compare the ranking of features in the beeswarm SHAP plot from Part C with the permutation importances you found in Question #1. Do they generally agree with each other? In the case of disagreements, which do you prefer? Briefly explain your preference, noting that there isn't necessarily a correct answer.
• Part E: A SHAP dependence plot is similar to a PDP, but instead of displaying the average effect of a feature, it displays the SHAP value for every observation. The plot's x-axis shows the feature value and the y-axis shows the corresponding SHAP value, with larger vertical spread indicating a feature's effect varies across observations, typically due to interactions with other predictors. For this question, you should consult the SHAP library documentation and create a SHAP dependence plot for the feature you identified in Part C - ii.