琥珀色呑んだくれ備忘録

メモとか備忘録とか

XGBoostExplainerが何をやっているか調べる(1.とりあえず使う)

目的

XGBoostの予測を分解するツールXGBoostExplainerは、あるインスタンスについて得られたXGBoostによる予測結果が、どのように構成されているか可視化してくれる。

f:id:kato-satoshi-0:20181212125239p:plain

コンセプトとしては、randomforestにおけるforestfloorと同じく、feature contributionを算出する。ターゲット集団の変数の感度分析に使うのではなく*1、個々のインスタンスのxgboostによる予測結果の説明をWaterfall Chartで可視化してくれるのが特徴。

なお、同様のWaterfall Chartによる説明は、そのほかの手法についても提供されているようだ。

lightgbmの予測結果に適用してくれるパッケージ

lm/glmの予測結果に適用してくれるパッケージ

探索的なデータ分析に使うだけならおおむね問題ないと思うが、具体的に何をやっているか説明しようとして困った。論文などの詳細な資料も見つけられなかったため、XGBoostExplainerの実装を追いかけて、何をやっているか調べた。

関連シリーズ

  1. とりあえず使ってみる(この記事)
  2. 予測結果の可視化プロセスをstep-by-stepで実行する
  3. 学習したxgboostのルール抽出をstep-by-stepで実行する
  4. 予測結果のbreakdownをstep-by-stepで実行する

参考

開発元の紹介記事

medium.com

とりあえず使ってみる

すでに日本語の記事がある。

qiita.com

xgboostExplainerのマニュアルにあるexampleからコピペを眺める。

インストール

本家の記事に従ってgithubからインストール

install.packages("devtools") 
library(devtools) 
install_github("AppliedDataSciencePartners/xgboostExplainer")

XGBモデルの学習と予測

今回はxgboostパッケージ付属のサンプルデータで、定番の食えるキノコと毒キノコの2値分類。細かいチューニングは、必要に応じてautoxgbあたりでチューニングするとよいが、今回は省略。

library(xgboost)
require(tidyverse)
library(xgboost)
library(xgboostExplainer)

set.seed(123)

data(agaricus.train, package='xgboost')

X = as.matrix(agaricus.train$data)
y = agaricus.train$label
table(y)
train_idx = 1:5000

train.data = X[train_idx,]
test.data = X[-train_idx,]

xgb.train.data <- xgb.DMatrix(train.data, label = y[train_idx])
xgb.test.data <- xgb.DMatrix(test.data)

param <- list(objective = "binary:logistic")
xgb.model <- xgboost(param =param,  data = xgb.train.data, nrounds=3)
# 
# col_names = colnames(X)
# 
# pred.train = predict(xgb.model,X)
# nodes.train = predict(xgb.model,X,predleaf =TRUE)
# trees = xgb.model.dt.tree(col_names, model = xgb.model)

個別の予測結果の可視化

xgboostExplainerのマニュアルにあるexampleのコピペ(つづき)。高度にwrapされているためわずか3行でstep-by-stepが完了する。

STEP.1. 学習済みXGBモデルからルールセット(leafまでのパス)を列挙してテーブル化

base_scoreオプションはxgboostのオプションそのままで、ターゲット集団のクラス不均衡を表す事前確率。すなわち正例:負例=300:700を仮定できる対象であれば、base_score = 0.3となる(デフォルトは1:1を表す0.5)。

library(xgboostExplainer)

explainer = buildExplainer(xgb.model,xgb.train.data, type="binary", base_score = 0.5, trees = NULL)
#> 
#> Creating the trees of the xgboost model...
#> Getting the leaf nodes for the training set observations...
#> Building the Explainer...
#> STEP 1 of 2
#> 
#> Recalculating the cover for each non-leaf... 
#> 
  |=================================================================| 100%
#> 
#> Finding the stats for the xgboost trees...
#>                                                                        
  |=================================================================| 100%
#> 
#> STEP 2 of 2
#> 
#> Getting breakdown for each leaf of each tree...
#>                                                                 
  |=================================================================| 100%
#> 
#> DONE!

STEP.2. Get multiple prediction breakdowns from a trained xgboost model

マニュアルには step2とあるのだが、実はパッケージを使うだけならスキップできてしまう。

STEP.3. 予測対象(インスタンス)に適用される各treeのパスを集計して可視化

2値分類(binary:logistic)では、片側のクラスに属する確率p(左軸の数値)のロジット(対数オッズ;棒グラフ中の数値)が足し合わされている様子を表示する。

showWaterfall(xgb.model, explainer, xgb.test.data, test.data,  2, type = "binary")
#> 
#> 
#> Extracting the breakdown of each prediction...
#> 
 |=================================================================| 100%
#> 
#> DONE!
#> 
#> Prediction:  0.811208
#> Weight:  1.457879
#> Breakdown
#>        intercept  gill-size=broad        odor=foul        odor=none 
#>      -0.27084657       1.61423045      -0.67129347       0.45408751 
#>       odor=anise      odor=almond cap-color=yellow 
#>       0.13628094       0.13073006       0.06468987

f:id:kato-satoshi-0:20181212131040p:plain

(参考) binary:logisticの場合、base_scoreで設定した事前確率は ベースラインとしてinterceptに反映される。下記の例ではinterceptだけが下がっていることに注目されたい。

explainer = buildExplainer(xgb.model,xgb.train.data, type="binary", base_score = 0.2, trees = NULL)
showWaterfall(xgb.model, explainer, xgb.test.data, test.data,  2, type = "binary")
#> 
#> 
#> Extracting the breakdown of each prediction...
#> 
  |=================================================================| 100%
#> 
#> DONE!
#> 
#> Prediction:  0.5178885
#> Weight:  0.07158444
#> Breakdown
#>        intercept  gill-size=broad        odor=foul        odor=none 
#>      -1.65714093       1.61423045      -0.67129347       0.45408751 
#>       odor=anise      odor=almond cap-color=yellow 
#>       0.13628094       0.13073006       0.06468987

f:id:kato-satoshi-0:20181212131004p:plain

次回は、xgboostExplainerにより、xgboostのモデルと予測結果から何が取り出され、どう捌かれているかだけを詳細に見ていく。

*1:XGBoostExplainerの開発者による記事では、感度分析も行っている

[感想]『自然科学研究のためのR入門』

まずは、本書を御恵贈いただき、あらためて著者に感謝したい。あくまで感想文として、この本を教科書として、同僚(特に若い方)やカウンターパートと仕事を共有するうえで、自分ならどのようなポイントに注意するかに焦点を絞ってコメントしたい*1

全体

対象としない読者

  • ご自身ではデータ分析とレポート作成をしないポジションのかた
  • このご時世でも『研究の再現性は不要』と思ってるひと
  • 『再現可能なレポート』について精通しているRユーザー
  • Rによるデータ分析自体の習得を目指しているかた

上記以外は必読。ただし、研究責任者としてご自身ではデータ分析に従事しない立場であっても、こうした技術的に解決できる枠組みを検討せず、対策を講じていないままという態度は無責任といえよう。

www.kyoritsu-pub.co.jp

測定した、シミュレーションした、あるいはカウンターパートから預かったデータを分析するのが、(特に自然科学の)研究者の日常的な活動のひとつ*2

言い尽くされていることではあるが、報告した結果を元のデータから再現できることは、データ分析者にとっての最良の自衛手段であり義務でもある。その立場と目的にたって、本書は研究者や学生が『再現可能なレポート』を書くために、必要なステップと今の技術を最速で習得するための一冊になっている *3

COI

著書をご恵投いただく程度には交流がある。

1章:RmarkdownとRstudioを使う

Rによるデータ分析結果を報告するという目的上、報告以前のデータ分析の知識はある程度は持っていることが前提。

引用文献リストの自動取り込み方法も触れられている。引用文献リストの掲載がひとつもないレポートは研究報告には値しないので*4*5、この作業の定型化と省力化は研究報告の作成時には非常に重要なポイントとなる。類書ではあまり最初に出てこないが、本書ではこれが1章でカバーされており、ツールの選定・検討中の研究者には嬉しい。

本書の特色として、各章にsession_info()の出力が付けられている。報告相手が別の環境でレポートの再現を試みることを想定して、分析環境の情報共有も重要なポイントといえる。 自分以外の人に、結果(あるいはデータ分析がうまく行かない)を相談したり、報告内容を再現してもらう際にしばしばトラップとなる*6

2,3,4章:統計モデリング、実験計画法と分散分析を使った分析とレポーティング

これらの各章では、データの分析をstep-by-stepで実施しながら、一通りの分析が終わった時点でレポートが出来上がっているという流れになっており、Rmarkdownの嬉しさが最大限に発揮されるポイントのひとつといえるだろう。

これ以降の章では、一連の分析風景を追いかけてくれているので具体的なイメージを持ちやすい。 機械学習などによる分析や自身で提案した新しい分析手法を報告する際にも、既存の方法をベースラインとした対比が必要となる*7。 データ分析に際して、どのような順番でどのような方法でアプローチしていくか、あるいはそれぞれの分析において着目すべきポイントや勘所が随所にちりばめられており、著者の経験に基づいた実践的な知識の一端に触れることができるのも楽しい。

5章:機械学習(など)を使った分析のレポーティング

冒頭でごく簡単に、Bioconductorからのパッケージの導入にも触れている。本書の読者には、実験から解析までこなすバイオインフォマティシャンも多く想定され、分析風景がより身近になるのではないだろうか。

再現性の観点から、欲を言えば、機械学習の方法には乱数を利用するものも多いため*8set.seed()で乱数を固定する方法は紹介されても良かったかなと思う *9 *10

最近は、初期検討として機械学習の各種手法を横並びに比較をすることも多く、メタフレームワークのひとつであるcaretを使った学習パラメタのチューニング・性能評価と比較方法の紹介はありがたい *11 *12

また、xgboostのパラメータチューニングについては、全自動でやってくれるautoxgboostもあるので薦めておきたい。

kato-kohaku-0.hatenablog.com

AUCROCによるパフォーマンスの比較について

ROCが三角形になっていて気になった*13。これは予測ラベルを直接投入しているためと思われるが、各症例の予測値や検査値をもとにROCをみながらカットオフの検討をしてる業界などでは、この図で報告すると戸惑われる可能性がある *14

たとえば、random forestのラベルの得票率とlassoの回帰の値をそれぞれ投入して、両者を比較したらよかったのではないかと思う *15

require(randomForest)
require(glmnet)
require(pROC)

data(spam, package = "kernlab")
# glimpse(spam)

X <- spam[, -NCOL(spam)]
y <- as.factor(spam[,  NCOL(spam)])

train <- sample(1:NROW(spam), NROW(spam) * 0.7)
test  <- setdiff(1:NROW(spam), train)

# train learner + predict test with random forest
rf <- randomForest(X[train, ], y[train], ntree=500)
pred.rf <- predict(rf, X[test, ], type="prob")
head(pred.rf)
#>    nonspam  spam
#> 4    0.026 0.974
#> 24   0.296 0.704
#> 30   0.082 0.918
#> 32   0.078 0.922
#> 33   0.236 0.764
#> 34   0.844 0.156

# tune-train learner + predict test with lasso
cv <- cv.glmnet(x = as.matrix(X[train, ]), y = as.numeric(y[train]))
lasso <- glmnet(x = as.matrix(X[train, ]), y = as.numeric(y[train]), lambda = cv$lambda.1se, alpha = 1)
pred.lasso <- predict(lasso, as.matrix(X[test, ]))
head(pred.lasso)
#>          s0
#> 4  1.481061
#> 24 1.529340
#> 30 1.899370
#> 32 1.652855
#> 33 1.592823
#> 34 1.532178

# prot ROC
roc.rf <- roc(response = y[test],             predictor = pred.rf[,1])
roc.la <- roc(response = as.numeric(y[test]), predictor = c(pred.lasso))

mt <- sprintf("RF(red) vs lasso(blue) = %.3f vs =%.3f",auc(roc.rf),auc(roc.la))
plot.roc(roc.rf, col = "red",  main = mt)
plot.roc(roc.la, col = "blue", add = TRUE)

f:id:kato-satoshi-0:20181013115315p:plain

6章:集大成

この章では、レポートの提出を強く意識する*16 報告相手に次第では、コード部分や処理途中の出力などはノイズでしかなく、場合によってはレポートそのものを見てもらえないこともある *17。重要なポイントと思った。

結論:(Rで)データ分析をしている研究者には必携の書

  • 報告の必要がなくても、手元の分析結果をレポートの形で残す習慣を身に着けたい
  • R以外で分析をするユーザーを対象とした書籍が出てくれると良い。

www.kyoritsu-pub.co.jp

*1:献本を遣り取りする文化から転職して既に長いため、書評のようにきちんと精査したものではないことは、ご容赦されたい。

*2:純粋な理論系は違うかもしれないが..

*3:必要最低限の情報提供にとどめることで、常々忙しい忙しいとか多言する研究者に取捨選択をさせないスタイルになっているが、各章末に発展内容への適切なガイドを付けることでカバーされている

*4:少なくとも自然科学の研究者にとって

*5:ただの作業記録やポエムならともかく

*6:報告書の作成時には見落としがちであるが、OSが違ったり、分析に必要なパッケージが入ってなかったり、使っているパッケージがバージョンアップでオプションの指定方法が変わったり、出力されるデータの構造が変わったりすることがある

*7:情報学系なら事情が異なるかもしれないが、自分の業界では査読で必ず突っ込まれる

*8:例えば本書でも取り上げられているrandom forestなど

*9:特に初学者が再現を試みたときに結果の数字が元のテキストと違っていると戸惑うことが多い

*10:一意な乱数に固定することで、報告対象にチャンピオンデータを誘導したりネガティブな結果を隠蔽しかねない悪用や弊害があるのは承知している

*11:自分はmlRユーザーなのでcaretはあまり常用してないけど、事情はそれほど変わらないだろう

*12:ついでながら、expand.grid()は、caretじゃなくてbaseの関数だけど、実践的には困ることはないだろう

*13:情報系の論文では見るので、目的に応じるのだと思う

*14:誰に報告するかといった業界や事情の違いがあるので、あくまで個人的な体験に基づく

*15:手元のデータでパパっとやったので、本書のデータでの解析例ではないことはご容赦されたい

*16:include=FALSE や echo=FALSE といったチャンクオプションを使いながら、適宜、隠蔽する方法が示されており、分析者が手元で再現できれば良い、あるいは、付録として分けるといった説明がされている。

*17:このあたりは、どういう相手かに応じて分ける必要がある

Windows 10 で GPU(CUDA)を利用するLightGBM のR-packageをセットアップする

基本的には公式サイトのガイドに従ってセットアップするだけ。だけなのだが、ちょこちょこ入れ忘れとか落とし穴があったりして、次から手間取らないように備忘録にしておく。

全体の準備

32bit 版のRをインストールしない

わりと嵌まったポイント。LightGBMは64bit版しかサポートしないが、32bit 版のRが入っているとビルドの際に32bit版をビルドしようとして失敗するとのことで、解決方法は、Rのインストール時に32bit版を入れないようにする(ホントにそう書いてある)。

f:id:kato-satoshi-0:20180828234534p:plain

Install CMake

公式サイト からダウンロードアイコンをクリック → インストーラーのダウンロード → セットアップ

f:id:kato-satoshi-0:20180828234807p:plain

Install Rtools

公式サイト からRtools35.exeをダウンロードしてインストール。

f:id:kato-satoshi-0:20180828234841p:plain

環境変数Pathに、Rtools\binRtools\mingw_64\binを設定しておく。 f:id:kato-satoshi-0:20180828234915p:plain

Install Visual Studio

公式サイト からダウンロード。

インストールは、とりあえずC++だけでよい。

f:id:kato-satoshi-0:20180828235059p:plain

Install MS-API

MS-MPI SDK と setup.exe を両方公式サイトからダウンロード

msmpisdk.msimsmpisetup.exe を実行してインストール。

f:id:kato-satoshi-0:20180828235128p:plain

Install git

公式サイト からダウンロード。

設定は適当に。

Install OpenCL

インストールガイドの

に従って、今回はCUDA Toolkitをインストール。

f:id:kato-satoshi-0:20180828235410p:plain

今回は、特にカスタマイズしなくてもよかった。

Install Boost Binary

sourceforgeからダウンロード。

Visual Studio 2017 をインストールしたので、 msvc-14.1-64.exeをインストールする。

f:id:kato-satoshi-0:20180828235504p:plain

Boostの環境設定

BOOST_ROOTBOOST_LIBRALYDIRのパスの設定を行う

f:id:kato-satoshi-0:20180828235634p:plain f:id:kato-satoshi-0:20180828235640p:plain

環境変数の確認

「コントロール パネル > システムとセキュリティ > システム(Win+Pause/Break)」から「システムの詳細設定→環境変数」をチェック

f:id:kato-satoshi-0:20180829000335p:plain

以下の項目を確認:

  • CMake関連のパスがPath変数に追加されている
  • Rtools関連のパスがPath変数に追加されている
  • MSAPI関連の環境変数がが登録されている
  • openCL(CUDA)関連の環境変数が追加されている
  • boost関連の環境変数が追加されている

f:id:kato-satoshi-0:20180829000352p:plain

Windows を再起動する

超重要。忘れる。

LightGBMのインストール

lgbdlを使ってインストールする。

lgbdlのインストール

devtools::install_github("Laurae2/lgbdl")
## Downloading GitHub repo Laurae2/lgbdl@master
## from URL https://api.github.com/repos/Laurae2/lgbdl/zipball/master
 
## Installing lgbdl
## "C:/PROGRA~1/R/R-35~1.1/bin/x64/R" --no-site-file --no-environ --no-save  \
##   --no-restore --quiet CMD INSTALL  \
##   "C:/Users/MyName/AppData/Local/Temp/Rtmp0SHyaL/devtools43c0670829ab/Laurae2-lgbdl-afff103"  \
##   --library="C:/Users/MyName/Documents/R/win-library/3.5"  \
##   --install-tests --no-multiarch

インストール

use_gpu = TRUEを指定する。R 3.5以降のバージョンではR35 = TRUEをセットするとご利益があるらしい(ビルドの際の最適化に関する警告が違う)。

lgbdl::lgb.dl(commit = "master",
              compiler = "vs", # Remove this for MinGW + GPU installation
              repo = "https://github.com/Microsoft/LightGBM",
              cores = 4,
              R35 = TRUE,
              use_gpu = TRUE)
## Installing package into 'C:/Users/MyName/Documents/R/win-library/3.5'
## (as 'lib' is unspecified)

インストール中のメッセージが大量に表示される。

## [1] TRUE

当初、LoadLibrary failure: %1 is not a valid Win32 application.というエラーが出て、インストールに失敗したので調べたら、Rのインストールの際に32bit版を入れてしまうと、両方でビルドしようとしてこける(64bit版しかサポートしない)とのこと。

このほか、ビルド中にどこかで止まる場合、チェックポイントのどれかが不十分と思われるので、再確認。

動作確認

library(lightgbm)
## Loading required package: R6
data(agaricus.train, package='lightgbm')
train <- agaricus.train
dtrain <- lgb.Dataset(train$data, label=train$label)
## Loading required package: Matrix
params <- list(objective="regression", metric="l2")
model <- lgb.cv(params, dtrain, 10, nfold=5, min_data=1, learning_rate=1, early_stopping_rounds=10)
## [1]:  valid's l2:1.43524e-16+9.92589e-17 
## [2]: valid's l2:6.30955e-32+5.75747e-32 
## [3]: valid's l2:7.50672e-35+1.33422e-34 
## [4]: valid's l2:7.03798e-35+1.35156e-34 
## [5]: valid's l2:6.93275e-35+1.35631e-34 
## [6]: valid's l2:6.91914e-35+1.35695e-34 
## [7]: valid's l2:6.91914e-35+1.35695e-34 
## [8]: valid's l2:6.91914e-35+1.35695e-34 
## [9]: valid's l2:6.91914e-35+1.35695e-34 
## [10]:    valid's l2:6.91914e-35+1.35695e-34

(おまけ)install_github()を使って、非GPU版のインストール

GPU使わない場合はOpenCLBoostの準備が不要なので、こちらの方が楽。

library(devtools)
options(devtools.install.args = "--no-multiarch") # if you have 64-bit R only, you can skip this
install_github("Microsoft/LightGBM", subdir = "R-package")
## Downloading GitHub repo Microsoft/LightGBM@master
## from URL https://api.github.com/repos/Microsoft/LightGBM/zipball/master

## Warning: GitHub repo contains submodules, may not function as expected!

## Installing lightgbm

## "C:/PROGRA~1/R/R-35~1.1/bin/x64/R" --no-site-file --no-environ --no-save  \
##   --no-restore --quiet CMD INSTALL  \
##   "C:/Users/MyName/AppData/Local/Temp/Rtmp0SHyaL/devtools43c039486000/Microsoft-LightGBM-aac8e8f/R-package"  \
##   --library="C:/Users/MyName/Documents/R/win-library/3.5"  \
##   --install-tests --no-multiarch
## 

本当は、ビルド中にシステムメッセージが大量に表示される 。

autoxgboost を使ってみる

前回のTokyoRで@hoxo-mさんがつぶやいていたautoxgboostを使ってみる。

何?

autoxgboostは、mlrmlrMBOを使ってxgboostのモデルをベイズ最適化でチューニングするためのラッパー。気合い入れて検討する前にベースラインとる用途だけでも、かなり楽になるのではないかと思う。

(参考)

準備

インストールは github から。

install.packages("devtools") # if you have not installed "devtools" package
devtools::install_github("ja-thomas/autoxgboost")

今回は spam データセットを使う。

set.seed(0)
data(spam, package = "kernlab")

70%を教師データとして使う。このデータのラベルはtype列。

train.set <- sample(1:NROW(spam), NROW(spam) * 0.7)
table(spam$type[train.set])
#> 
#> nonspam    spam 
#>    1954    1266

残りをテストデータに使う。

test.set <- setdiff(1:NROW(spam), train.set)
table(spam$type[test.set])
#> 
#> nonspam    spam 
#>     834     547

autoxgboost のロード。必要な周辺パッケージも併せてロードされる。

require(autoxgboost)
#> Loading required package: autoxgboost
#> Loading required package: ParamHelpers
#> Warning: replacing previous import 'BBmisc::isFALSE' by
#> 'backports::isFALSE' when loading 'ParamHelpers'
#> Loading required package: mlr
#> Loading required package: mlrMBO
#> Loading required package: smoof
#> Loading required package: BBmisc
#> Loading required package: checkmate
#> Loading required package: mlrCPO

モデルのチューニング

一度、教師データでチューニングしてから、テストデータで評価することにする。ほぼマニュアル通りに書き下し。MBOcontrol()の細かな設定は、後日改めて確認しておきたい。

tune.task <- makeClassifTask(data = spam[train.set, ], target = "type")

ctrl <- makeMBOControl()
ctrl <- setMBOControlTermination(ctrl, iters = 1L) #Speed up Tuning by only doing 1 iteration

チューニングの実行。これもほぼマニュアル通り。mlrの分類タスクは デフォルトの評価指標がthe mean misclassification error (mmce) なので、これをaucに変更する。このあたりはタスクにあった指標をお好みで選べばよい。

design.sizeとかtune.thresholdとか、まだ理解していない引数もあるので、これもおいおい確認。

res <- autoxgboost(tune.task, control = ctrl, measure = auc, nthread = 4, tune.threshold = FALSE)
#> Computing y column(s) for design. Not provided.
#> [mbo] 0: eta=0.0992; gamma=18.4; max_depth=19; colsample_bytree=0.521; colsample_bylevel=0.529; lambda=5.39; alpha=35.8; subsample=0.733; scale_pos_weight=0.00616 : y = 0.968 : 0.7 secs : initdesign
#> [mbo] 0: eta=0.0715; gamma=0.109; max_depth=9; colsample_bytree=0.641; colsample_bylevel=0.952; lambda=1.46; alpha=480; subsample=0.865; scale_pos_weight=0.173 : y = 0.5 : 0.3 secs : initdesign
#> [mbo] 0: eta=0.181; gamma=25.6; max_depth=13; colsample_bytree=0.928; colsample_bylevel=0.618; lambda=10.6; alpha=7.48; subsample=0.569; scale_pos_weight=0.606 : y = 0.97 : 0.8 secs : initdesign
#> [mbo] 0: eta=0.138; gamma=0.0292; max_depth=16; colsample_bytree=0.728; colsample_bylevel=0.596; lambda=452; alpha=0.343; subsample=0.618; scale_pos_weight=0.00176 : y = 0.969 : 1.2 secs : initdesign
#> [mbo] 0: eta=0.0208; gamma=2.48; max_depth=4; colsample_bytree=0.744; colsample_bylevel=0.752; lambda=0.0143; alpha=0.859; subsample=0.635; scale_pos_weight=0.0286 : y = 0.965 : 0.4 secs : initdesign
#> [mbo] 0: eta=0.0567; gamma=36.3; max_depth=15; colsample_bytree=0.542; colsample_bylevel=0.769; lambda=30.5; alpha=153; subsample=0.979; scale_pos_weight=396 : y = 0.944 : 0.6 secs : initdesign
#> [mbo] 0: eta=0.13; gamma=5.39; max_depth=17; colsample_bytree=0.785; colsample_bylevel=0.809; lambda=0.0852; alpha=0.192; subsample=0.672; scale_pos_weight=5.51 : y = 0.981 : 0.7 secs : initdesign
#> [mbo] 0: eta=0.17; gamma=1.11; max_depth=8; colsample_bytree=0.569; colsample_bylevel=0.985; lambda=2.84; alpha=0.0632; subsample=0.785; scale_pos_weight=3.19 : y = 0.986 : 0.6 secs : initdesign
#> [mbo] 0: eta=0.158; gamma=0.0599; max_depth=19; colsample_bytree=0.671; colsample_bylevel=0.693; lambda=275; alpha=0.0142; subsample=0.525; scale_pos_weight=131 : y = 0.954 : 0.3 secs : initdesign
#> [mbo] 0: eta=0.0861; gamma=0.441; max_depth=5; colsample_bytree=0.818; colsample_bylevel=0.905; lambda=0.00104; alpha=0.00136; subsample=0.935; scale_pos_weight=0.0513 : y = 0.962 : 0.3 secs : initdesign
#> [mbo] 0: eta=0.0414; gamma=0.0134; max_depth=3; colsample_bytree=0.999; colsample_bylevel=0.65; lambda=0.0032; alpha=19.7; subsample=0.914; scale_pos_weight=614 : y = 0.969 : 0.5 secs : initdesign
#> [mbo] 0: eta=0.0284; gamma=0.0255; max_depth=11; colsample_bytree=0.943; colsample_bylevel=0.724; lambda=0.0355; alpha=1.83; subsample=0.547; scale_pos_weight=1.06 : y = 0.969 : 0.8 secs : initdesign
#> [mbo] 0: eta=0.0812; gamma=0.185; max_depth=14; colsample_bytree=0.878; colsample_bylevel=0.565; lambda=0.115; alpha=0.0313; subsample=0.828; scale_pos_weight=0.0112 : y = 0.985 : 1.1 secs : initdesign
#> [mbo] 0: eta=0.191; gamma=8.07; max_depth=7; colsample_bytree=0.851; colsample_bylevel=0.894; lambda=134; alpha=174; subsample=0.877; scale_pos_weight=19.5 : y = 0.945 : 0.3 secs : initdesign
#> [mbo] 0: eta=0.12; gamma=0.542; max_depth=11; colsample_bytree=0.615; colsample_bylevel=0.846; lambda=0.534; alpha=0.00287; subsample=0.764; scale_pos_weight=28 : y = 0.981 : 0.7 secs : initdesign
#> Loading required package: rgenoud
#> ##  rgenoud (Version 5.8-2.0, Build Date: 2018-04-03)
#> ##  See http://sekhon.berkeley.edu/rgenoud for additional documentation.
#> ##  Please cite software as:
#> ##   Walter Mebane, Jr. and Jasjeet S. Sekhon. 2011.
#> ##   ``Genetic Optimization Using Derivatives: The rgenoud package for R.''
#> ##   Journal of Statistical Software, 42(11): 1-26. 
#> ##
#> [mbo] 1: eta=0.125; gamma=0.0108; max_depth=5; colsample_bytree=0.978; colsample_bylevel=0.618; lambda=0.00409; alpha=0.137; subsample=0.799; scale_pos_weight=0.625 : y = 0.983 : 0.5 secs : infill_cb

チューニング後の結果を表示すると、おススメのパラメータセットを喋ってくれる。

print(res)
#> Autoxgboost tuning result
#> 
#> Recommended parameters:
#>               eta: 0.170
#>             gamma: 1.111
#>         max_depth: 8
#>  colsample_bytree: 0.569
#> colsample_bylevel: 0.985
#>            lambda: 2.836
#>             alpha: 0.063
#>         subsample: 0.785
#>  scale_pos_weight: 3.189
#>           nrounds: 25
#> 
#> 
#> Preprocessing pipeline:
#> dropconst(rel.tol = 1e-08, abs.tol = 1e-08, ignore.na = FALSE)
#> 
#> With tuning result: auc = 0.986

結果オブジェクトの構成は以下の通り。

str(res, 1)
#> List of 5
#>  $ optim.result    :List of 9
#>   ..- attr(*, "class")= chr [1:2] "MBOSingleObjResult" "MBOResult"
#>  $ final.learner   :List of 11
#>   ..- attr(*, "class")= chr [1:3] "CPOLearner" "BaseWrapper" "Learner"
#>  $ final.model     :List of 8
#>   ..- attr(*, "class")= chr [1:3] "CPOModel" "BaseWrapperModel" "WrappedModel"
#>  $ measure         :List of 10
#>   ..- attr(*, "class")= chr "Measure"
#>  $ preproc.pipeline:List of 24
#>   ..- attr(*, "class")= chr [1:2] "CPOPrimitive" "CPO"
#>  - attr(*, "class")= chr "AutoxgbResult"

チューニングしたモデルの評価

ベストのlearnermodelも取っておいてくれるので、これを使う。

res$final.learner
#> Learner classif.xgboost.custom.dropconst from package xgboost
#> Type: classif
#> Name: ; Short name: 
#> Class: CPOLearner
#> Properties: numerics,missings,twoclass,multiclass,prob,featimp
#> Predict-Type: prob
#> Hyperparameters: nrounds=25,verbose=0,objective=binary:logistic,eta=0.17,gamma=1.11,max_depth=8,colsample_bytree=0.569,colsample_bylevel=0.985,lambda=2.84,alpha=0.0632,subsample=0.785,scale_pos_weight=3.19
res$final.model
#> Model for learner.id=classif.xgboost.custom.dropconst; learner.class=CPOLearner
#> Trained on: task.id = spam[train.set, ]; obs = 3220; features = 57
#> Hyperparameters: nrounds=25,verbose=0,objective=binary:logistic,eta=0.17,gamma=1.11,max_depth=8,colsample_bytree=0.569,colsample_bylevel=0.985,lambda=2.84,alpha=0.0632,subsample=0.785,scale_pos_weight=3.19

取り分けておいたテストデータを使って予測性能を評価する。

test.task <- makeClassifTask(data = spam[test.set, ], target = "type")

pred <- predict(res$final.model, task = test.task)

performance(pred, measures = auc)
#>       auc 
#> 0.9851731

このデータではテストセットの性能も、チューニングのときのものとあまり変わらない。

(修正あり)arulesの結果をdata.frameで探索する

R標準のデータフレーム+dplyr等で、探索的にルール抽出⇔眺めるのにパッケージを作った。本家のarulesだとちょっとやりにくいなあと思っている人むけ。

使いかた

arules::inspect() の代わりに inspectDF() するだけ。

arules::DATAFRAME() でデータフレームとして出力できるとのこと。

自作する前によく探すの、大事。

require(dplyr)
require(arules)
require(inspectDF)

data(Groceries)
params <- list(confidence=0.001, support=0.001, maxlen=7, minlen=2)
glo.apriori <- apriori(Groceries, parameter = params, control = list(verbose=FALSE))
glo.inspectDF  <- inspectDF(glo.apriori)

glo.inspectDF %>% str
#> 'data.frame':    40943 obs. of  8 variables:
#>  $ rule      : chr  "Rule 1" "Rule 2" "Rule 3" "Rule 4" ...
#>  $ LHS       : chr  "whole milk" "other vegetables" "other vegetables" "rolls/buns" ...
#>  $ RHS       : chr  "sparkling wine" "artif. sweetener" "bathroom cleaner" "nuts/prunes" ...
#>  $ n.lhs     : num  1 1 1 1 1 1 1 1 1 1 ...
#>  $ support   : num  0.00102 0.00102 0.00102 0.00102 0.00102 ...
#>  $ confidence: num  0.00398 0.00525 0.00525 0.00553 0.00583 ...
#>  $ lift      : num  0.712 1.615 1.914 1.647 1.024 ...
#>  $ count     : num  10 10 10 10 10 10 10 10 10 10 ...

取り出したルールは、グラフ描画用のユーティリティ関数で眺める。

set.seed(0)
glo.inspectDF %>% 
  arrange(support, confidence) %>%
  head(60) %>% 
  plotRuleGraph()

f:id:kato-satoshi-0:20180724000155p:plain

インストール

GitHubから.

install.packages("devtools") # if you have not installed "devtools" package
devtools::install_github("katokohaku/inspectDF")

ソースコードも同じ場所。

モチベーション

本家のsubset(), sort()とかでなく、通常のdata.frameの操作としてやりたい。 本家の inspect() がコンソールに全部出力するのがツライ

print(glo.apriori)
#> set of 40943 rules

rules.lhs %>% summary()
#>      rule               LHS                RHS                  n.lhs        
#>  Length:8035        Length:8035        Length:8035        Min.   :2.000  
#>  Class :character   Class :character   Class :character   1st Qu.:2.000  
#>  Mode  :character   Mode  :character   Mode  :character   Median :3.000  
#>                                                           Mean   :2.739  
#>                                                           3rd Qu.:3.000  
#>                                                           Max.   :5.000  
#>     support           confidence           lift             count       
#>  Min.   :0.001017   Min.   :0.01815   Min.   : 0.4379   Min.   : 10.00  
#>  1st Qu.:0.001118   1st Qu.:0.14961   1st Qu.: 2.0673   1st Qu.: 11.00  
#>  Median :0.001322   Median :0.25175   Median : 2.6852   Median : 13.00  
#>  Mean   :0.001740   Mean   :0.30231   Mean   : 2.9595   Mean   : 17.11  
#>  3rd Qu.:0.001830   3rd Qu.:0.41935   3rd Qu.: 3.5531   3rd Qu.: 18.00  
#>  Max.   :0.022267   Max.   :1.00000   Max.   :22.1394   Max.   :219.00

arulesは非常によくできたパッケージですが、高度にwrapされ過ぎていて、気軽に弄ろうと思うとツライ。本家のprint()summary()も、探索的なルール分析にはちょっと不親切だと思いませんか。探索的なルール分析にはsubset(), sort()などが用意されているものの、表示はことごとく上記のそれです。ツライ。

また、本家の inspect()print()で返すため、コンソールにすべての出力が強制表示される。いつ終わるとも知れぬコンソール出力を眺めるの、ツライ。あと、標準では2番目の列名が無名なので、この後そのまま探索的に弄ろうとすると困る、というか矢印要る? ツライ。

glo.inspect  <- glo.apriori %>% head(10) %>% inspect()
#>      lhs               rhs            support     confidence  lift    
#> [1]  {honey}        => {whole milk}   0.001118454 0.733333333 2.870009
#> [2]  {whole milk}   => {honey}        0.001118454 0.004377238 2.870009
#> [3]  {soap}         => {whole milk}   0.001118454 0.423076923 1.655775
#> [4]  {whole milk}   => {soap}         0.001118454 0.004377238 1.655775
#> [5]  {tidbits}      => {soda}         0.001016777 0.434782609 2.493345
#> [6]  {soda}         => {tidbits}      0.001016777 0.005830904 2.493345
#> [7]  {tidbits}      => {rolls/buns}   0.001220132 0.521739130 2.836542
#> [8]  {rolls/buns}   => {tidbits}      0.001220132 0.006633499 2.836542
#> [9]  {cocoa drinks} => {whole milk}   0.001321810 0.590909091 2.312611
#> [10] {whole milk}   => {cocoa drinks} 0.001321810 0.005173100 2.312611
#>      count
#> [1]  11   
#> [2]  11   
#> [3]  11   
#> [4]  11   
#> [5]  10   
#> [6]  10   
#> [7]  12   
#> [8]  12   
#> [9]  13   
#> [10] 13

glo.inspect %>% str
#> 'data.frame':    10 obs. of  7 variables:
#>  $ lhs       : Factor w/ 7 levels "{cocoa drinks}",..: 2 7 4 7 6 5 6 3 1 7
#>  $           : Factor w/ 1 level "=>": 1 1 1 1 1 1 1 1 1 1
#>  $ rhs       : Factor w/ 7 levels "{cocoa drinks}",..: 7 2 7 4 5 6 3 6 7 1
#>  $ support   : num  0.00112 0.00112 0.00112 0.00112 0.00102 ...
#>  $ confidence: num  0.73333 0.00438 0.42308 0.00438 0.43478 ...
#>  $ lift      : num  2.87 2.87 1.66 1.66 2.49 ...
#>  $ count     : num  11 11 11 11 10 10 12 12 13 13

使い道

dplyrあたりと組み合わせて使うと、tidyに探索できるのではないか。

require(stringr)
rules.lhs  <- glo.inspectDF %>% 
  filter(str_detect(LHS, pattern = "yogurt|sausage")) %>%
  arrange(confidence, lift) %>%
  filter(n.lhs > 1) 

rules.lhs %>% head()
#>        rule               LHS                    RHS n     support
#> 1   Rule 97 whole milk,yogurt               UHT-milk 2 0.001016777
#> 2   Rule 98 whole milk,yogurt         red/blush wine 2 0.001016777
#> 3   Rule 99 whole milk,yogurt house keeping products 2 0.001016777
#> 4  Rule 100 whole milk,yogurt             liver loaf 2 0.001016777
#> 5 Rule 7175 whole milk,yogurt            chewing gum 2 0.001118454
#> 6 Rule 7176 whole milk,yogurt        cling film/bags 2 0.001118454
#>   confidence      lift count
#> 1 0.01814882 0.5425339    10
#> 2 0.01814882 0.9444108    10
#> 3 0.01814882 2.1767518    10
#> 4 0.01814882 3.5698730    10
#> 5 0.01996370 0.9485170    11
#> 6 0.01996370 1.7530626    11

プロットは hideaki さんのqiitaの記事をそのまま流用。エッジ上の円の面積はルールのsupport値、色の濃さはルールのconfidence値を表す。

adujust.support.sizeオプションでヒューリスティックなサイズの調整ができる。

require(stringr)
rules.lhs  <- glo.inspectDF %>% 
  filter(str_detect(RHS, pattern = "yogurt")) %>%
  arrange(confidence, lift) %>%
  filter(n.lhs > 1) %>% 
  head(15)

par(mfrow = c(1,2))
set.seed(0)
rules.lhs %>% plotRuleGraph(label = "default")
set.seed(0)
rules.lhs %>% plotRuleGraph(label = "adjusted rule size", 
                            adujust.support.size = 4000)

f:id:kato-satoshi-0:20180724000159p:plain

今後

{support, confidence, lift}{size, color}の対応は可変にして、サイズも自動調整したい。

既存のGithubのレポジトリをRStudioでcloneする

RStudioの操作だけで完結したい人向け。タイトル通りのことしたくて戸惑ったので備忘録。

(1) File > New Project... > Version Control > Git
f:id:kato-satoshi-0:20180126003129p:plain

(2) 既存のレポジトリのアドレスをコピー
f:id:kato-satoshi-0:20180126003414p:plain

(3) コピーしたアドレスを(1)に張り付け
f:id:kato-satoshi-0:20180126003449p:plain

(4) clone完了
f:id:kato-satoshi-0:20180126003609p:plain

あとは、ローカルの変更をpushしたり、別の端末からpushされた変更をpullしたり。(1)からVersion Controlを選べばいいというのに気が付くまで、時間がかかった。

元ネタ?はこちら:
qiita.com

random forestを使った解釈性にもとづく推薦システム

[1706.06691] Interpretable Predictions of Tree-based Ensembles via Actionable Feature Tweakingという論文が提案されている。KDD2017でYahoo Researchの研究者によって発表されたもの、とのこと。Ensemble treesによって、あるサンプルがラベルAと予測されたとき、別のラベルBと予測されたかったら元のサンプルのどの特徴量にどのような操作を加えればよいか「改善案」を算出する。手元にあると便利そうだったのでRで実装したので紹介したい。
どんなことができるかの概要は、先日のTokyoRで喋ってきたので、そちらのスライドも参照して欲しい。

www.slideshare.net
大雑把な手順は、下記のとおり:

  1. データのスケーリング
  2. 予測モデル(random forest)の構築
  3. ルール抽出
  4. 各事例への推薦の計算

なお、日本語の解説は、接点QBさんによる丁寧な記事がすでにあるので、そちらを参照されたい。
setten-qb.hatenablog.com

前準備

ライブラリとソースのロード。この記事で使用するコードはgithubに置いてあるので、適当にローカルにおいて読み込む。そのうち時間を作ってパッケージ化したい。
github.com

require(tidyverse)
require(magrittr)
require(randomForest)
devtools::install_github("hoxo-m/pforeach")
require(pforeach)

source("./R/tweak_feature.R")

pforeachに隅々まで依存しているので、別途インストールが必須。

訓練セットとテストセットの準備

例によってirisを使いたいところだが、今回はspamデータセットを使ってspamと分類されるサンプルをnonspamに分類されるように推薦を行う。irisは多クラスだしルールも最小限なんで良いんだけど、アヤメの花びらを切って貼ってもカキツバタにはならないので今回はスキップ。実用上、変更できない説明変数が予測モデルに含まれる場合には、禁止リストをかませるとか適切な層別化するといった工夫が別途必要になる。

# data preparation -------------------------------------------------------------
set.seed(777)

data(spam, package = "kernlab")
dataset <- sample_frac(spam)
n.test <- floor(NROW(dataset) *0.1)

dataset.train <- chop(dataset, n.test)
dataset.test  <- tail(dataset, n.test)

dim(dataset);dim(dataset.train);dim(dataset.test)

chop()はデータを後ろから指定したサイズだけ切り落とすための自作ユーティリティ。

random forestの学習とモデル選択

訓練セットで学習して重要度のプロットと弱学習器の数の評価を行う。この辺りはルーチンなので変わったことはしていないが、ある程度予測精度が担保されていることが前提なので、ここで予測精度がダメな場合は再検討。

# bilding randomForest -----------------------------------------
X <- dataset.train[, 1:(ncol(dataset.train)-1)]
true.y <- dataset.train[, ncol(dataset.train)]

forest.all <- randomForest(X, true.y, ntree=500)
forest.all
par(mfrow=c(1,2))
varImpPlot(forest.all) # to view varImp, x3 & x4 should be removed.
plot(forest.all)
par(mfrow=c(1,1))

f:id:kato-satoshi-0:20180122185241p:plain
この方法は総当たり探索的なアプローチなので、変数(ルール)の数とツリーの数が増えると計算量が爆発的に増える。したがって変数選択とツリーの絞り込みのほか、max_depthを減らす工夫が重要。
今回は上位12変数と100ツリーくらいあれば良さそうなので、適当に切り落とす。

# model shrinkage based on importance -------------------------
top.importance <- forest.all$importance %>% data.frame %>%
  tibble::rownames_to_column(var = "var") %>% 
  arrange(desc(MeanDecreaseGini)) %>% 
  head(12)

dataset.train.fs <- dataset.train %>% select(top.importance$var)
dataset.test.fs  <- dataset.test %>% select(top.importance$var)

データのスケーリングと縮約済みモデルの学習

この手法はハイパーパラメータを一つだけにするため、各変数を標準化する必要がある。実用上、テストセットは訓練セットのスケーリングパラメータで標準化される必要があるのと、推薦後の改善は標準化前の尺度に復元する必要があり、毎回attributeを取ってくるのは面倒なので、それぞれrescale()とdescale()として、ユーティリティ関数を用意した。

# scaling feature-selected data  ---------------------------
X.train <- scale( dataset.train.fs )
X.test  <- rescale( dataset.test.fs, scaled = X.train )

dataset.test.fs[1:6, 1:6]
descale(X.test, scaled = X.train)[1:6, 1:6]
descale(X.test, scaled = X.test)[1:6, 1:6]

ちゃんと復元されていることを確認して欲しい。スケーリングした訓練データでモデルの学習を行う。

forest.rf <- randomForest(X.train, true.y, ntree=100)

forest.all
forest.rf
plot(forest.rf)

f:id:kato-satoshi-0:20180122191903p:plain
変数を減らしたので予測精度が下がっているが、「予測ミス(擬陽性)を起こすサンプルは本質的には改善されて良いのでは」という立場をとって許容範囲と思うことにして目をつぶることにする。

feature tweaking

ここからが提案手法の処理。手順は、random forest予測モデルを学習したのち、(1)ルールを抽出する+変更ルールを算出する(2)個別のサンプルの予測と推薦案を算出する(3)標準化の復元、の1+3ステップとなる。

# (0)random forest予測モデルを学習する
forest.rf <- randomForest(X.train, true.y, ntree=100)

# (1)ルールを抽出する+変更ルールを算出する
es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3, resample = TRUE)

# (2)個別のサンプルの予測と推薦案を算出する
tweaked <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam", .dopar = TRUE)

# (3)標準化の復元
dt <- descale.tweakedFeature(tweaked, X.test)
(1)ルールを抽出する+変更ルールを算出する

random forestではgetTree()で弱学習器を取ってこれる。これをつかって各ラベルに分類されるルールセットは、ラベルごとにdata.frameのlistとして整理しなおして抽出する。

# test sampling (for demo)
ep <- getRules.randomForest(forest.rf, k=2, label.to = NULL)
ep %>% str(2)
ep[["spam"]]$path[[1]]

f:id:kato-satoshi-0:20180122194147p:plain
xgboostにも同様にルールセットを取り出すxgb.model.dt.tree()があるので、そのうち対応したい。

抽出されたルールセットと変更ルール(候補案)のセットは対応しているため、抽出と同時に計算してしまう。元論文のアルゴリズムは各サンプルの評価時にこれを毎回計算する構造になっており、素直に実装してしまうと計算量が大変になるので、今から実装しようと思う人は注意。
epsiron は変更の強さを指定するハイパーパラメタ。原理上、どうやっても予測を覆せないサンプルがあり、この量は変更可能な提案が可能なサンプル数に影響する。

es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3)

f:id:kato-satoshi-0:20180122200255p:plain
ntreeはルールセットを抽出する弱学習器の数で、省略するとすべてのツリーを使用する。デフォルトでは先頭からntreeで指定したツリーを解析対象にするが、データ不均衡に対処するために訓練セットをアンダーサンプリングして学習した幾つかのモデルをrandomForest::combine(...)で結合したい場合がある。この場合は、先頭から順番に解析してしまうと困るので、resampleオプションをTRUEに指定することでランダムに弱学習器を評価することができる。

es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3, resample = TRUE)

f:id:kato-satoshi-0:20180122200705p:plain

(2)個別のサンプルの予測と推薦案を算出する

予測対象のデータセットが label.from と予測されたとき、label.to と予測されるための改善案を計算する。学習データと同じスケールで標準されていないとマトモな推薦が出来ないので注意。

tweaked_ <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam",
                 .dopar = FALSE)

f:id:kato-satoshi-0:20180122202613p:plain
各サンプルについてルールを全部評価するので結構な計算時間がかかるため、並列化オプション.doparはデフォルトではTRUEになっている。

tweaked_ <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam",
                 .dopar = TRUE)
(3)標準化の復元

元のサンプルの操作方向とその量が知りたいのでスケーリングを解除する。元のサンプルの特徴量(original)、予測を変えるための変更後の特徴量(suggest)、予測の変更に必要な努力量=両者の差分(diff)のlistを返す。
>|r|
dt <- descale.tweakedFeature(tweaked, X.test)
dt %>% str(1)
|

f:id:kato-satoshi-0:20180122204026p:plain

個別のサンプルへの推薦パターンの可視化

個別の推薦結果を可視化する。操作量を相対的に比較したいため、標準化されたままのものを可視化している。サンプルによって違う推薦パターンが見えるので面白い。

# plot suggestion for each instance
which(tweaked$predict == "spam")
plot.suggest(tweaked, 4)
plot.suggest(tweaked, 11)

f:id:kato-satoshi-0:20180122203751p:plain
f:id:kato-satoshi-0:20180122203745p:plain
操作量の大きい順に、0でない変数だけ可視化したければ、それぞれ .ordered = TRUEと、.nonzero.only = TRUEを指定する。

plot.suggest(tweaked, 15)
plot.suggest(tweaked, 15, .ordered = TRUE, .nonzero.only = TRUE)

f:id:kato-satoshi-0:20180122203928p:plain
f:id:kato-satoshi-0:20180122203934p:plain

変数重要度の可視化

論文中ではユーザーによる納得感(アンケート調査?)をもとに評価しているが、それだと機械的には評価できないので、「予測サンプル集団において操作回数が多かった変数は重要なんじゃないか」という考えと、「操作された変数の操作の方向性が揃っている変数は、その意味合いを解釈できそう」と考えて、これらを各変数の重要度と思うことにして、それぞれ可視化する。

# Plot population importances
pp <- plot.tweakedPopulation(tweaked, type = "absoluteSum")
pp <- plot.tweakedPopulation(tweaked, type = "direction")
pp <- plot.tweakedPopulation(tweaked, type = "frequency")

type = "frequency"は、サンプル集団において各変数が操作をされた集計値を可視化する。
f:id:kato-satoshi-0:20180122205329p:plain
type = "direction"は、各サンプルの推薦で各変数がどのような操作をされたかを集計して可視化する。相対的な重要性を比較したいので、標準化された状態の推薦量を使用する。
f:id:kato-satoshi-0:20180122210838p:plain
変更量の総和を可視化する"absoluteSum"オプションも作ってみたが、これはあまり有意義と感じなかった。

おわりに

カテゴリ変数を含むデータセットの取り扱い、変更できない変数に対する操作の提案、変更の実現コスト(お金とか労力とかの実リソース)に関するドメイン知識で重みづけ補正する、など、いろいろ工夫の余地があって実際の分析の際にはアイデアの見せどころと思う。
また、この種の感度分析的なアプローチ全般に言えることだが、因果関係の示唆となる関連性を補足しているかもしれないが、かならずしも因果関係そのものを見ているとは限らない点も気を付けたい。