Assignment 8: Explainable Machine Learing (Part 1)

Assignment 8: Explainable Machine Learing (Part 1)
Objective:
Interpreting how a machine learning model works is crucial from many aspects, such as debugging the model, improving the model, and derive new insights from the model. There are many model explanation techniques that we can use for interpreting a black box model itself or its predictions. Also, some models are transparent that we can directly get insights from its structure or learned parameters. In this assignment, you are going to train different models on a dataset and try different approaches to explain the model and get some insights. After completing this assignment, you should be able to answer the following questions:

How to explain transparent models?
How to implement the permutation method to explain black-box models?
How to create and interpret partial dependence plots?
How to implement the global surrogate method to explain black-box models?
How to use SHAP to explain model predictions and interpret its explanation results?
How to use LIME to explain model predictions and interpret its explanation results?
The data can be downloaded from A8-1-data.zip

0. Preparation
Import relevant libraries and load the dataset:

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from dataprep.eda import plot

from sklearn.model_selection import train_test_split
from sklearn.inspection import PartialDependenceDisplay
from sklearn.metrics import log_loss
from sklearn.preprocessing import MinMaxScaler
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.neural_network import MLPClassifier
The dataset we are going to use is a binary forest covertype dataset and we will predict the forest cover type from cartographic variables only. The first ten features are numerical features, there are also two categorical features in a one-hot encoding fashion (4 and 40 vector length each), the last column is the target forest cover type (binary). More details about the dataset can be found at https://archive.ics.uci.edu/ml/datasets/covertype. Please make sure you understand the meaning of each feature.

# load data and take a look at data distribution
data = pd.read_csv('bforest_sample.csv', delimiter=',')
plot(data)
Stats and Insights
Number of plots per page:  
9
 
Elevation
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Aspect
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Slope
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Horizontal_Distance_To_Hydrology
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Vertical_Distance_To_Hydrology
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Column Insights
Vertical_Distance_To_Hydrology is skewed
Vertical_Distance_To_Hydrology has 3242 (6.55%) zeros
Vertical_Distance_To_Hydrology has 4895 (9.89%) negatives
Horizontal_Distance_To_Roadways
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Hillshade_9am
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Hillshade_Noon
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Hillshade_3pm
'hist.bins': 50
Number of bins in the histogram
'hist.yscale': 'linear'
Y-axis scale ("linear" or "log")
'hist.color': '#aec7e8'
Color
Copy All Parameters
Prev1234567Next
# split data into training and test dataset with respect to ratio 0.8:0.2
train, test = train_test_split(data, test_size=0.2, random_state=733)
X_train, y_train = train.iloc[:, :-1], train.iloc[:, -1]
X_test, y_test = test.iloc[:, :-1], test.iloc[:, -1]
feature_names = list(X_train)

# rescale
scaler = MinMaxScaler()
X_train = pd.DataFrame(data=scaler.fit_transform(X_train), columns=feature_names)
X_test = pd.DataFrame(data=scaler.transform(X_test), columns=feature_names)
Task 1. Transparent Model
In this task, you are going to train a logistic regression model and interpret it. For logistic regression, since the effect of each feature is simply added together, we can interpret it directly by looking at the coefficient of each feature. Please follow the comment to finish the code:

# train a logistic regression model
lr = LogisticRegression(solver='liblinear').fit(X_train, y_train)
# set default figure size
plt.rcParams['figure.figsize'] = [16, 8]
# show the coefficient value for each feature by a bar chart
def explain_logistic_regression(lr, feature_names):
    x = np.arange(len(feature_names))
    y = lr.coef_[0]
    plt.bar(x,y)
    plt.xticks(x, feature_names, rotation = "vertical")
    plt.grid(True)
    plt.xlabel("feature name")
    plt.ylabel("coefficient")
    plt.title("Logistic Regression Coefficient Value")
    plt.show()
explain_logistic_regression(lr, feature_names)

For a specific prediction, we can get a more concrete effect of each feature by the product of the coefficient and input feature value. Please follow the comment to finish the code:

# show the effect from each input feature by a bar chart
def explain_logistic_regression_prediction(lr, feature_names, sample):
   x = np.arange(len(feature_names))
   y = sample[0] * lr.coef_[0]
   plt.bar(x, y)
   plt.xticks(x, feature_names, rotation = "vertical")
   plt.grid(True)
   plt.xlabel("feature name")
   plt.ylabel("coefficient * feature")
   plt.title("Logistic Regression Coefficient * sample feature")
   plt.show()
explain_logistic_regression_prediction(lr, feature_names, (X_test.iloc[0, :], y_test.iloc[0]))

What can you get from the above plots? Please write down two findings:

Findings:

The Elevation feature contributes the most positive effect on the outcome, while the Hillshade_Noon gives the most negative influence to the result.
Sample 0 belongs to soil type 29, which gives a negative infulunce to the result.
Task 2. Post-hoc Explanation (Global Model)
In this task, you are going to build a gradient boosting tree model and a neural network, and use some techniques we introduced in class to interpret these models. First, let's train the models using training data:

# train a gradient boosting classifier
gb = GradientBoostingClassifier().fit(X_train.values, y_train.values)
# train a multi-layer perceptron classifier
mlp = MLPClassifier(learning_rate_init=0.1).fit(X_train.values, y_train.values)
2.1 Permutation
# get importance score (E' - E) for each feature by permutation, use log_loss as error
def permutation_importance(model, feature_names, X, y):
    #get log error for non-permutated features
    pred_y = model.predict(X)
    E1 = log_loss(y, pred_y)
    
    scores = [0 for i in range(len(feature_names))]

    #get permutation score for each feature
    for index, feature in enumerate(feature_names):
        # copy array
        tmp = np.copy(X)
        #shuffle current feature values
        np.random.shuffle(tmp[:,index])
        #predict permutated label
        pred_y_permu = model.predict(tmp)
        #get importance score
        E2 = log_loss(y, pred_y_permu)
        scores[index] = E2-E1

    # show the top 5 most important features
    print("-------------Top 5 most important features-------------")
    scores_index = np.argsort(np.array(scores))[::-1]
    scores.sort(reverse=True)
    for i in range(5):
        print(feature_names[scores_index[i]], scores[i])
    
permutation_importance(gb, feature_names, X_test.to_numpy(), y_test.to_numpy())
permutation_importance(mlp, feature_names, X_test.to_numpy(), y_test.to_numpy())
-------------Top 5 most important features-------------
Elevation 8.444034723089143
Horizontal_Distance_To_Roadways 0.7315736979917755
Horizontal_Distance_To_Hydrology 0.7060960070169386
Hillshade_Noon 0.6260232639531607
Soil_Type_22 0.32393064239436864
-------------Top 5 most important features-------------
Elevation 6.711551736800165
Wilderness_Area_1 1.5941755209970045
Wilderness_Area_3 1.1319374133106583
Hillshade_Noon 0.8808001737015418
Soil_Type_32 0.7970877604985018
2.2 Partial dependence plots
In this section, you are going to use the plot_partial_dependence provided by sklearn to see the marginal effect of each single numerical features:

# plot partial dependence for numerical features (first 10 features) for gradient boosting classifier
plt.rcParams['figure.figsize'] = [16, 10]
PartialDependenceDisplay.from_estimator(gb, X_test, features = np.arange(0,10) ,feature_names = feature_names[0:10],n_cols=4)
plt.show()

# plot partial dependence for numerical features (first 10 features) for neural network
PartialDependenceDisplay.from_estimator(mlp, X_test.values, features = np.arange(0,10) ,feature_names = feature_names[0:10],n_cols=4)
<sklearn.inspection._plot.partial_dependence.PartialDependenceDisplay at 0x7f91cf5651f0>

2.3 Global Surrogate
Now, let's train a simple logistic regression based on the gradient boosting tree model and neural network we built before. And use the method we implemented previously to interpret the derived transparent model:

# train a logistic regression model on gb and mlp and explain using explain_logistic_regression()
gb_pred_y = gb.predict(X_train.values)
mlp_pred_y = mlp.predict(X_train.values)

lr_gb = LogisticRegression(max_iter=500)
lr_gb.fit(X_train, gb_pred_y)

lr_mlp = LogisticRegression(max_iter=500)
lr_mlp.fit(X_train, mlp_pred_y)

explain_logistic_regression(lr_gb, feature_names)
explain_logistic_regression(lr_mlp, feature_names)


Task 3. Post-hoc Explanation (Single Prediction)
3.1 Attribution
From now, let's focus on interpreting single predictions. SHAP is an efficient method to approximatly calculate the shapely value we mentioned in class. Please install the library and take a look at the doc. Explain the first prediction made by the gradient boosting tree model in the test set:

!pip install shap
import shap

# load JS visualization code to notebook
shap.initjs()

# use Tree SHAP explainer to explain the gradient boosting tree model
# you only need to explain and plot the first explaination
explainer = shap.Explainer(gb)
shap_values = explainer(X_test)

# visualize the first prediction's explanation
shap.plots.force(shap_values[0])

What can you get from the above representation? Please write down two findings:

Findings:

The Elevation feature has the biggest influence on the prediction as we have seen from Task 1.
The SHAP base value for this dataset is -0.4003, which indicates that the average or reference prediction value for the training data set is negative.
3.2 LIME
LIME is a library implemented by the authors of the paper. Please install the library and take a look at the doc and tutorial. Use it to explain the first prediction made by the neural network model in the test set:

!pip install lime
import lime
import lime.lime_tabular

# use LimeTabularExplainer to explain the neural network model
# you only need to explain and plot the first explaination
explainer = lime.lime_tabular.LimeTabularExplainer(X_test.values, feature_names=feature_names, kernel_width=5, verbose=False)
exp = explainer.explain_instance(X_test.values[0,:],mlp.predict_proba, num_features=54)
exp.show_in_notebook()
Prediction probabilities
0.41
0
0.59
1
0
1
Elevation > 0.66
0.40
Soil_Type_34 <= 0.00
0.36
Soil_Type_03 <= 0.00
0.36
Soil_Type_18 <= 0.00
0.26
Soil_Type_02 <= 0.00
0.25
Soil_Type_38 <= 0.00
0.25
Soil_Type_04 <= 0.00
0.24
Soil_Type_32 <= 0.00
0.22
Soil_Type_39 <= 0.00
0.21
Soil_Type_28 <= 0.00
0.18
Soil_Type_22 <= 0.00
0.18
Soil_Type_10 <= 0.00
0.16
Soil_Type_07 <= 0.00
0.15
Soil_Type_35 <= 0.00
0.14
Soil_Type_21 <= 0.00
0.14
Soil_Type_40 <= 0.00
0.12
Soil_Type_24 <= 0.00
0.12
Soil_Type_17 <= 0.00
0.11
Horizontal_Distance_...
0.11
Soil_Type_20 <= 0.00
0.10
Soil_Type_13 <= 0.00
0.09
Soil_Type_25 <= 0.00
0.09
0.00 < Wilderness_Are...
0.08
Soil_Type_26 <= 0.00
0.08
Soil_Type_23 <= 0.00
0.07
Soil_Type_19 <= 0.00
0.07
Soil_Type_33 <= 0.00
0.06
Soil_Type_09 <= 0.00
0.06
Soil_Type_08 <= 0.00
0.06
Wilderness_Area_4 <=...
0.06
Slope <= 0.13
0.06
Wilderness_Area_3 <=...
0.05
Soil_Type_27 <= 0.00
0.05
Soil_Type_36 <= 0.00
0.05
Soil_Type_12 <= 0.00
0.05
Soil_Type_16 <= 0.00
0.04
Wilderness_Area_2 <=...
0.04
0.83 < Hillshade_Noo...
0.04
Vertical_Distance_To...
0.04
Horizontal_Distance_T...
0.03
Soil_Type_06 <= 0.00
0.03
Horizontal_Distance_T...
0.02
Soil_Type_31 <= 0.00
0.02
0.85 < Hillshade_9am ...
0.02
0.16 < Aspect <= 0.36
0.01
Soil_Type_29 > 0.00
0.01
Soil_Type_30 <= 0.00
0.01
0.48 < Hillshade_3pm ...
0.01
Soil_Type_11 <= 0.00
0.00
Soil_Type_14 <= 0.00
0.00
Soil_Type_37 <= 0.00
0.00
Soil_Type_01 <= 0.00
0.00
Soil_Type_05 <= 0.00
0.00
Soil_Type_15 <= 0.00
0.00
Feature    Value
Elevation    0.74
Soil_Type_34    0.00
Soil_Type_03    0.00
Soil_Type_18    0.00
Soil_Type_02    0.00
Soil_Type_38    0.00
Soil_Type_04    0.00
Soil_Type_32    0.00
Soil_Type_39    0.00
Soil_Type_28    0.00
Soil_Type_22    0.00
Soil_Type_10    0.00
Soil_Type_07    0.00
Soil_Type_35    0.00
Soil_Type_21    0.00
Soil_Type_40    0.00
Soil_Type_24    0.00
Soil_Type_17    0.00
Horizontal_Distance_To_Hydrology    0.34
Soil_Type_20    0.00
Soil_Type_13    0.00
Soil_Type_25    0.00
Wilderness_Area_1    1.00
Soil_Type_26    0.00
Soil_Type_23    0.00
Soil_Type_19    0.00
Soil_Type_33    0.00
Soil_Type_09    0.00
Soil_Type_08    0.00
Wilderness_Area_4    0.00
Slope    0.10
Wilderness_Area_3    0.00
Soil_Type_27    0.00
Soil_Type_36    0.00
Soil_Type_12    0.00
Soil_Type_16    0.00
Wilderness_Area_2    0.00
Hillshade_Noon    0.88
Vertical_Distance_To_Hydrology    0.36
Horizontal_Distance_To_Fire_Points    0.16
Soil_Type_06    0.00
Horizontal_Distance_To_Roadways    0.74
Soil_Type_31    0.00
Hillshade_9am    0.91
Aspect    0.33
Soil_Type_29    1.00
Soil_Type_30    0.00
Hillshade_3pm    0.54
Soil_Type_11    0.00
Soil_Type_14    0.00
Soil_Type_37    0.00
Soil_Type_01    0.00
Soil_Type_05    0.00
Soil_Type_15    0.00
Run LIME multiple times, (Q1) what do you think of the stability of LIME? (Q2) Can you briefly explain the reason?

Your Answer:

The outcome of LIME is not the same each time. To be more specific, the weight of each feature changed.
The instability of LIME is due to the random sampling process from the original dataset and perturbations performed on the input features. This random nature leads to different outcomes when repeat the explanation on the same dataset with the same model.
Submission
Complete the code in this notebook, and submit it to the CourSys activity Assignment 8.