SHAPでモデルを解釈してみた

はじめに
SHAPとは
ライブラリについて
インストール
データセット
モデル作成
- LightGBM
- Xgboost
SHAP Value
- Violin Plot
- Dependence Plot
Monotonic Constraints
まとめ

はじめに

XAI（Explainable AI）という言葉を聞いたことはありますでしょうか．
日本語では「説明可能なAI」と呼ばれていて，構築した学習モデルが入力に対してどのように出力に起因したのかというモデルの解釈性を示す分野として注目を集めております．機械学習のモデルはブラックボックスとなりがちなため，なぜそのような結果になったのかという解釈は一般的には難しい側面があります．こう言った課題はアカデミックな領域ではあまり問題とならない場合もあるのですが，ビジネスの現場では機械学習のアルゴリズム導入の妨げとなる要因になります．というのもシステムを構築してサービスを提供する企業において，想定していない挙動をとることに対して信頼性や安全性と言った観点で難色を示す場合が往往にしてあるからです．
そのため，今回は機械学習の解釈性という非常に重要なトピックを題材に，モダンな機械学習アルゴリズムと親和性が高く扱いやすいSHAPについて構築したモデルにおける活用方法を含めて説明します．

SHAPとは

SHAPとはブラックボックスしがちな予測モデルの各変数の寄与率を求めるための手法で，各特徴量が予測モデルの結果に対して正負のどちらの方向に対してどれくらい寄与したかを把握することによって予測モデルの解釈を行えるようになります． SHAPの理論は主に協力ゲーム理論におけるShaply Valueが由来となっており，協力ゲーム理論におけるプレイヤーを特徴量に置き換えて特徴量毎の貢献度(Shapley Value)を求めています．shapley valueの計算式は以下のように定義されます．

$\displaystyle{ \phi_i = \sum_{S ⊂ F\ \{i\}} \frac{|S|!(|F| - |S| - 1)!}{N!} [ f_{S∪\{i\}}\left(x_{S∪\{i \}}\right)-f_{S}\left(x_{S}\right)] }$

この式の $\frac{|S|!(|F| - |S| - 1)!}{N!}$ の部分が組み合わせ数の逆数となっており，後半部分の $[ f_{S∪{i}}(x_{S∪{i}})-f_S(x_S)$ ]はある組み合わせにおけるプレイヤーiの貢献度合いを示していることより，全ての組み合わせにおける貢献度合いの平均を算出していることになります．このようにして全ての組み合わせにおける貢献度を求めることでShapley Valueが計算されることがわかります．
ただし，全ての場合でshapely valueを計算することができれば良いのですが，現実的な問題として特徴量が増えていった際に，特徴量同士の組み合わせ数が膨大になり有限な時間で計算が終わらないといった事象が発生します．そのため，実際には近似値を計算するなどの工夫が行われています．

より詳細な中身のロジックについては本家の論文を参照してみてください．

ライブラリについて

ここまで述べてきたSHAPの理論についてですが，この理論を実装したSHAPと言うライブラリがあります．

このライブラリでは最適化された計算手法やモデル解釈のための可視化メソッドが用意されており，我々が生成したモデルに対して簡単にSHAP Valueを計算することをサポートしてくれます．モデルにおける各特徴量の貢献度を計算する仕組みが内部実装されている関係で，サポートされていないアルゴリズムやライブラリでのモデルは扱うことができませんが，有名なアルゴリズムは既にサポートしているので特殊な場合を除きすぐに実装して計算することができます．

SHAPライブラリで解釈を行うことが可能な対応しているアルゴリズム(ライブラリ)について紹介します．
まず，木系のアルゴリズムでは以下のモデルがSHAPで扱うことができます．

XGBoost
LightGBM
CatBoost
scikit-learn tree model
pyspark tree model

SHAPライブラリでは上記の木系アルゴリズムを高速に処理するためにC++ライブラリとして実装されています．その他にもニューラルネットワークのフレームワークではTensorflowとKerasをサポートしており，一部の機能ではPytorchも利用できるようになっています．

実際にSHAPライブラリの具体的な利用方法を知りたい場合は本家のGithubを参照することをお勧めします．というのも実装例として多くのJupyter NotebookによるサンプルをREADMEに用意していることから参照することで実装と出力イメージが両方掴めるため，これからSHAPで初めて解釈を実施してみようと言う方は一読してみてください．

github.com

インストール

SHAPライブラリのインストール方法を実行環境の情報と併せて以下に記載します．

ProductName:  Mac OS X
ProductVersion: 10.14.6
Python:  3.6.9

今回はLightGBMとXgboostのモデルに対してSHAPによる解釈を実施していくため，まず下準備としてLIghtGBMとXgboostのライブラリをインストールします．

$ pip install lightgbm==2.3.1 xgboost==1.0.0

続いて，本命のSHAPをインストールです．

$ pip install shap==0.35.0

インストール自体はpipで非常に簡単にインストールすることができるので導入しやすいと思います．

データセット

扱うデータセットはUCI machine learning repositoryで公開されているデータセット一覧の中から「Wine Quality Data Set」を利用します．このデータセットはワインの品質に関する情報を含んでおり，白ワインと赤ワイン別にサンプルデータが分割されているためフラグを自身で付与することで2値分類用のデータセットを作成することができます．

https://archive.ics.uci.edu/ml/datasets/Wine+Quality

本記事ではこのデータセットを用いて品質に関する情報から白ワインと赤ワインの2値分類の問題としてモデルを作成し，SHAPによるモデルの解釈を行います．

具体的なカラムを参考情報として以下にまとめておきます．

カラム名	データ型	説明
sulfur dioxide	float64	二酸化硫黄
chlorides	float64	塩化化合物
density	float64	濃度
sulphates	float64	亜硫酸塩
volatile acidity	float64	揮発酸(酢酸)
pH	float64	水素イオン指数
residual sugar	float64	残糖量
alcohol	float64	アルコール
fixed acidity	float64	酒石酸濃度
free sulfur dioxide	float64	遊離亜硫酸濃度
citric acid	float64	クエン酸濃度
quality	int64	ワインの味

モデル作成

続いてモデル作成を実際のコードと併せて説明していきます．

今回はLightGBMとXgboostの2つのそれぞれでモデルを作成します．
この2つのモデルを選択した深い意味は特にはないのですが，私がベースラインのモデルを作成する際に頻繁に利用するアルゴリズムということで2つのモデルを利用しました．

まずは各種必要なライブラリのインポートです．
基本的な分析関連のライブラリをインポートしておきます．可視化ライブラリはmatplotlib，データ操作に関してはpandas，numpy，scikit-learnを利用します．

import numpy as np
import pandas as pd
import xgboost as xgb
import lightgbm as lgb
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score, roc_curve, auc
from sklearn.metrics import classification_report

先ほど解説したUCIのワインに関するデータをpandasを利用してDataFrame形式で読み込みます．白ワインと赤ワインを識別するカラムはデフォルトのcsvファイルには含まれていないため，自分で作成する必要があります．そのため，今回は新たに class というカラムを用意して赤ワインと白ワインを識別するフラグを用意しています．

# データの読み込み
wine_red = pd.read_csv("../data/winequality-red.csv", sep=";")
wine_white = pd.read_csv("../data/winequality-white.csv", sep=";")

# classカラムを追加( 赤ワインが 0, 白ワインが 1 )
wine_red["class"] = 0
wine_white["class"] = 1

# 2つのデータセットを結合
wine = pd.concat([wine_red, wine_white], axis=0)

作成したデータを学習データセットとテストデータに分割します．
全体の2割をテストデータとして利用します．

# 固定パラメータの指定
SEED = 42
TEST_SIZE = 0.2
target = "class"
predictors = [i for i in wine.columns if i not in target]

# データの分割
X_train, X_test, y_train, y_test = train_test_split(
    wine[predictors],
    wine[target],
    test_size=TEST_SIZE,
    random_state=SEED
)

LightGBM

ここからはLightGBMでのモデル学習と精度確認をしていきます．
ハイパーパラメータは2値分類で学習するための値を指定して学習します．イテレーションの回数は5000回として，early stoppingとして100を設定してモデルを学習します．

# Lightgのハイパーパラメータ
lgb_params = {
    'task': 'train',
    'max_depth': 10,
    'metric': 'auc',
    'objective': 'binary',
    'boosting': 'gbdt',
    'learning_rate': 0.005,
    'seed': SEED,
    'num_threads': 8,
    'early_stopping_round': 100,
    'subsample': .8,
    'min_gain_to_split': .1, 
    'lambda_l1': 1,
    'lambda_l2': 0,
    'min_child_weight': .8,
    'feature_fraction': .8,
    'num_leaves': 31
}

# モデルの学習
lgb_train = lgb.Dataset(X_train, y_train)
lgb_eval = lgb.Dataset(X_test, y_test)
lgb_model = lgb.train(
    params=lgb_params,
    train_set=lgb_train,
   num_boost_round=5000,
    valid_sets=lgb_eval
)

# テストデータによる推論
y_pred = lgb_model.predict(X_test, num_iteration=lgb_model.best_iteration)

# 評価指標の計算
fpr, tpr, thresholds = roc_curve(y_test, y_pred)
roc_auc = auc(fpr, tpr)

# 重要変数とROC-Curveの可視化
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(13,6))
feat_imp = pd.Series(lgb_model.feature_importance(), index=X_test.columns).sort_values(ascending=False)
feat_imp = feat_imp.head(30)
feat_imp = feat_imp.sort_values(ascending=True)
feat_imp.plot(kind='barh', title='Feature Importance', ax=axes[0], color=sns.color_palette()[0])
axes[0].set_xlabel('Feature Importance Score')
axes[1].plot(fpr, tpr, label='ROC curve (area = %0.3f)' % roc_auc)
axes[1].plot([0,1], [0,1], 'k--')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve')
legend = axes[1].legend(frameon=True, loc='lower right', fontsize = 'medium') # 凡例
frame = legend.get_frame()
frame.set_facecolor('white')
frame.set_edgecolor('gray')

学習したモデルの重要変数とROCを可視化したグラフは以下に示す通りになりました．ROC-AUCは99.7%と高精度なモデルを学習できているようです．続いて，重要変数を見てみると最も予測結果に寄与している変数はtotal sulfur dioxideということでワインの酸化防止剤として活用されて二酸化硫黄の含有率がワインの識別に最も関係があることが結果として見て取れます．逆に quality はワインの識別には関係ないことがわかりますね．

f:id:takaherox:20200701232345p:plain — LightGBM Feature Importance

Xgboost

続いてXgboostによるモデル作成と精度確認の部分になります．
ハイパーパラメータなどはLightGBMとほぼ同様になるように設定しています．

# xgboostのハイパーパラメータ
xgb_params = {
    "booster": "gbtree",
    "max_depth": 10,
    "objective": "binary:logistic",
    'eval_metric': 'auc',
    "eta": 0.005,
    "seed": SEED,
    "nthread": 8,
    "subsample": 0.8,
    "gamma": 0,
    "alpha": .1,
    "lambda": 0,
    "min_child_weight": .8,
    "colsample_bytree": .8,
    "max_leaves": 31
}

# モデルの学習
xgb_train = xgb.DMatrix(X_train, y_train)
xgb_test = xgb.DMatrix(X_test, y_test)
evals_result = {}
xgb_model = xgb.train(
    xgb_params,
    xgb_train,
    num_boost_round=5000,  
    early_stopping_rounds=500,
    evals = [(xgb_test, 'eval'), (xgb_train, 'train')],
    evals_result=evals_result
)

# テストデータによる推論
y_pred_xgb = xgb_model.predict(xgb_test)

# 評価指標の計算
fpr, tpr, thresholds = roc_curve(y_test, y_pred_xgb)
roc_auc = auc(fpr, tpr)

# 重要変数とROC-Curveの可視化
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(13,6))
feat_imp = pd.Series(xgb_model.get_fscore()).sort_values(ascending=False)
feat_imp = feat_imp.head(30)
feat_imp = feat_imp.sort_values(ascending=True)
feat_imp.plot(kind='barh', title='Feature Importance', ax=axes[0], color=sns.color_palette()[0])
axes[0].set_xlabel('Feature Importance Score')
axes[1].plot(fpr, tpr, label='ROC curve (area = %0.3f)' % roc_auc)
axes[1].plot([0,1], [0,1], 'k--')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve')
legend = axes[1].legend(frameon=True, loc='lower right', fontsize = 'medium') # 凡例
frame = legend.get_frame()
frame.set_facecolor('white')
frame.set_edgecolor('gray')

当然ですが，こちらもLightGBMの結果と同様にほとんど同じ結果になっています．精度は誤差の範囲ですが99.9%とLightGBMより高い精度のモデルになりました．

f:id:takaherox:20200701232426p:plain — Xgboost Feature Importance

SHAP Value

ここからが本題のSHAPになります．
それぞれのモデルでSHAP Valueの計算と可視化をしてみます． TreeExplainerクラスにてモデルを入力してインスタンスを生成し，shap_valuesメソッドを実行するだけで簡単にSHAP Valueを計算することができます． SHAPライブラリ内ではLightGBMのpredictメソッドを実行しており，その際にパラメータの pred_contrib を指定してSHAP Valueを求めています．また，Xgboostも同様の方法でSHAP Valueを計算しています．

# LightGBMのSHAP Valueの計算
lgb_explainer = shap.TreeExplainer(lgb_model, X_train)
lgb_shap_values = lgb_explainer.shap_values(X_train)

# XgboostのSHAP Valueの計算
xgb_explainer = shap.TreeExplainer(xgb_model)
xgb_shap_values = xgb_explainer.shap_values(X_train)

SHAP Valueを計算した後は可視化をします．
可視化にはいくつかの種類がありますので自身で適切な可視化を選択する必要があります．試しにViolinPlotとDependencePlotで結果を見てみます．SHAPライブラリでは複数の可視化パターンを用意しているためGithubを参考にして実施したい解釈に合わせて適切な可視化パターンを選択してみてください．

Violin Plot

# LightGBMでのViolin Plotの可視化
shap.summary_plot(lgb_shap_values, X_train)

# XgboostでのViolin Plotの可視化
shap.summary_plot(xgb_shap_values, X_train)

f:id:takaherox:20200704234439p:plain — Violin Plot(左がLightGBM, 右がXgboost)

f:id:takaherox:20200705105807p:plain — Violin Plot(左がLightGBM, 右がXgboost)

Violinプロットを見てみるとパラメータの値が結果に対してどのように寄与しているのかを把握することができます．最上位にきている total sulfer dioxide を見てみると，この値が大きいほど予測結果を白ワインである確度が高くなるように作用していることがわかります．一方で chlorides を見てみると逆に作用しており，値が大きいほど赤ワインである方向に作用しているようです．このように属性情報の値の大きさがどちらに作用しているのかを確認することでモデルの入力と出力の関係を理解することに繋げられます．

Dependence Plot

# LightGBMでのPartial Dependence Plotの可視化
shap.dependence_plot("rank(0)", lgb_shap_values, X_train)

# XgboostでのPartial Dependence Plotの可視化
shap.dependence_plot("rank(0)", xgb_shap_values, X_train)

f:id:takaherox:20200704234449p:plain — Partial Dependency Plot(右がLightGBM, 左がXgboost)

f:id:takaherox:20200705105852p:plain — Partial Dependency Plot(右がLightGBM, 左がXgboost)

Dependence Plotを見てみると，単一の特徴量に対して全てのパラメータがどのように予測結果に影響しているのかを理解することができます．また，別のパラメータとの交互作用も合わせて確認することができます．交互作用を確認したいパラメータを指定する場合は interaction_index を引数で指定してあげれば簡単に可視化できます．
今回の結果では total sulfer dioxide の値が大きくなるとSHAP Valueが大きくなっていることが見て取れます．LightGBMもXgboostでも同様の結果になっていますね．このように単一の特徴量における予測結果への影響が線型的に単調増加している場合は解釈が容易なのですが，必ずしもそういう結果になるとは限りません．次の項目では解釈をしやすくするためのMonotonic Constraints(単調制約性)について簡単に記載しておきます．

Monotonic Constraints

Monotonic Constraintsとは単調制約性のことを指します．つまり，一般的に特徴量が与えるモデルへの影響は複雑な作用をする場合が往々にしてありますが，特徴量がモデルへの与える影響を単調増加(減少)になるように制約を与えるということになります．

今回の例では単調制約性を与えてないない場合は citric acid は以下のように0.4まではSHAP Value値が大きくなりますが，0.4以降ではSHAP Valueが小さくなっていきます．このような場合にMonotonic Constraintsを与えた場合にどのような効果が得られるのかを確認して見ましょう．

f:id:takaherox:20200705202243p:plain — Xgboost Dependence Plot - cirtic acid

Monotonic Constraintsを設定する方法は簡単で学習時のパラメータに制約項目を指定するだけです．事前に特徴量が目的変数に対してどのような相関があるのかを調べておき，負の相関，無相関，正の相関のそれぞれに応じて -1, 0, 1 を指定します．今回は正の相関を指定するため1を指定します．

# Xgboostのハイパーパラメータ
xgb_params = {
    "booster": "gbtree",
    "max_depth": 10,
    "objective": "binary:logistic",
    'eval_metric': 'auc',
    "eta": 0.005,
    "seed": SEED,
    "nthread": 8,
    "subsample": 0.8,
    "gamma": 0,
    "alpha": .1,
    "lambda": 0,
    "min_child_weight": .8,
    "colsample_bytree": .8,
    "max_leaves": 31,
    "monotone_constraints": "(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 )"  #  単調制約性
}

# モデルの学習
xgb_model_mono = xgb.train(
    xgb_params,
    xgb_train,
    num_boost_round=5000,  
    early_stopping_rounds=500,
    evals = [(xgb_test, 'eval'), (xgb_train, 'train')],
    evals_result=evals_result
)

# SHAP Valueの可視化
xgb_mono_explainer = shap.TreeExplainer(xgb_model_mono)
xgb_mono_shap_values = xgb_mono_explainer.shap_values(X_train)
shap.dependence_plot("rank(9)", xgb_mono_shap_values, X_train)

先ほどの図と比べても変化は一目瞭然ですね．
パラメータで指定した正の相関のように値の増加に伴ってSHAP Valueが単調増加していることが分かります．このように解釈をしやすいようにモデルに制約を与えることもできます．この制約を与えてモデルを作成した場合には予測精度が低下する場合もあるため，モデル精度を下げたくない場合は注意してください．

f:id:takaherox:20200705202252p:plain — Xgboost Dependence Plot - cirtic acid monotonic

まとめ

今回はSHAPの簡単な理論と木系のアルゴリズムにてSHAPの扱い方をモデルの作成からSHAP Valueの可視化までを通じて説明しました．XAIはデータ分析においても非常に重要なトピックでありますし，ビジネスサイドでも役に立つということで是非活用してみてください．
SHAPライブラリはDeepLearningのフレームワークでも扱うことができるため，今後はそちらも試して見たいと思います．

tkherox blog

データサイエンスおよびソフトウェア開発、たまに育児についての話を書いています