imlパッケージによる機械学習モデルの解釈を試みた in R

##1.はじめに
説明可能なAIとして、機械学習モデルの解釈が最近いろいろと出てきているようです。
[Shapを用いた機械学習モデルの解釈説明](https://qiita.com/shin_mura/items/cde01198552eda9146b7)

SHAP（SHapley Additive　exPlnation）値は、特に決定木ベースのモデルから結果を解釈するための最新の方法で、ゲーム理論のShaply値がベースになっています。
[なぜ決定木だとSHAP値を正確に計算できるの?](https://qiita.com/hyt-sasaki/items/49355cb6e9775fe998a6)


Rでは、imlパッケージは、その他にも機械学習の解釈に用いられるPDP,ALEが使えますので早速使ってみたいと思います。
##2.データの準備
下記の続きから始めます。
[機械学習パッケージmlrを使った順序関係があるカテゴリ分類（その４）](https://qiita.com/kenkenvw/items/f330a4413540ae26e6f3)

BostonHousingの目的変数であるmedvをクラス分けし、３カテゴリ分類(medvC)を行いました。

```{r}
print(task_BH_train)

$>Supervised task: BostonHousing
$>Type: classif
$>Target: medvC
$>Observations: 354
$>Features:
$>   numerics     factors     ordered functionals 
$>         14           0           0           0 
$>Missings: FALSE
$>Has weights: FALSE
$>Has blocking: FALSE
$>Has coordinates: FALSE
$>Classes: 3
$>  0   1   2 
$>144 181  29 
$>Positive class: NA
```
この訓練データを使って、ハイパーパラメータを求めた結果が、次のとおりになったところから始めたいと思います。
(その４でtask_BH_trainまで終了したところから)

ベイズ最適化により求めたハイパーパラメータの結果
eta=0.103; colsample_bytree=0.957; min_child_weight=2; nrounds=147; max_depth=4

imlパッケージを使うために学習器のpredict.type(予測表示)を"prob"(確率)に変更します。

```{r}
classif.lrn = makeLearner("classif.xgboost", 
                          objective= "multi:softprob",
                          eta=0.103,
                          colsample_bytree=0.957,
                          min_child_weight=2,
                          nrounds=147,
                          max_depth=4,
                          predict.type = 'prob'                       
                         )
fit = train(classif.lrn, task_BH_train) 　
fit

$>Model for learner.id=classif.xgboost; learner.class=classif.xgboost
$>Trained on: task.id = BostonHousing; obs = 354; features = 14
$>Hyperparameters: $>nrounds=147,verbose=0,objective=multi:softprob,eta=0.103,colsample_bytree=0.957,min_child_weight=2,max_depth=4
```
##3.imlパッケージの適用
predictorを適用させます。データは説明変数と目的変数を分離させます。

```{r}
library(iml)
dataX=task_BH_train$env$data[,-13]　#目的変数の除去
dataY=as.integer(task_BH_train$env$data[,13]) #factorを整数に変換
predictor = Predictor$new(fit, data = dataX, y = dataY, type = "prob")
```

##4.Partial Dependence Plots (PDPs) 
PDPは興味のある変数以外の影響を周辺化して消してしまうことで、インプットとアウトカムの関係を単純化します。
ある特徴量に対するICE曲線は各観測値に対して特徴量のユニークな値を強制的に取らせたときの各観測値に対する予測値を示します。
PDP曲線はすべての観測値の平均予測値を示します。

特徴量"lstat"を例にPDPとICEを求めます。

```{r}
pdp = Partial$new(predictor, feature = "lstat")
pdp$plot()
```
![a.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/7c582673-f34d-be15-ff8c-ab64bddff9f7.png)
PDPは黄色、ICEは黒色です。
PDPは説明変数間に相関関係が強くなるとうまく働かないようです。

##5.Accumulated Local Effects (ALE) Plot
特徴量が機械学習モデルの予測に平均的にどのように影響するかを調べます。
ALEプロットは、部分依存プロット（PDP）に比べて、より高速で偏りのない代替方法です。
・*メリット* 
ALEはPDPと比べ、特徴量が相関している場合でも機能します。
ALEプロットの解釈は明確で、与えられた値に条件付きで、特徴量を変更した場合の予測に対する相対的な効果をALEプロットから読み取ることができます。
ALEプロットは、ゼロを中心にしています。これは、ALE曲線の各ポイントでの値が平均予測値との差だからです。

・*デメリット*
ALEプロットは、値の間隔が大きくなると少し不安定になることがあります（小さな浮き沈みが多い）。この場合、間隔を小さくすると、推定値がより安定しますが、間隔が小さすぎると、ALEプロットは正確性がそこなわれる可能性があります。
PDPとは異なり、ALEプロットにはICE曲線は付きません。ICE曲線が優れているのは、特徴効果の不均一性、ある特徴量の効果が、データの部分集合によっては異なって見えることが明らかになるからです。

```{r}
ale = FeatureEffect$new(predictor, feature = "lstat")
ale$plot()
```
![b.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/2f7c5eca-8733-c692-e850-b9cc31d5cd74.png)

##6.特徴量間の相互作用
ある特徴量が変化したとき、他の特徴量にどの程度影響を与えるのかを知ることができます。0が相互作用なし、1が100%作用していることを表します。

```{r}
interact = Interaction$new(predictor)
plot(interact)
```
![c.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/fd4ee1cf-7f79-cd28-b25d-df4a9cdb0b54.png)

また、ある特徴量を指定して、他の特徴量との双方向の作用を調べることもできます。


```{r}
interact = Interaction$new(predictor, feature = "rm")
plot(interact)
```
![d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/763f81b5-3840-2031-0ec1-3075f679591d.png)
##7.Shapley Value
Shapley　Valueを使って各特徴量の貢献度を計算します。
簡単に表現すると、予測値が得られる過程を特徴量同士の協力ゲームとみなし、協力の有無に応じて報酬（予測値）を分配するということです。
まずテストデータによる予測をします。予測結果を見てみます。

```{r}
testTask =subsetTask(task_BostonHousing ,subset =test.set )
pred = predict(fit, testTask)
pred$data[30:35,]

$>id truth       prob.2     prob.1       prob.0         response
$>79  30     1 0.0010077378 0.92985988 0.0691323802        1
$>88  31     1 0.0023708523 0.98703450 0.0105946120        1
$>95  32     1 0.0004544009 0.98309523 0.0164503735        1
$>97  33     1 0.0009404384 0.98599362 0.0130659556        1
$>99  34     2 0.9919853210 0.00780226 0.0002123969        2
$>100 35     1 0.0670449585 0.92902768 0.0039274190        1
```
id99がクラス３(response=2)と判断した根拠を見てみます。

```{r}
testTask =subsetTask(task_BostonHousing ,subset =test.set )  #（その４）からの続き
testx=testTask$env$data[,-13]
shapley = Shapley$new(predictor, x.interest = testx[34,])
shapley$plot()
```
![e.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/d83c935c-48b2-045b-38a1-0c7b7f9fb5e3.png)

```{r}
shapley$y.hat.average
$>         0          1          2 
$>0.40666482 0.51124503 0.08209015 

shapley$results
feature class           phi      phi.var feature.value
$>
$>29    crim     2 -1.826819e-03 2.210042e-03  crim=0.08187
$>30      zn     2  6.909858e-05 2.313188e-07          zn=0
$>31   indus     2  3.004585e-02 3.351996e-03    indus=2.89
$>32     nox     2  1.022677e-02 6.606550e-04     nox=0.445
$>33      rm     2  6.380180e-01 1.204818e-01       rm=7.82
$>34     age     2  1.216921e-02 3.808842e-03      age=36.9
$>35     dis     2 -1.483182e-02 6.444198e-03    dis=3.4952
$>36     rad     2 -1.831554e-03 7.526614e-05         rad=2
$>37     tax     2 -1.462168e-03 2.373234e-03       tax=276
$>38 ptratio     2 -1.541350e-03 5.509440e-04    ptratio=18
$>39       b     2 -1.875290e-04 2.351367e-04      b=393.53
$>40   lstat     2  2.267409e-01 7.286554e-02    lstat=3.57
$>41  chas.0     2  2.150447e-03 8.535962e-05      chas.0=1
$>42  chas.1     2  1.374929e-04 1.814638e-07      chas.1=0
```
表の見方ですが、クラス３のベースの予測値が8%(y.hat.average)になり、説明変数のrmが+64％、lstatが+23％とそれぞれ説明変数の寄与分(phi)の割合を足していきます。
rmとlstatだけで、すでに予測値は95%になっています。

#8.参考
[Rのimlで複雑なモデルの予測結果を簡潔に説明する -eXplainable AI, XAI-](https://qiita.com/gnbrganchan/items/de452da1f19666e8d132)
[変数重要度とPartial Dependence Plotで機械学習モデルを解釈する](https://dropout009.hatenablog.com/entry/2019/01/07/124214)
[Interpretable Machine Learning](https://christophm.github.io/interpretable-ml-book/index.html)