琥珀色呑んだくれ備忘録

メモとか備忘録とか

XGBoostExplainerが何をやっているか調べる(4.モデルから予測ルールを抽出する)

目的

今回は、xgboostExplainerによって、xgboostの学習済みモデルからルールがどうやって抽出されているかにフォーカスし、適宜xgboostの資料を見ながら追いかける。

(参考資料)

関連シリーズ

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

準備:XGBモデルの学習と予測

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

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)

学習したxgboostのルール抽出

buildExplainer()の中身を抜き書きしながら、step-by-stepで眺める

# explainer = buildExplainer(xgb.model,xgb.train.data, type="binary", base_score = 0.5, trees = NULL)
# function (xgb.model, trainingData, type = "binary", base_score = 0.5, trees_idx = NULL) 
# {
trainingData = xgb.train.data
type = "binary" 
base_score = 0.5
trees_idx = NULL

xgb.model.dt.tree()によるパスの抽出

xgboost::xgb.model.dt.tree()を使う。

col_names = attr(trainingData, ".Dimnames")[[2]]
col_names %>% head()
#> [1] "cap-shape=bell"    "cap-shape=conical" "cap-shape=convex" 
#> [4] "cap-shape=flat"    "cap-shape=knobbed" "cap-shape=sunken"

cat("\nCreating the trees of the xgboost model...")
#> 
#> Creating the trees of the xgboost model...
trees = xgb.model.dt.tree(col_names, model = xgb.model, trees = trees_idx)

trees %>% 
  mutate(Feature = str_trunc(Feature, width=12, side="left")) %>% 
  select(-Missing)
#>    Tree Node   ID      Feature Split  Yes   No      Quality       Cover
#> 1     0    0  0-0    odor=foul   0.5  0-1  0-2 2711.3557100 1250.000000
#> 2     0    1  0-1 ...ize=broad   0.5  0-3  0-4 1263.3979500  901.000000
#> 3     0    2  0-2         Leaf    NA <NA> <NA>    0.5982857  349.000000
#> 4     0    3  0-3    odor=none   0.5  0-5  0-6  264.4693910  202.250000
#> 5     0    4  0-4 ...lor=green   0.5  0-7  0-8  203.9635310  698.750000
#> 6     0    5  0-5 ...or=yellow   0.5  0-9 0-10  121.6305770  157.500000
#> 7     0    6  0-6 ...ing=silky   0.5 0-11 0-12   74.7532349   44.750000
#> 8     0    7  0-7         Leaf    NA <NA> <NA>   -0.5991260  685.500000
#> 9     0    8  0-8         Leaf    NA <NA> <NA>    0.5578948   13.250000
#> 10    0    9  0-9  odor=almond   0.5 0-13 0-14   65.5069656  148.000000
#> 11    0   10 0-10         Leaf    NA <NA> <NA>   -0.5428572    9.500000
#> 12    0   11 0-11 ...?=bruises   0.5 0-15 0-16   26.8929482   38.250000
#> 13    0   12 0-12         Leaf    NA <NA> <NA>    0.5200000    6.500000
#> 14    0   13 0-13   odor=anise   0.5 0-17 0-18   62.2878838  143.250000
#> 15    0   14 0-14         Leaf    NA <NA> <NA>   -0.4956522    4.750000
#> 16    0   15 0-15         Leaf    NA <NA> <NA>   -0.5838926   36.250000
#> 17    0   16 0-16         Leaf    NA <NA> <NA>    0.4000000    2.000000
#> 18    0   17 0-17         Leaf    NA <NA> <NA>    0.5957143  139.000000
#> 19    0   18 0-18         Leaf    NA <NA> <NA>   -0.4857143    4.250000
#> 20    1    0  1-0    odor=foul   0.5  1-1  1-2 1489.0952100 1145.302000
#> 21    1    1  1-1 ...ize=broad   0.5  1-3  1-4  691.4842530  825.759888
#> 22    1    2  1-2         Leaf    NA <NA> <NA>    0.4634756  319.542114
#> 23    1    3  1-3    odor=none   0.5  1-5  1-6  142.5185550  186.003708
#> 24    1    4  1-4 ...lor=green   0.5  1-7  1-8  114.8126300  639.756165
#> 25    1    5  1-5  odor=almond   0.5  1-9 1-10   69.8917542  144.674194
#> 26    1    6  1-6 ...ing=silky   0.5 1-11 1-12   42.7613907   41.329517
#> 27    1    7  1-7         Leaf    NA <NA> <NA>   -0.4640480  627.485962
#> 28    1    8  1-8         Leaf    NA <NA> <NA>    0.4361763   12.270212
#> 29    1    9  1-9   odor=anise   0.5 1-13 1-14   75.4965286  135.787827
#> 30    1   10 1-10         Leaf    NA <NA> <NA>   -0.4301576    8.886355
#> 31    1   11 1-11 ...?=bruises   0.5 1-15 1-16   16.6019783   35.249847
#> 32    1   12 1-12         Leaf    NA <NA> <NA>    0.4107886    6.079669
#> 33    1   13 1-13         Leaf    NA <NA> <NA>    0.4617254  127.362411
#> 34    1   14 1-14         Leaf    NA <NA> <NA>   -0.4283619    8.425421
#> 35    1   15 1-15         Leaf    NA <NA> <NA>   -0.4537036   33.327763
#> 36    1   16 1-16         Leaf    NA <NA> <NA>    0.3296103    1.922086
#> 37    2    0  2-0    odor=foul   0.5  2-1  2-2  934.9316410  956.514771
#> 38    2    1  2-1 ...ize=broad   0.5  2-3  2-4  434.1681820  689.965515
#> 39    2    2  2-2         Leaf    NA <NA> <NA>    0.4022448  266.549286
#> 40    2    3  2-3    odor=none   0.5  2-5  2-6   87.7709885  156.324432
#> 41    2    4  2-4 ...lor=green   0.5  2-7  2-8   73.6840591  533.641052
#> 42    2    5  2-5  odor=almond   0.5  2-9 2-10   44.8767204  121.285667
#> 43    2    6  2-6 ...ing=silky   0.5 2-11 2-12   27.4686470   35.038773
#> 44    2    7  2-7         Leaf    NA <NA> <NA>   -0.4028375  523.192078
#> 45    2    8  2-8         Leaf    NA <NA> <NA>    0.3751199   10.448951
#> 46    2    9  2-9   odor=anise   0.5 2-13 2-14   48.7504120  113.642014
#> 47    2   10 2-10         Leaf    NA <NA> <NA>   -0.3680061    7.643652
#> 48    2   11 2-11 ...?=bruises   0.5 2-15 2-16   11.4850969   29.765743
#> 49    2   12 2-12         Leaf    NA <NA> <NA>    0.3515948    5.273030
#> 50    2   13 2-13         Leaf    NA <NA> <NA>    0.4004391  106.384308
#> 51    2   14 2-14         Leaf    NA <NA> <NA>   -0.3663105    7.257702
#> 52    2   15 2-15         Leaf    NA <NA> <NA>   -0.3922864   28.009958
#> 53    2   16 2-16         Leaf    NA <NA> <NA>    0.2832851    1.755784
# %>%
#   mutate_at(.vars = vars("Quality","Cover"), .funs = round)

LeafノードのQualityは0じゃないの?」と思うかもしれないが、マニュアルに

Quality: either the split gain (change in loss) or the leaf value

と書いてあり、LeafノードではQualityのセルが予測結果の格納場所として流用されていることに注意。(あとからこれを利用する)*1

cat("\nGetting the leaf nodes for the training set observations...")
#> 
#> Getting the leaf nodes for the training set observations...

なお、関数内で呼ばれているpredict(..., predleaf = TRUE)は、predleaf = TRUEを指定することで、予測値の代わりに訓練データのインスタンスが各treeで所属するLeafのノード番号を取得できる。 今回はNROW(trainingData)=5000, nrounds = 3なので、5000行3列の所属Leafの行列が得られる。

nodes.train %>% dim()
#> [1] 5000    3
nodes.train <- NULL

ただし、取得されたのち削除しても特に問題はなく、学習に使われたインスタンスの予測結果の情報が、buildExplainer()のどこかで使われている形跡は見当たらなかった*2

予測値の再分配

xgboostExplainerは、Leafの予測値を親ノードに再分配することで、予測結果を各ノードにおけるImpactとして分解し、そこを通るインスタンスの予測結果を説明する。

以下では、buildExplainer()のエンジン部分であるxgboostExplainer:::getStatsForTrees()の一連のステップをトレースする。

# tree_list = xgboostExplainer:::getStatsForTrees(trees, nodes.train, type = type, base_score = base_score)
# function (trees, nodes.train, type = "binary", base_score = 0.5) 
# {

なお、関数内部でdata.table::copy()'により、- attr(*, ".internal.selfref")=`で、データの更新が参照渡しとして行われているのに注意*3

tree_datatable = data.table::copy(trees)

Cover (H)の再計算

xgb.model.dt.tree()の情報から、各ノードまでの二次の勾配の和を取り出す。これはCover列の情報そのもののはずなのだが、xgb.model.dt.tree()の出力するCoverは精度に問題がある*4とのことで、きちんととした精度で計算しなおす。

f:id:kato-satoshi-0:20181231135410p:plain:w480

LeafノードのCoverから、逆向きにたどりながら足し合わせていけばよい

type = "binary"
base_score = 0.5

tree_datatable[, leaf := (Feature == "Leaf")]
non.leaves = which(tree_datatable[, leaf] == F)

tree_datatable[, H := Cover]

cat("\n\nRecalculating the cover for each non-leaf... \n")
#> 
#> 
#> Recalculating the cover for each non-leaf...
print.counter = 1
for (i in rev(non.leaves)) {
  left = tree_datatable[i, Yes]
  right = tree_datatable[i, No]
  
  tree_datatable[i, H := (tree_datatable[ID == left, H] + 
                            tree_datatable[ID == right, H])]
  
  if(print.counter < 5){
    print(i - 1)
    
    bind_rows(tree_datatable[i, ],
              tree_datatable[ID == left, ], 
              tree_datatable[ID == right,]) %>% 
      select(-Tree,-Node,-Feature,-Split,-Missing) %>% 
      print()
    print.counter <- print.counter + 1
  }
}
#> [1] 47
#>      ID  Yes   No    Quality     Cover  leaf         H
#> 1: 2-11 2-15 2-16 11.4850969 29.765743 FALSE 29.765742
#> 2: 2-15 <NA> <NA> -0.3922864 28.009958  TRUE 28.009958
#> 3: 2-16 <NA> <NA>  0.2832851  1.755784  TRUE  1.755784
#> [1] 45
#>      ID  Yes   No    Quality      Cover  leaf          H
#> 1:  2-9 2-13 2-14 48.7504120 113.642014 FALSE 113.642010
#> 2: 2-13 <NA> <NA>  0.4004391 106.384308  TRUE 106.384308
#> 3: 2-14 <NA> <NA> -0.3663105   7.257702  TRUE   7.257702
#> [1] 42
#>      ID  Yes   No    Quality    Cover  leaf        H
#> 1:  2-6 2-11 2-12 27.4686470 35.03877 FALSE 35.03877
#> 2: 2-11 2-15 2-16 11.4850969 29.76574 FALSE 29.76574
#> 3: 2-12 <NA> <NA>  0.3515948  5.27303  TRUE  5.27303
#> [1] 41
#>      ID  Yes   No    Quality      Cover  leaf          H
#> 1:  2-5  2-9 2-10 44.8767204 121.285667 FALSE 121.285663
#> 2:  2-9 2-13 2-14 48.7504120 113.642014 FALSE 113.642010
#> 3: 2-10 <NA> <NA> -0.3680061   7.643652  TRUE   7.643652

CoverとHを比べてみると、たしかに足し算の結果がずれている..。

勾配(G)とweightの再分配

xgboost::xgb.model.dt.tree()の出力のうち、LeafノードではQualityのセルが予測結果の格納場所として流用されているので、これを起点にして、各親ノードにweightを再分配する。

定義上、 w^*_j=-\frac{G_j}{H_j+\lambda} なのだが、これ以降 \lambdaが入っていないので、G_j = \sum_{i \in I_j}^{} g_i とは厳密な意味では違うのかもしれないが、ここは読み取れなかった。

base_weight = log(base_score/(1 - base_score))

tree_datatable[, previous_weight := base_weight]
tree_datatable[1, previous_weight:=0]

tree_datatable[leaf==T, weight := base_weight + Quality]

tree_datatable[leaf==T, G := -weight * H]

以降は、round(Tree)単位に分けて処理を進める。

tree_list = split(tree_datatable,as.factor(tree_datatable$Tree))
tree_list %>% str(1)
#> List of 3
#>  $ 0:Classes 'data.table' and 'data.frame':  19 obs. of  15 variables:
#>   ..- attr(*, ".internal.selfref")=<externalptr> 
#>  $ 1:Classes 'data.table' and 'data.frame':  17 obs. of  15 variables:
#>   ..- attr(*, ".internal.selfref")=<externalptr> 
#>  $ 2:Classes 'data.table' and 'data.frame':  17 obs. of  15 variables:
#>   ..- attr(*, ".internal.selfref")=<externalptr>

各roundのLeaf以外の親ノードjの G_jを計算する。定義から、 I_j = I_{left, j}+I_{right, j} なので、 G_j = \sum_{i \in I_{left,j}}^{} g_i + \sum_{i \in I_{right,j}}^{} g_i = G_{L,j} + G_{R,j} の足し算を計算すればよい*5

f:id:kato-satoshi-0:20181231140712p:plain:w480

num_tree_list = length(tree_list)
treenums =  as.character(0:(num_tree_list-1))
t = 0
cat('\n\nFinding the stats for the xgboost trees...\n')
#> 
#> 
#> Finding the stats for the xgboost trees...
# pb <- txtProgressBar(style=3)
for (tree in tree_list){
  t=t+1
  num_nodes = nrow(tree)
  non_leaf_rows = rev(which(tree[,leaf]==F))
  for (r in non_leaf_rows){
    left = tree[r,Yes]
    right = tree[r,No]
    leftG = tree[ID==left,G]
    rightG = tree[ID==right,G]
    
    tree[r,G:=leftG+rightG]
    w=tree[r,-G/H]
    
    tree[r,weight:=w]
    tree[ID==left,previous_weight:=w]
    tree[ID==right,previous_weight:=w]
    
    # bind_rows(tree[r,], tree[ID==left,], tree[ID==right,]) %>% print()
  }
  
  tree[,uplift_weight:=weight-previous_weight]
  # setTxtProgressBar(pb, t / num_tree_list)
}

tree_list %>% str(1)
#> List of 3
#>  $ 0:Classes 'data.table' and 'data.frame':  19 obs. of  16 variables:
#>   ..- attr(*, ".internal.selfref")=<externalptr> 
#>   ..- attr(*, "index")= int(0) 
#>   .. ..- attr(*, "__ID")= int [1:19] 1 2 11 12 13 14 15 16 17 18 ...
#>  $ 1:Classes 'data.table' and 'data.frame':  17 obs. of  16 variables:
#>   ..- attr(*, ".internal.selfref")=<externalptr> 
#>   ..- attr(*, "index")= int(0) 
#>   .. ..- attr(*, "__ID")= int [1:17] 1 2 11 12 13 14 15 16 17 3 ...
#>  $ 2:Classes 'data.table' and 'data.frame':  17 obs. of  16 variables:
#>   ..- attr(*, ".internal.selfref")=<externalptr> 
#>   ..- attr(*, "index")= int(0) 
#>   .. ..- attr(*, "__ID")= int [1:17] 1 2 11 12 13 14 15 16 17 3 ...

各ルールのインパクトの集計(Tree Breakdown)

buildExplainerFromTreeList()を実行すると、

  • 各treeでgetTreeBreakdown()が呼ばれ、
    • LeafgetLeafBreakdown()が呼ばれ、
      • findPath()Leaf → root node までのパスを取得
      • Leaf → root nodeまでのupliftに基づいて、各パスのImpactを集計
      • root nodeのupliftをInterceptとする

f:id:kato-satoshi-0:20181231140052p:plain:w480

upliftの計算自体は、weightの再分配の際にいっしょに処理されている。

#> STEP 2 of 2
explainer = xgboostExplainer:::buildExplainerFromTreeList(tree_list, col_names)
#> 
#> 
#> Getting breakdown for each leaf of each tree...
#> 
  |=================================================================| 100%

explainer %>% str(0)
#> Classes 'data.table' and 'data.frame':   28 obs. of  129 variables:
#>  - attr(*, ".internal.selfref")=<externalptr>

本質的な部分ではないので、step-by-stepは省略。 取り出したImpactを該当列(とIntercept)に格納する。最終的には [ルール数+1]列 ×Leaf数の行列がexplainerとして返る。 explainerの中身は以前見た通り。

kato-kohaku-0.hatenablog.com

可視化もこれまで通り。

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

*1:論文を読んでからコードを読み進めていて、ここではまった。

*2:この実装、全般的にレガシーが残っている傾向があって、追いかけていったら途中で行き止まりだった、みたいなことが、あちこちである

*3:メモリの節約のためと思われるが、今回はトレースのため、元コードから変数名を変更した

*4:途中の丸め方に問題があるのかな

*5:論文中には、GR ← G - GLと書いてある(Algrithm.1)

XGBoostExplainerが何をやっているか調べる(3.予測結果の再構成プロセスを眺める)

目的

今回は、インスタンスの予測結果が再構成されるプロセスを、explainPredictions()の中身を抜き書きしながら、何から何が取り出されているか見ていく。

関連シリーズ

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

準備

XGBモデルの学習と予測

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

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)

予測ルールを抽出する

buildExplainer()で、学習したxgboostのモデルから予測ルールを抽出する。

library(xgboostExplainer)

explainer = buildExplainer(xgb.model,xgb.train.data, type="binary", base_score = 0.5, trees = NULL)

予測値を再構成する

overview

指定したインスタンスの予測結果を、buildExplainer()が学習したモデルから再構成した結果(breakdown)を眺める。

再構成の大まかな手順は以下の通り:

  1. 指定したインスタンスが各roundでどのleafに落ちるか予測結果を得る
  2. 各roundで該当するleafの各特徴量の予測値への貢献量(予測ルール)を得る
  3. 2を集計し、round全体での予測ルールを得る
  4. 3を集計し、予測値を得る

対象の予測

今回は省略のため、1インスタンスだけを対象として指定するが、下記idxは複数のインスタンスを指定するベクトルでもよい。

require(data.table)
# breakdown = explainPredictions(xgb.model, explainer, slice(DMatrix, as.integer(idx)))
# function (xgb.model, explainer, data) 

DMatrix = xgb.test.data
idx = 2

data = slice(DMatrix, as.integer(idx))

nodes = predict(xgb.model, data, predleaf = TRUE)
print(nodes)
#>      [,1] [,2] [,3]
#> [1,]   17   13   13

xgboost:::predict.xgb.Booster()は、predleaf = TRUEオプションを指定することで、あるインスタンスの予測時に、それぞれのroundのtreeでどのleafに落ちたか?を返してくれる

初期化

インスタンスの数(行)× Interceptを含むすべてのルール(列)からなるゼロ行列を準備する。

colnames = names(explainer)[1:(ncol(explainer) - 2)]
  
preds_breakdown = data.table(matrix(0, nrow = nrow(nodes), ncol = length(colnames)))
setnames(preds_breakdown, colnames)
  
num_trees = ncol(nodes)

cat("\n\nExtracting the breakdown of each prediction...\n")
#> 
#> Extracting the breakdown of each prediction...

preds_breakdown.init <- preds_breakdown

対象のtreeの取り出し

今後のトレースを楽にするために「すべての行が0の列を削除し、残った列の列名を短縮する」関数を自作した。

selectNonzeroToShort <- function(data, w=12){
  select_if(data,
            .predicate = function(x){ sum(abs(x)) > 0} , 
            .funs = function(x){ str_trunc(x, width=w, side="left") }) 
}

最初のtreeだけ少し丁寧にトレースする(残りは繰り返しなので省略)

x=1

nodes_for_tree = nodes[, x]
nodes_for_tree
#> [1] 17

buildExplainer()でとりだしたLeafのうち、このroundで対象となる木のLeafを列挙。

tree_breakdown = explainer[tree == x - 1]
tree_breakdown %>% selectNonzeroToShort()
#>     ...or=yellow ...?=bruises odor=almond  odor=anise  odor=foul
#>  1:   0.00000000   0.00000000  0.00000000  0.00000000  0.7088975
#>  2:   0.00000000   0.00000000  0.00000000  0.00000000 -0.2745896
#>  3:   0.00000000   0.00000000  0.00000000  0.00000000 -0.2745896
#>  4:  -1.00780012   0.00000000  0.00000000  0.00000000 -0.2745896
#>  5:   0.00000000   0.00000000  0.00000000  0.00000000 -0.2745896
#>  6:   0.06468987   0.00000000 -1.02528500  0.00000000 -0.2745896
#>  7:   0.00000000  -0.05144537  0.00000000  0.00000000 -0.2745896
#>  8:   0.00000000   0.93244731  0.00000000  0.00000000 -0.2745896
#>  9:   0.06468987   0.00000000  0.03399723  0.03208427 -0.2745896
#> 10:   0.06468987   0.00000000  0.03399723 -1.04934438 -0.2745896
#>      odor=none ...ize=broad ...ing=silky ...lor=green  intercept leaf
#>  1:  0.0000000    0.0000000    0.0000000   0.00000000 -0.1106117    2
#>  2:  0.0000000   -0.1919848    0.0000000  -0.02193993 -0.1106117    7
#>  3:  0.0000000   -0.1919848    0.0000000   1.13508088 -0.1106117    8
#>  4:  0.1868594    0.6632849    0.0000000   0.00000000 -0.1106117   10
#>  5: -0.6576614    0.6632849    0.8995779   0.00000000 -0.1106117   12
#>  6:  0.1868594    0.6632849    0.0000000   0.00000000 -0.1106117   14
#>  7: -0.6576614    0.6632849   -0.1528694   0.00000000 -0.1106117   15
#>  8: -0.6576614    0.6632849   -0.1528694   0.00000000 -0.1106117   16
#>  9:  0.1868594    0.6632849    0.0000000   0.00000000 -0.1106117   17
#> 10:  0.1868594    0.6632849    0.0000000   0.00000000 -0.1106117   18

対象のleafのとりだし

round=1のtreeのなかのLeaf ==17が、今回の予測ルール(予測値に各特徴量が寄与する変化量)。これを取り出して、preds_breakdownに足し合わせていく。

preds_breakdown_for_tree = 
  tree_breakdown[match(nodes_for_tree, tree_breakdown$leaf), ]

print("preds_breakdown_for_tree")
#> [1] "preds_breakdown_for_tree"
preds_breakdown_for_tree %>%  
  selectNonzeroToShort(w=20) %>% 
  glimpse()
#> Observations: 1
#> Variables: 8
#> $ `cap-color=yellow` <dbl> 0.06468987
#> $ `odor=almond`      <dbl> 0.03399723
#> $ `odor=anise`       <dbl> 0.03208427
#> $ `odor=foul`        <dbl> -0.2745896
#> $ `odor=none`        <dbl> 0.1868594
#> $ `gill-size=broad`  <dbl> 0.6632849
#> $ intercept          <dbl> -0.1106117
#> $ leaf               <int> 17

preds_breakdown = 
  preds_breakdown + 
  preds_breakdown_for_tree[,colnames, with = FALSE]

print("preds_breakdown")
#> [1] "preds_breakdown"
preds_breakdown %>%  
  selectNonzeroToShort(w=20) %>% 
  glimpse()
#> Observations: 1
#> Variables: 7
#> $ `cap-color=yellow` <dbl> 0.06468987
#> $ `odor=almond`      <dbl> 0.03399723
#> $ `odor=anise`       <dbl> 0.03208427
#> $ `odor=foul`        <dbl> -0.2745896
#> $ `odor=none`        <dbl> 0.1868594
#> $ `gill-size=broad`  <dbl> 0.6632849
#> $ intercept          <dbl> -0.1106117

すべてのroudsをまとめて実行すると、以下の通り。

preds_breakdown <- preds_breakdown.init
for (x in 1:num_trees) {
  print(x)
  
  nodes_for_tree = nodes[, x]
  tree_breakdown = explainer[tree == x - 1]

  preds_breakdown_for_tree = 
    tree_breakdown[match(nodes_for_tree, tree_breakdown$leaf), ]
  
  print("preds_breakdown_for_tree")
  preds_breakdown_for_tree %>%  
    selectNonzeroToShort(w=20) %>% 
    glimpse()
  
  preds_breakdown = 
    preds_breakdown + 
    preds_breakdown_for_tree[,colnames, with = FALSE]
  
  print("preds_breakdown")
  preds_breakdown %>%  
    selectNonzeroToShort(w=20) %>% 
    glimpse()
}
#> [1] 1
#> [1] "preds_breakdown_for_tree"
#> Observations: 1
#> Variables: 8
#> $ `cap-color=yellow` <dbl> 0.06468987
#> $ `odor=almond`      <dbl> 0.03399723
#> $ `odor=anise`       <dbl> 0.03208427
#> $ `odor=foul`        <dbl> -0.2745896
#> $ `odor=none`        <dbl> 0.1868594
#> $ `gill-size=broad`  <dbl> 0.6632849
#> $ intercept          <dbl> -0.1106117
#> $ leaf               <int> 17
#> [1] "preds_breakdown"
#> Observations: 1
#> Variables: 7
#> $ `cap-color=yellow` <dbl> 0.06468987
#> $ `odor=almond`      <dbl> 0.03399723
#> $ `odor=anise`       <dbl> 0.03208427
#> $ `odor=foul`        <dbl> -0.2745896
#> $ `odor=none`        <dbl> 0.1868594
#> $ `gill-size=broad`  <dbl> 0.6632849
#> $ intercept          <dbl> -0.1106117
#> [1] 2
#> [1] "preds_breakdown_for_tree"
#> Observations: 1
#> Variables: 8
#> $ `odor=almond`     <dbl> 0.05139002
#> $ `odor=anise`      <dbl> 0.05522851
#> $ `odor=foul`       <dbl> -0.2125787
#> $ `odor=none`       <dbl> 0.1433645
#> $ `gill-size=broad` <dbl> 0.5101908
#> $ intercept         <dbl> -0.08586974
#> $ leaf              <int> 13
#> $ tree              <dbl> 1
#> [1] "preds_breakdown"
#> Observations: 1
#> Variables: 7
#> $ `cap-color=yellow` <dbl> 0.06468987
#> $ `odor=almond`      <dbl> 0.08538725
#> $ `odor=anise`       <dbl> 0.08731278
#> $ `odor=foul`        <dbl> -0.4871683
#> $ `odor=none`        <dbl> 0.3302238
#> $ `gill-size=broad`  <dbl> 1.173476
#> $ intercept          <dbl> -0.1964815
#> [1] 3
#> [1] "preds_breakdown_for_tree"
#> Observations: 1
#> Variables: 8
#> $ `odor=almond`     <dbl> 0.04534281
#> $ `odor=anise`      <dbl> 0.04896816
#> $ `odor=foul`       <dbl> -0.1841252
#> $ `odor=none`       <dbl> 0.1238637
#> $ `gill-size=broad` <dbl> 0.4407547
#> $ intercept         <dbl> -0.0743651
#> $ leaf              <int> 13
#> $ tree              <dbl> 2
#> [1] "preds_breakdown"
#> Observations: 1
#> Variables: 7
#> $ `cap-color=yellow` <dbl> 0.06468987
#> $ `odor=almond`      <dbl> 0.1307301
#> $ `odor=anise`       <dbl> 0.1362809
#> $ `odor=foul`        <dbl> -0.6712935
#> $ `odor=none`        <dbl> 0.4540875
#> $ `gill-size=broad`  <dbl> 1.61423
#> $ intercept          <dbl> -0.2708466

# return(preds_breakdown)

検算

元実装の結果と合致してるか検算。分解再構成した結果の可視化。

preds_breakdown %>% selectNonzeroToShort(w=20)
#>    cap-color=yellow odor=almond odor=anise  odor=foul odor=none
#> 1:       0.06468987   0.1307301  0.1362809 -0.6712935 0.4540875
#>    gill-size=broad  intercept
#> 1:         1.61423 -0.2708466

weight <- sum(preds_breakdown)
weight
#> [1] 1.457879

prediction <- exp(weight)/(1+exp(weight))
prediction
#> [1] 0.811208

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

次回は、buildExplainer()が、どのような手続きで学習したxgboostのモデルから予測ルールを抽出しているのか詳細に見ていく。

XGBoostExplainerが何をやっているか調べる(2.可視化プロセスを眺める)

目的

今回は、xgboostExplainerによって、xgboostのモデルと予測結果から何が取り出され、どう捌かれているかにフォーカスして追いかける。

関連シリーズ

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

予測結果の可視化プロセスをstep-by-stepで眺める

showWaterfall(..., type = "binary")の中身を抜き書きしながら、都度、何が取り出されているか見ていく。*1

準備

XGBモデルの学習と予測

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

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)

可視化の手続き

overview

あるインスタンスの予測結果を可視化する手順は:

  1. 学習したxgboostのモデルから予測ルールを抽出する
  2. showWaterfall()のなかで
    1. 指定したインスタンスの予測結果を予測ルールから分解再構成する
    2. 分解再構成した予測結果をwatarfall chartで可視化する

なお、

  • buildExplainer()どうやってモデルから予測ルールを抽出しているか
  • explainPredictions()どうやって分解再構成するか

については今回は触れない。

学習したxgboostのモデルから予測ルールを抽出する

buildExplainer()が抽出した予測ルール(ノード)を眺める。

library(xgboostExplainer)

explainer = buildExplainer(xgb.model,xgb.train.data, type="binary", base_score = 0.5, trees = NULL)

今回はルールの数が多いので、一度も登場しないものを除外してみる:

nonzero.cols <- abs(explainer) %>% 
  colSums() %>% 
  is_greater_than(0.0) %>% 
  which()

explainer %>% 
  select(nonzero.cols %>% rev()) %>% 
  set_colnames(
    colnames(.) %>% 
      str_trunc(width=10, side="left"))
#>     tree leaf   intercept  ...r=green ...g=silky ...e=broad  odor=none
#>  1:    0    2 -0.11061173  0.00000000  0.0000000  0.0000000  0.0000000
#>  2:    0    7 -0.11061173 -0.02193993  0.0000000 -0.1919848  0.0000000
#>  3:    0    8 -0.11061173  1.13508088  0.0000000 -0.1919848  0.0000000
#>  4:    0   10 -0.11061173  0.00000000  0.0000000  0.6632849  0.1868594
#>  5:    0   12 -0.11061173  0.00000000  0.8995779  0.6632849 -0.6576614
#>  6:    0   14 -0.11061173  0.00000000  0.0000000  0.6632849  0.1868594
#>  7:    0   15 -0.11061173  0.00000000 -0.1528694  0.6632849 -0.6576614
#>  8:    0   16 -0.11061173  0.00000000 -0.1528694  0.6632849 -0.6576614
#>  9:    0   17 -0.11061173  0.00000000  0.0000000  0.6632849  0.1868594
#> 10:    0   18 -0.11061173  0.00000000  0.0000000  0.6632849  0.1868594
#> 11:    1    2 -0.08586974  0.00000000  0.0000000  0.0000000  0.0000000
#> 12:    1    7 -0.08586974 -0.01726586  0.0000000 -0.1483337  0.0000000
#> 13:    1    8 -0.08586974  0.88295840  0.0000000 -0.1483337  0.0000000
#> 14:    1   10 -0.08586974  0.00000000  0.0000000  0.5101908  0.1433645
#> 15:    1   12 -0.08586974  0.00000000  0.7008943  0.5101908 -0.5018480
#> 16:    1   13 -0.08586974  0.00000000  0.0000000  0.5101908  0.1433645
#> 17:    1   14 -0.08586974  0.00000000  0.0000000  0.5101908  0.1433645
#> 18:    1   15 -0.08586974  0.00000000 -0.1208858  0.5101908 -0.5018480
#> 19:    1   16 -0.08586974  0.00000000 -0.1208858  0.5101908 -0.5018480
#> 20:    2    2 -0.07436510  0.00000000  0.0000000  0.0000000  0.0000000
#> 21:    2    7 -0.07436510 -0.01523278  0.0000000 -0.1291144  0.0000000
#> 22:    2    8 -0.07436510  0.76272457  0.0000000 -0.1291144  0.0000000
#> 23:    2   10 -0.07436510  0.00000000  0.0000000  0.4407547  0.1238637
#> 24:    2   12 -0.07436510  0.00000000  0.5980809  0.4407547 -0.4287505
#> 25:    2   13 -0.07436510  0.00000000  0.0000000  0.4407547  0.1238637
#> 26:    2   14 -0.07436510  0.00000000  0.0000000  0.4407547  0.1238637
#> 27:    2   15 -0.07436510  0.00000000 -0.1059506  0.4407547 -0.4287505
#> 28:    2   16 -0.07436510  0.00000000 -0.1059506  0.4407547 -0.4287505
#>     tree leaf   intercept  ...r=green ...g=silky ...e=broad  odor=none
#>      odor=foul  odor=anise  ...=almond  ...bruises  ...=yellow
#>  1:  0.7088975  0.00000000  0.00000000  0.00000000  0.00000000
#>  2: -0.2745896  0.00000000  0.00000000  0.00000000  0.00000000
#>  3: -0.2745896  0.00000000  0.00000000  0.00000000  0.00000000
#>  4: -0.2745896  0.00000000  0.00000000  0.00000000 -1.00780012
#>  5: -0.2745896  0.00000000  0.00000000  0.00000000  0.00000000
#>  6: -0.2745896  0.00000000 -1.02528500  0.00000000  0.06468987
#>  7: -0.2745896  0.00000000  0.00000000 -0.05144537  0.00000000
#>  8: -0.2745896  0.00000000  0.00000000  0.93244731  0.00000000
#>  9: -0.2745896  0.03208427  0.03399723  0.00000000  0.06468987
#> 10: -0.2745896 -1.04934438  0.03399723  0.00000000  0.06468987
#> 11:  0.5493453  0.00000000  0.00000000  0.00000000  0.00000000
#> 12: -0.2125787  0.00000000  0.00000000  0.00000000  0.00000000
#> 13: -0.2125787  0.00000000  0.00000000  0.00000000  0.00000000
#> 14: -0.2125787  0.00000000 -0.78526446  0.00000000  0.00000000
#> 15: -0.2125787  0.00000000  0.00000000  0.00000000  0.00000000
#> 16: -0.2125787  0.05522851  0.05139002  0.00000000  0.00000000
#> 17: -0.2125787 -0.83485874  0.05139002  0.00000000  0.00000000
#> 18: -0.2125787  0.00000000  0.00000000 -0.04271214  0.00000000
#> 19: -0.2125787  0.00000000  0.00000000  0.74060173  0.00000000
#> 20:  0.4766099  0.00000000  0.00000000  0.00000000  0.00000000
#> 21: -0.1841252  0.00000000  0.00000000  0.00000000  0.00000000
#> 22: -0.1841252  0.00000000  0.00000000  0.00000000  0.00000000
#> 23: -0.1841252  0.00000000 -0.67413428  0.00000000  0.00000000
#> 24: -0.1841252  0.00000000  0.00000000  0.00000000  0.00000000
#> 25: -0.1841252  0.04896816  0.04534281  0.00000000  0.00000000
#> 26: -0.1841252 -0.71778146  0.04534281  0.00000000  0.00000000
#> 27: -0.1841252  0.00000000  0.00000000 -0.03984976  0.00000000
#> 28: -0.1841252  0.00000000  0.00000000  0.63572177  0.00000000
#>      odor=foul  odor=anise  ...=almond  ...bruises  ...=yellow

それぞれの木*2のなかでleafにたどり着くためのルールの組み合わせと、更新される予測の量がわかる。

指定したインスタンスの予測結果を予測ルールから分解再構成する

指定したインスタンスの予測結果を、explainPredictions()が再構成した結果(breakdown)を眺める*3

# showWaterfall(xgb.model, explainer, xgb.test.data, test.data,  2, type = "binary")

DMatrix = xgb.test.data
data.matrix = test.data
idx = 2
type = "binary"
threshold = 1e-04
limits = c(NA, NA)

breakdown = explainPredictions(xgb.model, explainer, slice(DMatrix, as.integer(idx)))

ルールが多いので非ゼロのものだけ眺める。

breakdown_summary = as.matrix(breakdown)[1, ]
breakdown_summary[which(abs(breakdown_summary) > 0.0)]
#> cap-color=yellow      odor=almond       odor=anise        odor=foul 
#>       0.06468987       0.13073006       0.13628094      -0.67129347 
#>        odor=none  gill-size=broad        intercept 
#>       0.45408751       1.61423045      -0.27084657

今回は、objective = "binary:logistic"なので、予測値(対数オッズ)は逆ロジット変換により予測結果(クラス所属確率)に変換される

(weight = rowSums(breakdown))
#> [1] 1.457879

(pred = 1/(1 + exp(-weight)))
#> [1] 0.811208

分解再構成した予測結果を可視化する

プロットの方針は以下の通り。

  1. ベースライン(Intercept)を先頭に表示、残りをインパクト(貢献量)が大きい順に表示
  2. インパクトが微小(< threshold)のルールは、その他(other_impact)にまとめる
data_for_label = data.matrix[idx, ]
i = order(abs(breakdown_summary), decreasing = TRUE)
breakdown_summary = breakdown_summary[i]
data_for_label = data_for_label[i]

intercept = breakdown_summary[names(breakdown_summary) == "intercept"]
data_for_label = data_for_label[names(breakdown_summary) != "intercept"]
breakdown_summary = breakdown_summary[names(breakdown_summary) != "intercept"]

thresholdの値はオプション指定で変更でき、インパクトの大きなルールだけに絞り込んで表示できる。やたら細かいルールがいっぱい生えちゃったとかでもない限り*4thresholdはデフォルトのままで良さそう。

i_other = which(abs(breakdown_summary) < threshold)
other_impact = 0
if (length(i_other > 0)) {
  other_impact = sum(breakdown_summary[i_other])
  names(other_impact) = "other"
  breakdown_summary = breakdown_summary[-i_other]
  data_for_label = data_for_label[-i_other]
}
if (abs(other_impact) > 0) {
  breakdown_summary = c(intercept, breakdown_summary, other_impact)
  data_for_label = c("", data_for_label, "")
  labels = paste0(names(breakdown_summary), " = ", data_for_label)
  labels[1] = "intercept"
  labels[length(labels)] = "other"
} else {
  breakdown_summary = c(intercept, breakdown_summary)
  data_for_label = c("", data_for_label)
  labels = paste0(names(breakdown_summary), " = ", data_for_label)
  labels[1] = "intercept"
}

breakdown_summary
#>        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

頑張って並び替えた。

watarfall chartによる予測結果の描画

算出した値の出力。getinfoの箇所はラベルの実名か何かを取ろうとしてる? キノコのデータだとスキップされる。

if (!is.null(getinfo(DMatrix, "label"))) {
  cat("\nActual: ", getinfo(slice(DMatrix, as.integer(idx)), "label"))
}
cat("\nPrediction: ", pred)
#> 
#> Prediction:  0.811208
cat("\nWeight: ", weight)
#> 
#> Weight:  1.457879
cat("\nBreakdown")
#> 
#> Breakdown
cat("\n")
print(breakdown_summary)
#>        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

予測結果(クラス所属確率)で表示するための工夫をする。inverse_logit_labelslogitは、y軸を対数オッズからクラス所属確率に読み替えるための関数。*5

inverse_logit_trans <- scales::trans_new("inverse logit", transform = plogis, inverse = qlogis)

inverse_logit_labels = function(x) {
  return(1/(1 + exp(-x)))
}

logit = function(x) {
  return(log(x/(1 - x)))
}

ybreaks <- logit(seq(2, 98, 2)/100)

waterfalls::waterfall(
  values = breakdown_summary,
  rect_text_labels = round(breakdown_summary, 2),
  labels = labels, 
  total_rect_text = round(weight, 2),
  calc_total = TRUE,
  total_axis_text = "Prediction") + 
  scale_y_continuous(labels = inverse_logit_labels,
                     breaks = ybreaks, 
                     limits = limits) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

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

waterfalls::waterfall()は、CRANのパッケージをそのまま利用する。

cran.r-project.org

次回は、explainPredictions()が、どのような手続きで指定したインスタンスの予測結果を分解再構成しているのか詳細に見ていく。

*1:なお、この記事以降では、binary:logisticを扱って進めるが、reg:linearと本質的には変わらない。前者がわかれば後者は相対的にシンプルな手続きであり自然に理解できる。

*2: nrounds = 3 なので、決定木(0番のtree)と、対応するエラーを予測する木(1と2)

*3:インスタンスはベクトルで複数指定できるが、先頭のインスタンスだけがプロット時に使用される

*4:tree depthのチューニングに失敗してoverfitしてそうではあるし、一つのインスタンスの説明用途としては失敗しているともいえる。

*5:inverse_logit_trans`はコメントアウトしても動作する。レガシーだろうか?

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. 予測結果を分解再構成するプロセスをstep-by-stepで実行する
  4. 学習したxgboostのルール抽出を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

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