LightGBMの出力結果を解析したい！(SHAPのススメ)

1.はじめに

SIGNATEで「日本取引所グループ ファンダメンタルズ分析チャレンジ」というコンペが開催されています。私も参加していますが、その中で出てくる知識に関して基礎部分をまとめよう！という動機で書いた記事の第2弾です。
コンペのチュートリアルでは、「XGBoost」での予測結果に関しての根拠説明として本記事で紹介する「SHAP」を使用していますが、本記事では現在一番メジャーなLightGBMを使用してこのSHAPを紹介していきます。

• 動作環境
  • OS : Windows10 pro
  • Python : 3.8.3 // Miniconda 4.9.1
  • LightGBM: 3.1.1
  • SHAP : 0.38.1
  • jupyter notebook

2.LightGBMでとりあえず学習して、特徴量重要度を図示する

今回は、回帰分析を例として行うことにする。
サンプルデータはscikit-learnに添付されているボストンの住宅価格データを使用する。
LightGBMは本記事のメインではないのでサラッと流す。もしもっと知りたい等あればコメントしてください。別記事で掘り下げます。

import numpy as np
from sklearn.datasets import load_boston
import matplotlib.pyplot as plt
import pandas as pd
import lightgbm as lgb
from sklearn.model_selection import train_test_split

#ボストンのデータセットを読み込む
boston = load_boston()

#読み込んだデータセットをデータフレーム形式に変更
dataset = pd.DataFrame(data = boston['data'], columns = boston['feature_names'])
dataset['price'] = boston['target'] #住宅価格がtargetという名前なのでわかりやすくpriceに変更
dataset.head(2) #データフレームの最初2行を表示する

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT	price
0	0.00632	18.0	2.31	0.0	0.538	6.575	65.2	4.0900	1.0	296.0	15.3	396.9	4.98	24.0
1	0.02731	0.0	7.07	0.0	0.469	6.421	78.9	4.9671	2.0	242.0	17.8	396.9	9.14	21.6

一応おさらいすると、予測対象は上表テーブルの「price(住宅価格)」である。

#説明変数と目的変数に分割する
X = dataset[dataset.columns[dataset.columns != 'price']] #price以外の変数を抜き出してXとする
y = dataset['price']

#train_test_splitを使用して「学習8割,テスト2割」にデータを分割する
X_train, X_test, y_train, y_test = train_test_split(
                                    X, y, test_size=0.2, random_state=0
                                )

#LightGBMに与えるパラメータ。回帰なので「regression」それ以外適当です
params = {
    'boosting_type': 'gbdt',
    'objective': 'regression',
    'metric': 'rmse',
    'random_state': 0,
    'verbose' : -1,
}

#入力データセット
lgb_train = lgb.Dataset(X_train, y_train)
lgb_val = lgb.Dataset(X_test, y_test, reference=lgb_train)

#学習
gbm = lgb.train(params,
                lgb_train,
                valid_sets=[lgb_train, lgb_val],
                num_boost_round = 10000,
                early_stopping_rounds=100,
                verbose_eval=100)

これで「gbm」に学習したモデルが入っている状態ある。

最後に何が学習に貢献したのか？を図示する(importance)

lgb.plot_importance(gbm, figsize=(8,4), max_num_features=5, importance_type='gain')

3.SHAPで判断根拠を可視化(結果解釈)する

今回はSHAPの理論には触れない。詳細はここに詳しく書いてあるので参照してほしい。
機械学習モデルの予測値を解釈する「SHAP」と協力ゲーム理論の考え方

簡単に言うと、とある特徴量があった場合となかった場合にどれくらいスコアが変わるのか？というゲーム理論の考えを用いている

3-1. Feature importanceとの違い

最初に、「Feature importance」と何が違うの？？と思われる方がいると思うので触れておく。
Feature importanceとは、学習において以下の2パターンの改善に役立った指標のランキングであり、
・木の分岐回数が多かった特徴量(importance_type = 'split')
・gini係数の改善に役立った特徴量(importance_type = 'gain') ※上図のinportanceはこれ

要するに簡単に言えば「学習時」のlossの改善に何が貢献したか？を表しているといえる。
ただし、例えば未知のテストデータをこの学習モデルに与えたときに出力した結果が「想定と違うんだけど・・」と言われたときにはこのimportanceでは説明がつかないのである。
※出力結果がこうなった理由は＊＊の特徴量の数字が大きく、その特徴量のFeature importanceが高いから・・みたいに直接的にはは説明はできない。
そんな時こそSHAPの出番なのである(コンペよりも実務とかでこそこういう質問は多いはず）

3-2. SHAPの導入～SHAP値出力まで

SHAPはpip install shapで簡単に導入できる。

※(注)もし以下のコードでうまくいかない場合は5.追記(2021/12/30) も参考ください

#pip install shap

import shap

#notebook内でJavascriptを動かすためのおまじない
shap.initjs()

"""
shap.TreeExplainer:決定木用(XGBoost、lightBGM等含む)
shap.LinearExplainer :線形モデル用
shap.DeepExplainer :Deeplearning用
"""
#TreeExplainerは、決定木系のモデルのSHAP値を取得するもの。
explainer = shap.TreeExplainer(model=gbm)

print(explainer.expected_value)

実行結果

22.66837213

この「explainer.expected_value」が学習モデルの基礎スコアなのである。
この基礎スコアに各特徴量の寄与が重なった結果、その予測値が計算されたという考え方なのである

"""
modelと解釈したいデータを渡す　
※なお、indexのresetは別に不要ですがtestデータのindexが何番かをいちいち確認するのが面倒なのでresetしてます。
"""
X_test_shap = X_test.copy().reset_index(drop=True)
shap_values = explainer.shap_values(X=X_test_shap)

print(X_test_shap.shape)
print(shap_values.shape)
print(shap_values[0])#テストデータの0番目の要素を出力

実行結果

(76, 13) 　※テストデータは76個、説明変数13個
(76, 13)　 ※SHAP値も同じ形をしている
[ 1.32674474e-01  5.77572033e-02  2.66389567e-01  1.35371898e-03
  2.65023634e-01 -1.82737578e+00  7.54980695e-01 -3.17615038e-01
  4.55483300e-02 -9.91601699e-01  9.43307157e-01  3.22479589e-01
  1.98247682e+00]　※テストデータの各特徴量(カラム)別に算出されたSHAP値

例としてテストデータの0番目要素をSHAPに入力すると0番目テストデータにおける各特徴量別のSHAP値(貢献スコア)が得られることがわかると思う。※shap_valuesの出力順番は元のカラムの並び順(X_test_shap.columnsで調べればわかる)

3-3. SHAPの可視化

さて、求めたSHAP値をどう使ってどう図示するか？だが色々な方法がある。

(A) summary_plot

summary_plotでは結果出力にどの特徴量が大きく影響していたか？に関して図示してくれる
左図の赤い部分が「説明変数別の数字が大きかった部分」を指していて、逆に青い部分は「説明変数別の数字が小さかった部分」を示している。
LSTATの場合、赤いplotほどSHAP値が負になっていて青いplotは逆にSHAP値が高い。
つまり、LSTATが大きいほどスコアには負の影響、小さいほどスコアには正の影響(逆相関)だとわかる
また、'bar'を指定すれば、平均の貢献度順に右図のように図示をしてくれる。

そして先程のfeature_importanceと見比べると1位と2位の変数が入れ替わっていることがわかる。

shap.summary_plot(shap_values, X_test_shap) #左側の図
shap.summary_plot(shap_values, X_test_shap, plot_type='bar') #右側の図

(B) force_plot、waterfall_plot

force_plot(waterfall_plot)では、それぞれ個々のテストデータに対する具体的な貢献度を可視化できる。今回２つ例を出しているが、見やすい方を選べばいい(私のケースでは下図のwaterfall_plotの方が見やすいと言ってくれるケースが多かったので今回は載せていますが、一般的には上図の方が多そうです)。
うっすら「base_value=22.67」という記述が図で確認できるが、これは先程紹介した「学習モデルの基礎スコア」で、全データ共通の指標である。
その基本スコアに対して、各特徴量が正/負に貢献した結果最終的に「24.3」という予測をしたよ！という過程が一目でわかるようになっている。
※model予測結果 = model基礎スコア(explainer.expected_value) + 特徴量毎のSHAP値(shap_values)

※force_plotは１つずつではなく、まとめて複数サンプル見る方法もあるが今回は紹介しない
参考：SHAPのgithub

n = 0#テストデータのn番目の要素を指定

shap.force_plot(explainer.expected_value, shap_values[n,:], X_test_shap.iloc[n,:])#上の図

#waterfall_plotは私の環境ではエラーになるので、代わりにwaterfall_legacyを使用している。
shap.plots._waterfall.waterfall_legacy(explainer.expected_value, 
                                       shap_values[n,:], X_test_shap.iloc[n,:]) #下の図

例えば、このテストデータが「想定より予測値が低いけどなんで？」という質問が来た場合に「SHAPの指標を参考にすると、「TAX」と「RM」がスコアにブレーキをかけてることがわかります。「RM」のスコアがプラスになる条件はsummary_plotも参考にすればRM自体の特徴量が大きくないといけなそうなので・・・」みたいに説明できるかもしれません。

また、参考までに一応実モデルの予測結果も確認しておく

n = 0
print(gbm.predict(X_test_shap)[n])

実行結果

24.303770801258906  ※当然だがSHAPから算出した結果と一致している！

(C) decision_plot

decision_plotは基礎スコアから、各テストデータに対して決定木のように予測の過程を可視化することができる
下図では全データを図示しているが、個々に表示も可能。
ここで着目すべきは、どの説明変数あたりからスコアに差が出始めているのか？である。

shap.decision_plot(explainer.expected_value, shap_values,X_test_shap)

(D) dependence_plot

dependence_plotでは、変数間の関係性や、変数と予測値との関係性をより詳細にとらえられる。
y=axのグラフで、縦軸yがSHAP値、横軸xが特徴量というグラフで表される

LSTATの値が大きくなるほどShapley Valueが小さくなることが見て取れる。
色付けは「交互作用」が一番はっきりと表れる変数が自動で選ばれる。※今回は(CRIM)
CRIMが大きい数字なほど赤くなっていて、逆に小さいほど青い点が打たれる。
グラフから言えることは「同じLSTATの数字であった場合、CRIMが小さいほどLSTATのSHAP値が上昇する」ということである。

shap.dependence_plot("LSTAT", shap_values, X_test)

4.おわりに

どうだったであろうか？コンペや実務で、「なんでこの数字を予測したのか？モデルは馬鹿なのか？」と思ったときにこのツールを判断材料に使うと、判断材料の目安になりそうなのがわかっていただけたと思う。
実際に私も実務で顧客から「私はこのデータはこうなると思ってるんだけど、なんでこの予測結果なの？」という問いに対して、このSHAPをうまく活用して回答した経験があるので皆さんも気軽に使ってみてほしい。なお、回帰だけでなく分類問題にもこれは使用できるので調べてみてほしい。
※勿論SHAPだけで全ての現象を語れるわけではない。ただ参考値だとしてもブラックボックスよりマシなのである。

5.追記(2021/12/30)

この記事のコードだとうまくいかない例があるようです。(カテゴリがある時？)
※私のケースだとexplainer.expected_value=Noneになる
その場合は以下でお試しください。

import shap
shap.initjs()

explainer = shap.Explainer(gbm) #TreeExplainerでも多分変わらないが・・
shap_values = explainer(X_test_shap) #これで全確認データの数字を一気に入手
n = 0 #0番目のテストデータのSHAPを確認したい
print('shap値：', shap_values[n].values) #n番目テストデータのSHAP値
print('基礎スコア：', shap_values[n].base_values) #n番目テストデータの基礎スコア
print('説明変数値：', shap_values[n].data) #n番目テストデータの説明変数数値

実行結果

shap値： [ 1.32674474e-01  5.77572033e-02  2.66389567e-01  1.35371898e-03
  2.65023634e-01 -1.82737578e+00  7.54980695e-01 -3.17615038e-01
  4.55483300e-02 -9.91601699e-01  9.43307157e-01  3.22479589e-01
  1.98247682e+00]
基礎スコア： 22.66837212769739
説明変数値： [6.7240e-02 0.0000e+00 3.2400e+00 0.0000e+00 4.6000e-01 6.3330e+00
 1.7200e+01 5.2146e+00 4.0000e+00 4.3000e+02 1.6900e+01 3.7521e+02
 7.3400e+00]

上からSHAP値、基礎スコア、元データの順に数字が格納されています。
後は色んなshap命令に以下のように適応させていけばいい ※shap_values[n]を与えるだけでOK

#n番目のテストデータのwaterfall_plot
shap.plots.waterfall(shap_values[n]) #記事では_waterfall.waterfall_legacyだったが・・