How to calculate shapley values from scratch

The shapley value is a popular and useful method to explain machine learning models. The shapley value of a feature is the average contribution of a feature value to the prediction. In this article I’ll show you how to compute shapley values from scratch. This might help your understanding of the method.

Introduction

The exact calculation of shapley values becomes problematic if more features are added, since the computation time grows exponentially. So Strumbelj et al.(2014) propose an approximation with Monte-Carlo sampling to calculate the shapley value for the $j$-th feature:

$$\phi_{j}=\frac{1}{M}\sum_{m=1}^{M}(\hat{f}(x_{+j}^{m})-\hat{f}(x_{-j}^{m}))$$

where $\hat{f}(x_{+j}^m)$ is the prediction for $x$, but with a random number of feature values replaced by feature values from a random data point $z$, except for the respective value of feature $j$.

The algorithm: Approximate Shapley estimation for single feature value

Output: Shapley value for the value of the $j$-th feature

Required: Number of iterations $M$, instance of interest $x$, feature index $j$, data matrix $X$, and machine learning model $f$

For all $m = 1,…,M$:

• Draw random instance $z$ from the data matrix $X$
• Pick a random subset of feature column indices $o$ (with $j$ not in $o$).
• Construct two new instances
  • With $j$ from $x$: $x_{+j}$, where all values in $x$ with index in $o$ are replaced by the respective values in $z$.
  • Without $j$ from $x$: $x_{−j}$, where all values in $x$ with index in $o$ are replaced by the respective values in $z$ and also the value for $j$ is replaced by the value in $z$.
• Compute marginal contribution: $\phi_j^m=\hat{f}(x_{+j})−\hat{f}(x_{−j})$

Compute Shapley value as the average: $\phi_j(x)=\frac{1}{M}\sum_{m=1}^M\phi_j^m$

Let’s code it

First we get a simple dataset. I picked the breast cancer wisconsin dataset from scikit-learn for this. Then we train a logistic regression model on the task.

import random
import warnings
import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline

warnings.filterwarnings("ignore")

# Load the dataset
features, target = load_breast_cancer(return_X_y=True, as_frame=True)

# Make a train/test split using 30% test size
X_train, X_test, y_train, y_test = train_test_split(
    features, target, test_size=0.30, random_state=2022
)

# prepare a model
f = make_pipeline(StandardScaler(), LogisticRegression())

# fit the model
f.fit(X_train, y_train)

features.columns

Index(['mean radius', 'mean texture', 'mean perimeter', 'mean area',
       'mean smoothness', 'mean compactness', 'mean concavity',
       'mean concave points', 'mean symmetry', 'mean fractal dimension',
       'radius error', 'texture error', 'perimeter error', 'area error',
       'smoothness error', 'compactness error', 'concavity error',
       'concave points error', 'symmetry error', 'fractal dimension error',
       'worst radius', 'worst texture', 'worst perimeter', 'worst area',
       'worst smoothness', 'worst compactness', 'worst concavity',
       'worst concave points', 'worst symmetry', 'worst fractal dimension'],
      dtype='object')

We investigate the feature ‘compactness error’ here, so $j$ will be $15$. The procedure is the same for every other feature as well. Also, we pick the number of iterations $M = 1000$. Next we also have to pick a sample instance to explain.

sample_idx = 0
x = X_test.iloc[sample_idx]

# calculate the shaply value for feature j
j = 15
M = 1000
n_features = len(x)
marginal_contributions = []
feature_idxs = list(range(n_features))
feature_idxs.remove(j)
for _ in range(M):
    z = X_train.sample(1).values[0]
    x_idx = random.sample(feature_idxs, min(max(int(0.2*n_features), random.choice(feature_idxs)), int(0.8*n_features)))
    z_idx = [idx for idx in feature_idxs if idx not in x_idx]
    
    # construct two new instances
    x_plus_j = np.array([x[i] if i in x_idx + [j] else z[i] for i in range(n_features)])
    x_minus_j = np.array([z[i] if i in z_idx + [j] else x[i] for i in range(n_features)])
    
    # calculate marginal contribution
    marginal_contribution = f.predict_proba(x_plus_j.reshape(1, -1))[0][1] - f.predict_proba(x_minus_j.reshape(1, -1))[0][1]
    marginal_contributions.append(marginal_contribution)
    
phi_j_x = sum(marginal_contributions) / len(marginal_contributions)  # our shaply value
print(f"Shaply value for feature j: {phi_j_x:.5}")

Shaply value for feature j: -0.026152

Compare to shap values

We use the python package shap to compare the shapley values we estimated to the estimate of a well-established software. Note, that the shap package actually uses a different method to estimate the shapley values.

import shap

# explain the model's predictions using SHAP
explainer = shap.KernelExplainer(f.predict_proba, X_train)
shap_values = explainer.shap_values(X_test.iloc[sample_idx,:])

Let’s see how we did!

print(f"Shaply value calulated from shap: {shap_values[1][j]:.5}")

Shaply value calulated from shap: -0.021179

So, it looks like our method did really well for that feature (and the specific instance)! Go ahead and try it for other features and samples. For sure this implementation is not the fastest and most elegant solution, but it helped my understanding of the method and might benefit yours as well. Also have a detailed look at a different, popular method for explaining machine learning models called LIME in this article.

Further reading:

• I highly recommend checking out the awesome, free book on interpretable machine learning, and especially the chapter on shapley values: https://christophm.github.io/interpretable-ml-book/shapley.html