決定木ベースの機械学習手法におけるSHAPとfeature importanceの違い

# はじめに
閲覧ありがとうございます。今回は決定木ベースの機械学習手法で用いることのできる、結果を解釈する指標となるSHAPと呼ばれる値の算出方法とその結果の解釈方法を共有できればと考えています。
最後までお付き合い頂けますと嬉しいです！

# 決定木ベースの機械学習手法
まずはそもそも「決定木ベースの機械学習」とは何か？について説明します。
決定木は以下の図で示すような手法です。
```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 25
}
</style>
center header 図1 決定木モデルの例
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/01dce6b1-2c77-db5c-b124-f9bb87c2ea60.png>
@enduml
```

訓練データをある特徴量の値の大小で分割し、ラベリングを行う方法です。この分割は分割後のデータ集合に含まれるラベルのばらつきが最も小さくなるように行われます。
例えば、「家の価格」を予測するとき、土地の広さが50m^2以上or未満であったり、その地域の人口が1000人以上or未満であったりといった観点でデータを分割していきます。

この決定木から派生した「決定木ベース」の機械学習が存在します。大きく「バギング」、「ブースティング」のいずれかに分類されます。
詳細な説明は割愛しますが、イメージは以下です。

### バギング
独立した複数の決定木を作り、それらの出力を平均する手法。
```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図２　バギングの例（ランダムフォレスト）
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/34c2282e-9aa6-3468-eb75-c6ac5966fa36.png>
@enduml
```
### ブースティング
単独の決定木の出力を次の出力として学習するということを繰り返し行う手法。
```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図3　ブースティングの例（XGBoost, LightGBM）
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/56a41fca-c0da-e170-fc99-120090699ab6.png>
@enduml
```
今回はブースティングの手法であるLightGBMを学習器として用い、二値分類の問題に対しSHAP及びfeature importanceを比較します。


# SHAPとは
協力ゲーム理論に基づく機械学習モデルの解釈手法。
SHAPは「予測用データと学習済みモデル」を用いて、特徴量ごとに算出される値です。
モデルが予測値の算出に大きく寄与するような（モデルの出力を大きく変化させるような）特徴量を表現できます。

SHAPと協力ゲーム理論の関係性については[DataRobotのサイト](URL "https://www.datarobot.com/jp/blog/explain-machine-learning-models-using-shap/")が非常に参考になります。

算出方法については[「SHAPの論文を読んでみた」](URL "https://qiita.com/ground0state/items/ff8dd27edd210ae4468c")の内容及び、記載の文献が参考になります。

# feature importanceとは
決定木系モデルにおいてfeature importance(特徴量重要度)は一般的に用いられる解釈指標です。
この指標は「学習済みモデル単体」から計算され、特徴量ごとに算出される値です。
学習時のサンプル分割に計算された**ジニ不純度の減少**を用いて計算されます。

ジニ不純度はサンプル集合(ノード)内に存在する分類ラベルの混ざり具合を表します。
あるサノードが、「すべて同じラベルの場合」は0となり、「複数のラベルが混在する」ほど1に近づく値で以下の式で計算されます。

```math
G(r) = \sum_{i=1}^k p(i) × (1 - p(i))
```
| 変数 | 説明 |
|:-----------|------------:|
| G(r)  | ジニ不純度 | 
| k     | 分類ラベルの種類数  |
|p(i)   | ターゲットとなるラベルの出現頻度 |


https://yolo-kiyoshi.com/2019/09/16/post-1226/#outline__1

# 実験
訓練データで学習した学習済みモデルに対し、「モデルのfeature importance」及び「モデルとテストデータから算出されるSHAP」を比較します。

今回は再現性の確保のためにsklearnのサンプルデータセットである、lord_wine()データセットを用いて、多値(3値)分類モデルにおいて実験を行いました。

環境は以下です。
・　macOS　13.2.1
・　Python 3.11.7

### データに関する説明
元のデータは13個の説明変数と1つの目的変数からなります。
目的変数はclass_0,class_1,class_2の3クラスからなります。

目的変数別のサンプル数は(class_0,class_1,class_2)=(59,71,48)となっており、サンプル数に若干の偏りがあるデータとなっています。

今回のコードでは、訓練データは(class_0,class_1,class_2)=(40,50,34)、テストデータでは(class_0,class_1,class_2)=(19,21,14)となっています。

# 実装
### データ読み込み ~ モデリング
```Python
# 実行環境 Python3.11.7

# import modules
import matplotlib.pyplot as plt
import pandas as pd
import shap
from lightgbm import LGBMClassifier as lgbmc
from sklearn.datasets import load_wine as wine
from sklearn.metrics import classification_report as cr
from sklearn.model_selection import train_test_split as tts

# データセットの読み込み
data = wine()
X = data.data
y = data.target

# データサイズを確認
print(X.shape, y.shape)

# データ分割
X_train, X_test, y_train, y_test = tts(X, y, test_size=0.3, random_state=42)

# モデルの学習
model = lgbmc()
model.fit(X_train, y_train)

# モデルの評価
y_pred = model.predict(X_test)
print(cr(y_test, y_pred))
```
モデルの予測精度は以下です。どの分類ラベルにおいてもテストデータに対し、高い予測精度を持っていることがわかります。
```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図4 学習済みモデルのテストデータに対する予測精度
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/28c46afa-71ab-3d58-130e-fc54fda964c6.png>
@endum
```

```Python
# SHAP値の計算
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# SHAP値のプロット（3クラスのbar_plot）
shap.summary_plot(shap_values, X_test, feature_names=data.feature_names)

# SHAP値のプロット(特定クラスのbesswarm_plot)
# こちらを適用する場合コメントアウトを外す
# また、shap_valuesの対象クラスを希望のものに変更する
# shap.summary_plot(shap_values[0], X_test, feature_names=data.feature_names)

```
SHAPの結果は以下の図で出力されます。

```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図5 学習済みモデルとテストデータで算出したSHAP値のbar_plot
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/f334096a-6afb-e442-269d-d8e9d0efa6c2.png>
@endum
```

分類ラベル・特徴量ごとにSHAPの値が計算されます。
各分類ラベルにおけるSHAPの総和で降順にソートされていますが、その内訳は一定ではありません。

では次にfeature_importanceを算出します。

```Python

# feature_importanceのプロット
importance = pd.DataFrame(model.feature_importances_, index=data.feature_names)
importance.sort_values(ascending=True,by = 0, inplace=True)
plt.barh(importance.index, importance[0])

```

```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図6 学習済みモデルから算出したfeature_importance
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/b2597224-f0fe-e047-0bf2-9402575355eb.png>
@endum
```
feature importanceでは「分類ラベルごとの予測に重要な特徴量」ではなく、「学習時のサンプルの分割に重要な特徴量」が表現されることになります。
そのためSHAPとは上位の特徴量の順位が異なっていることがわかります。
一方で上位の特徴量の種類自体は大きな変化がありません。

# 解釈
では、これらの性質の違いを比較します。
今回は各指標で最も重要であると判断された「flavanoids」と「proline」の２つの特徴量について、より詳細にデータを見ていきましょう。

まずは縦軸に「proline」、横軸に「flavanoids」の値をプロットし、点の色をラベルで色分けした図を出力します。

``` Python
# ２つのカラムの値で散布図をプロット（クラスで色分け）
data = wine()
c1 = "flavanoids"
c2 = "proline"

data = pd.concat([pd.DataFrame(X, columns=data.feature_names), pd.DataFrame(y, columns=["target"])], axis=1)
data_0 = data[data["target"]==0]
data_1 = data[data["target"]==1]
data_2 = data[data["target"]==2]

plt.figure()
plt.scatter(data_0[c1], data_0[c2], label="class_0", c="red")
plt.scatter(data_1[c1], data_1[c2], label="class_1", c="blue")
plt.scatter(data_2[c1], data_2[c2], label="class_2", c="green")
plt.xlabel(c1)
plt.ylabel(c2)
plt.legend()
plt.show()
```

```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図7 proline-flavanoidsの散布図
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/702101fb-4f08-807d-1c2a-f8756b4c801d.png>
@endum
```
### 散布図から分かること
まずは「proline」について注目してみましょう。
図を見ると、class_0とclass_1,class_2の分割には寄与していますが、class_1とclass_2の分割にはあまり寄与していないことが見て取れます。

対して「flavanoids」は、class_2とclass_0,class_1の分割には寄与していますが、class_0とclass_1の分割にはあまり寄与していないことが見て取れます。
SHAPの結果の図（図5）を見ても、class_2の割合が大きくなっていますね。

### 分類ラベル別のSHAPの結果
上のコードでコメントアウトで示した部分を利用して、クラス別にSHAPの値を見てみます。
以下のコードではbesswarm_plotというものが出力され、サンプルごとのSHAPを見ることができます。
```Python
# SHAP値の計算
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# SHAP値のプロット(特定クラスのbesswarm_plot)
# こちらを適用する場合コメントアウトを外す
# また、shap_valuesの対象クラスを希望のものに変更して実行する
shap.summary_plot(shap_values[0], X_test, feature_names=data.feature_names)
# shap.summary_plot(shap_values[1], X_test, feature_names=data.feature_names)
# shap.summary_plot(shap_values[2], X_test, feature_names=data.feature_names)

```
class_0の結果から見ていきます。
```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図8 class_0に関するSHAPのbasswarm_plot
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/7929e905-0749-c4f6-c64a-af321b406ac9.png>
@endum
```
縦軸は特徴量の名前、横軸は各サンプルのSHAPを表します。各点は1サンプルを表します。
各点は青から赤のグラデーションで表され、赤色が強いほどそのサンプルでの特徴量は大きな値をとります。
右方に赤い点、左方に青い点が集中している場合、「その特徴量が大きいほど、そのクラスの予測確率が増加する」、右方に青い点、左方に赤い点が集中している場合、「その特徴量が大きいほど、そのクラスの予測確率が減少する」と解釈ができます。理想的な場合には、赤い点と青い点の距離が離れることが期待されます。

「proline」についてみると、SHAPは値の大小によってよく分離しており、大きな値を取るほどclass_0の予測確率が増加することが見て取れます。また「alcohol」や「flavanoids」も同様の傾向が見られます。
class_0の予測には複数の要素が重要なようです。

次にclass_1の結果です。

```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図9 class_1に関するSHAPのbasswarm_plot
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/5ffb519d-3365-348d-ac62-44aaf271af1f.png>
@endum
```

class_1では「color_intensity」が特に重要な特徴量であり、特に小さい値を取るときに予測確率を増加させる傾向があるようです。
class_0でも重要な特徴量であった「alcohol」,「proline」も予測にある程度寄与するようです。
最後にclass_2についての結果です。


```plantuml
@startuml
<style>
root {
  Padding 0
  Margin 10
}
header {
  FontColor black
  FontSize 17
}
</style>
center header 図10 class_2に関するSHAPのbasswarm_plot
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603842/3ba9001e-b4ed-7186-533d-0063b091d66b.png>
@endum
```
「flavanoids」について、class_0とは反対に値が小さくなるほど、予測確率が増加する傾向にあるようです。
また他のクラスではあまり重要でなかった「hue」が予測に寄与するようです。


# 最後に
今回はSHAPとfeature_importanceの違いについて、標準的なデータセットを用いて実験を行いました。
feature_importanceと比較して、SHAPは「分類クラスごとの重要度が算出できる」、「テストデータに対する予測の解釈が行える」などの長所があり、予測についてのより詳細な情報が得られると考えられます。一方でfeature_importanceはラベル全体を考慮した上で重要な特徴量を選択できると言えます。
モデルの解析の目的によって、どちらを見るのかは変わってきそうです。

SHAP・feature_importanceを利用したより良い解釈や、より優れた解釈指標がありましたら教えていただけますと幸いです。

# 参考文献

この記事の作成に際し、以下の記事を参考にさせていただきました。
[入門]初心者の初心者による初心者のための決定木分析
https://qiita.com/3000manJPY/items/ef7495960f472ec14377
SHAPで因果関係を説明
https://qiita.com/s1ok69oo/items/0bf92b84e565789a2191
Shapを用いた機械学習モデルの解釈説明
https://qiita.com/shin_mura/items/cde01198552eda9146b7

ここまでお付き合いいただきありがとうございました。