Rのimlで複雑なモデルの予測結果を簡潔に説明する -eXplainable AI, XAI-

#はじめに
近年のAI・機械学習ブームの一方で、一部の予測モデルのブラックボックス的な性質ゆえに、AIへの過剰な期待・無責任なAI利用への警鐘を鳴らす声も頻繁に聞かれるようになりました。
[AIの判断、企業に説明責任](https://www.nikkei.com/article/DGXMZO3818677026112018MM8000/)という日経新聞の記事が話題になりましたが、特に経営に関わる判断にAIを活用する際にはその説明能力を問われるようになってきています。
このような、ただ予測精度を求めるだけでなく、「なぜそのような判断をしたのか」を説明できるAIは一般的に説明可能なAI（eXplainable AI, XAI）と呼ばれ、今、多くの企業から注目を集めています。
今回の記事では、Rのiml (Interpretable Machine Learning) パッケージを使って様々な角度からブラックボックスモデルの中身を説明する方法をご紹介していきます。
本記事で紹介するソースコードは[github上で公開しています](https://github.com/gnbrganchan/Qiita/tree/master/iml)。

#imlパッケージの概要
imlは2019年に公開されたRパッケージ[CRAN - Package iml](https://cran.r-project.org/web/packages/iml/index.html)です。
Fisher et al. (2018)のFeature importanceや、 Apley (2018) のaccumulated local effects plotsといった、説明可能なAIに関する技術が実装されています。

##imlでできること
imlでできることを大まかに分類すると下記の通りになります。
1. モデルの性能や特性を評価する
 - Feature Importance
 - Tree Surrogate
2. 特徴量（変数）に対するモデルの応答をみる
 - Feature Effect (ALE)
 - Feature Interaction
3. データに対する予測がどのように得られたかを説明する
 - LIME
 - Shapley Value
以下でそれぞれについて提供されているメソッドを紹介します。

#imlを使ってみる
##ライブラリの呼出し
ライブラリ名はそのまま"iml"です。
また、ブラックボックスAIとしてよく使われるランダムフォレストのライブラリも呼び出しておきます。

```R:imlを呼出し
install.packages("iml") #初回読み込み前に要install
library(iml)
library(randomForest)
```
##使用するデータ
この記事では詳しくは触れませんが、MASSのBostonという有名なデータセット（Boston市の住宅価格データセット）を使用します。
人口1人あたりの犯罪発生率 (crim)や住居の平均部屋数 (rm)などの特徴量（変数）から住宅価格 (medv) を予測するのが目的となります。

```R:Bostonデータセットを読込
data("Boston", package  = "MASS")
head(Boston)
#     crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
#1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
#2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
#3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
#4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
#5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
#6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21
#  medv
#1 24.0
#2 21.6
#3 34.7
#4 33.4
#5 36.2
#6 28.7
```
##ランダムフォレストモデルの作成
ランダムフォレストモデルの作成は目的変数、使用データセット、木の深さを指定するだけで簡単にできます。
imlのPredictor\$new()を呼び出すことで、ランダムフォレストのモデルをデータを保持したオブジェクトを作成することができます。
以降、代表的なブラックボックスモデルであるこのランダムフォレストモデルを、imlのメソッドを使って様々な角度から紐解いていきます。

```R:RandomForestモデル⇒Predictorオブジェクト作成
rf = randomForest(medv ~ ., data = Boston, ntree = 50)
X = Boston[which(names(Boston) != "medv")]
predictor = Predictor$new(rf, data = X, y = Boston$medv)
```

##1. モデルの性能や特性を評価する
###FeatureImp
FeatureImp\$new()を呼び出すことで、どの特徴量（変数）が重要であるかを測ることができます。
Feature importanceは、それぞれの特徴量（変数）をシャッフルしたときにどのくらい精度が低下するかを調べることで、その変数の重要性を算出します。
例えば、ある変数列の値をシャッフルすると予測精度が大きく低下した場合、その変数は重要度が高いと言えます。
FeatureImp\$new()の結果をplotすることで可視化できますし、\$resultsによってdata.frame形式で出力することもできます。

```R:FeatureImportance
imp = FeatureImp$new(predictor, loss = "mse")
plot(imp)
```

![featureimportance.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/458253/a6922feb-60ff-2c02-176a-ea93571af17b.png)


```R:FeatureImportance$results
imp$results
#   feature importance.05 importance importance.95 permutation.error
#1    lstat    15.1501353  22.074149     27.975032         45.377962
#2       rm    10.9331383  11.489301     19.441368         23.618625
#3      dis     2.6690680   3.027744      4.525194          6.224152
#4      nox     2.8722394   3.003836      3.498187          6.175004
#5     crim     2.4794144   2.856576      2.971373          5.872280
#6    indus     1.8486110   2.099748      2.274380          4.316464
#7  ptratio     1.8601196   2.006714      2.251311          4.125215
#8      age     1.5103398   1.827092      2.179101          3.755965
#9      tax     1.5761642   1.771273      2.111650          3.641216
#10   black     1.3858946   1.471415      1.700118          3.024797
#11     rad     1.0344114   1.298408      1.396144          2.669146
#12    chas     1.0462975   1.126355      1.171892          2.315455
#13      zn     0.9673134   1.055457      1.219117          2.169709
```
Bostonデータセットを使った上記の例では、lstat（給与の低い職業に従事する人口の割合[%]）が住宅価格への一番影響力が大きいという結果が得られました。

###TreeSurrogate
TreeSurrogate\$new()は、ランダムフォレストのような複雑なモデルの入力と予測値から、よりシンプルなモデル（決定木）を生成します。
すなわち、ブラックボックス化したAIを人間が見ても理解できる形に簡易化して可視化してくれます。
maxdepthパラメータにより、生成する決定木の深さを指定します。
決定木が深いほどモデルとしては複雑になりますが、より元のモデルを正確に表現します。

```R:TreeSurrogate
tree = TreeSurrogate$new(predictor, maxdepth = 2)
plot(tree)
```
![treesurrogate.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/458253/41b6cf96-bf6d-ae0a-5deb-a3328ebc3081.png)

この例ではランダムフォレストを決定木に変換していますが、あらゆるブラックボックスモデルに適用可能です。
ここで注意しておかなければならないのは、これはブラックボックスモデルの構造を解釈しているわけでなく、単にブラックボックスモデルへの入力と出力のデータから決定木を構成しているということです。

ここで得られた決定木を使って予測を行うこともできます。

```R:TreeSurrogate-Predict
head(tree$predict(Boston))
#    .y.hat
#1 27.09989
#2 27.09989
#3 27.09989
#4 27.09989
#5 27.09989
#6 27.09989
# 警告メッセージ: 
# self$predictor$data$match_cols(data.frame(newdata)) で: 
#  Dropping additional columns: medv
```

##2. 特徴量（変数）に対するモデルの応答をみる
###FeatureEffect (ALE)
ALE (Accumulated Local Effect)を計算することで、ある特徴量（変数）の近傍で、平均的にモデルの予測値がどう変化するか？ ということを知ることができます。
imlではFeatureEffect\$new()として定義されています。

```R:FeatureEffect-'lstat'
ale = FeatureEffect$new(predictor, feature = "lstat")
ale$plot()
```
![featureeffect.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/458253/dba988c3-2f3c-36c0-ca3d-b53f58fb7450.png)

X軸のマークは'lstat'の分布を表しており、どの領域が関係性が高いかを示しています。
つまり、点が少ない領域にはあまり関心を払う必要がないことを示唆しています。

下記のように、注目する特徴量を切り替えることもできます。

```R:FeatureEffect-'rm'
ale$set.feature("rm")
ale$plot()
```
![featureeffect2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/458253/a576e1a3-babd-a80c-0c9d-8b25ffd4e7ce.png)


ALEと似た概念にPDP (Partial Dependence Plot)がありますが、PDPでは「あるデータに対して、注目する特徴量（変数）以外はすべて同じ値をもつデータがたくさんあったとき、対象の特徴量（変数）の値の増減によって予測値がどう変化するか？」ということを調べます。
ALEはPDPと比べて、
 - 特徴量（変数）同士が相関してもうまく動く
 - 必要な計算量が少ない
といった特徴があります。

###FeatureInteraction
また、Interaction\$new()を使うことで、特徴量（変数）間の相互作用を計算することもできます。
つまり、ある特徴量（変数）が変化したとき、他の特徴量（変数）にどの程度影響を与えるのかを知ることができます。
これはH-statisticを指標としており、0が相互作用なし、1が100%作用していることを表します。

```R:Interaction
interact = Interaction$new(predictor)
plot(interact)
```
![interaction.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/458253/b4b7366a-65b9-8cf0-c750-586ef791edfb.png)

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

```R:Interaction-2way
interact = Interaction$new(predictor, feature = "crim")
plot(interact)
```
![interaction2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/458253/0f581259-7eeb-5935-f39b-8fd08d3a6e33.png)

H-statisticの詳しい説明はここでは割愛しますが、H-statisticは[0,1]の無次元量であるため、異なる特徴量・モデル間の比較に使うことができます。


##3. データに対する予測がどのように得られたかを説明する
###LIME (Locally Interpretable Model-agnostic Explanations)
TreeSurrogateは入力値全体に対するモデルの挙動をシンプル化してくれるものでしたが、LIMEはモデルの局所的な挙動をシンプルに説明してくれます。
つまり、説明したいデータの周辺で線形近似を行い、近似モデルの重みによって、各変数の予測に対する説明を表します。
LIMEについては[LIMEで機械学習の予測結果を解釈してみる](https://qiita.com/fufufukakaka/items/d0081cd38251d22ffebf)に詳しい解説が書かれていますので、この記事では割愛します。

```R:LIME
library(glmnet)
library(Matrix)
library(foreach)
library(gower)
lime.explain = LocalModel$new(predictor,x.interest = X[1,])
lime.explain$results
#              beta x.recoded    effect x.original feature feature.value
#rm       4.4836483     6.575 29.479987      6.575      rm      rm=6.575
#ptratio -0.5244767    15.300 -8.024493       15.3 ptratio  ptratio=15.3
#lstat   -0.4348698     4.980 -2.165652       4.98   lstat    lstat=4.98
plot(lime.explain)
```
![lime.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/458253/c4f0ddcb-765b-aa9d-3752-2eb3899aed4c.png)

近似が単純な分、得られる説明が短く分かりやすいのが利点ですが、ブラックボックスモデルの複雑性に依存して出力が安定せず、サンプルプロセスを繰り返すと出力が変わる可能性がある点には注意が必要です。

###Shapley
ゲーム理論を使って各変数の貢献度を計算しよう、という少し今までとはテイストの違う考え方がShapley Valueです。
Shapleyの解説をするとそれだけで1記事できてしまうので割愛します（後日別記事で解説しようと思います）が、簡単に言えば、予測値が得られる過程を特徴量（変数）同士の協力ゲームとして、協力の有無に応じて報酬（予測値）を分配するということです。
Shapleyについては[SHAPでモデルの予測結果を説明する](https://orizuru.io/blog/machine-learning/shap/)に書かれている説明がわかりやすくて良いと思いました。
imlではShapley\$new()で計算できます。
すべてのデータとすべての特徴量（変数）の組み合わせを計算すると膨大な時間がかかっているため、モンテカルロ計算をしています。
sample.sizeパラメータでサンプリング回数を指定でき、サンプル回数が多いほど結果は安定します。（必要な計算量は多くなります）
LIMEと比べて、一部のデータに対する説明ができることが利点です。

```R:Shapley
shapley = Shapley$new(predictor, x.interest = X[1,])
shapley$plot()
```
![shapley.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/458253/acb06114-e0d9-d0cc-63b0-b98926488a00.png)

Shapley Valueをdata.frame形式で出力することもできます。

```R:Shapley-data.frame
results = shapley$results
head(results)
#  feature          phi      phi.var feature.value
#1    crim -0.317960865  2.046510169  crim=0.00632
#2      zn -0.036996333  0.024963557         zn=18
#3   indus  0.964100389  0.558548012    indus=2.31
#4    chas -0.006707333  0.002755066        chas=0
#5     nox -0.151510825  1.830613219     nox=0.538
#6      rm -0.450538504 15.837141138      rm=6.575
```

#おわりに
今後も説明可能なAIの需要は高まっていき、次々と新しい理論・実装が現れてくることが予想されます。
これまでもXAIを試みたパッケージはいくつか公開されていますが、今回ご紹介したimlは現時点での最先端と呼べるパッケージだと思います。
まだ公開されてから日の浅いパッケージですので、今後の更なる発展に期待しつつ、今あるメソッドを存分に使いこなしましょう。

本記事も随時アップデートしていきます。
本文中の誤記・不適切な記述にお気づきになりましたらお手数ですがコメントにてお知らせいただけると幸いです。

###参考
[Introduction to iml: Interpretable Machine Learning in R](https://cran.rstudio.com/web/packages/iml/vignettes/intro.html)
[How to use in R model-agnostic data explanation with DALEX & iml](https://www.slideshare.net/kato_kohaku/how-to-use-in-r-modelagnostic-data-explanation-with-dalex-iml)
[Wikipedia -Shapley値](https://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%A3%E3%83%BC%E3%83%97%E3%83%AC%E3%82%A4%E5%80%A4)