##1.はじめに
仕事では順序関係があるカテゴリ分類のデータを扱うことが多いので、Rのmlrパッケージを使った機械学習にチャレンジしてみました。
ちなみに、仕事はIT関係では全くないので、全くの素人分析屋です。

前回までは、mlrを使ってランダムフォレストによるモデルを作成し、予測してみました。
機械学習パッケージmlrを使った順序関係があるカテゴリ分類（その１）
機械学習パッケージmlrを使った順序関係があるカテゴリ分類（その２）
機械学習パッケージmlrを使った順序関係があるカテゴリ分類（その３）

今回は、勾配ブースティング(xgboost)によるモデルを使ってみたいと思います。
xgboostはKaggleでよく用いられるアルゴリズムで、上位入賞者によく用いられている高性能なアルゴリズムで有名です。

##2.データとタスクの定義
早速、タスクを定義していきます。
まずは、分類問題としてやってみます。
xgboostのクラス分類は０から始める必要があるため、1を引きます。

library(tidyverse)
library(xgboost)
library(mlr)
library(mlbench)
library(parallel)　#detectCores()を使ってCPUコア数を取得する
library(parallelMap)
data(BostonHousing, package = "mlbench")
df = BostonHousing
df$medvC = cut(df$medv,breaks =3,labels =c(1,2,3))
n = nrow(df)
train.set = scan("train.txt")　　　　　#第1回目で保存したもの
test.set = setdiff(1:n, train.set)

df$medv = NULL
df$medvC = as.numeric(df$medvC)-1　　#クラスを0baseにする
df$medvC = as.factor(df$medvC)       #目的変数のfactor化
task_BostonHousing = makeClassifTask(id = "BostonHousing", data = df, target = "medvC")
task_BostonHousing = createDummyFeatures(task_BostonHousing)
task_BH_train =subsetTask(task_BostonHousing ,subset =train.set )
task_BH_train

$>Read 354 items
$>Supervised task: BostonHousing
$>Type: classif
$>Target: medvC
$>Observations: 354
$>Features:
$>   numerics     factors     ordered functionals 
$>         14           0           0           0 
$>Missings: FALSE
$>Has weights: FALSE
$>Has blocking: FALSE
$>Has coordinates: FALSE
$>Classes: 3
$>  0   1   2 
$>144 181  29 
$>Positive class: NA

##3.学習器の構築
objective [default=reg:linear]
最小化する損失関数を定義します。主に使用される関数は次のとおりです。
・binary:logistic ：バイナリ分類のためのロジスティック回帰、予測された確率（クラスではなく）を返します。
・multi:softmax ：マルチクラス分類，予測されたクラス（確率ではなく）を返します。また、クラスの数を定義する num_class (クラス数)パラメータを追加で設定します。
・multi:softprob ：softmaxと同じですが、各データポイントが各クラスに属する予測確率を返します。

eval_metric [ default according to objective ]
デフォルト値は、回帰の場合はrmse、分類の場合はerrorです。
代表的な値は以下の通りです。
・rmse ：二乗平均平方根誤差　　　　　　　　　［回帰］
・mae ：平均絶対誤差　　　　　　　　　　　　　［回帰］
・logloss ：ログロス（交差エントロピー）　　　　［分類］
・error ：二値分類エラー率（0.5：閾値）　　　　［分類］
・merror ：多クラス分類エラー率　　　　　　　　［分類］
・mlogloss ：マルチクラス-ログロス　　　　　　　［分類］
・auc：曲線下面積　　　　　　　　　　　　　　　［分類］

makeLearner関数は、分類には"classif.xgboost"を指定します。

set.seed(123)
parallelStartSocket(detectCores())
classif.lrn = makeLearner("classif.xgboost",
                          booster = "gbtree",
                          objective = "multi:softmax",
                          eval_metric ="mlogloss")

$>Starting parallelization in mode=socket with cpus=8.

##3.1分類問題にXgboostを使う場合の注意点
Xgboostは通常、factorを扱えないため、目的変数であるクラスはintegerに変換する必要があるが、mlrのmakeClassifTask関数は目的変数がfactorであることが必要なため、Xgboostであっても目的変数はfactorにします。

##4.ハイパーパラメータの調整
4.1 General Parameters: 全体的な機能
　gblinear、gbtreeの2種類がありますが、gbtreeを選びます（決定木ベース）
4.2 Booster Parameters:
　分類木か回帰木を選びます
4.3　Learning Task Parameters:
　実行状況において最適化します

4.4 最適化するParametersの種類
eta [default=0.3]
・いわゆる学習率
・各ステップの重みを縮小することで，モデルをよりロバストなものにします
・使用される代表的な値　0.01-0.2

min_child_weight [default=1]
・過学習を制御するために使用されます。高い値は、選択された特定のsampleに非常に特異的になる関係を学習するのを防ぎます。しかし、値が高すぎると未学習につながるので交差検証法等を用いて調整します

max_depth [default=6]
・木の最大深さ
深いほど，モデルが特定のsampleに特異的な関係を学習できるので，過学習を制御するために利用されます
・交差検証法等を用いて調整します
・標準的な値　3-10

gamma [default=0]
・ノードが分割されるのは、結果として生じる分割が損失関数の正の減少をもたらす場合のみです。Gammaは分割に必要な最小損失削減量を指定します。この値は損失関数に依存して変化することがあり、調整する必要があります。

subsample [default=1]
・各木についてランダムにサンプルされる観測値の割合を示します。低い値は，アルゴリズムをより保守的にし，過学習を防ぎますが，小さすぎる値は未学習につながるかもしれません。
・標準的な値　0.5-1

colsample_bytree [default=1]
・各木に対してランダムにサンプルされるコラムの割合を示します。
・標準的な値　0.5-1

lambda [default=1]
・重みのL2正則化項。過学習を減らすために検討します。

alpha [default=0]
・重みのL1正則化項。非常に高次元の場合に使用され、実装時にアルゴリズムが高速に実行されます。

nrounds
・nroundsで繰り返し計算回数（木の数）を指定します。

調整するハイパーパラメータの範囲を設定します。
リサンプリングは、層化抽出法を使った交差検証法を用います。

ps = makeParamSet(
  makeNumericParam("eta", lower = 0.01, upper = 0.3),
  makeNumericParam("colsample_bytree", lower = 1, upper = 2, trafo = function(x) x/2),
  makeIntegerParam("min_child_weight", lower = 1, upper = 20),
  makeIntegerParam("nrounds", lower = 10,  upper = 300),
  makeIntegerParam("max_depth", lower = 1, upper = 8)
  )
rdesc = makeResampleDesc("CV",iters = 7, stratify = TRUE)

#4.1ベイズ最適化を用いたハイパーパラメータのチューニング
チューニングの方法は、前回まで一番単純なグリッド　サーチ法でしていましたが、パラメータが多くなると指数関数的に計算量が増加してしまいます。
そこで、今回はベイズ最適化の方法を使います。ガウス過程の回帰によりパラメータの最小値を逐次的に探す方法で、多数のパラメータの組み合わせを効率的に探すことができます。
参考：ベイズ最適化でハイパーパラメータを調整する
パッケージはmlrMBOを使います

library(mlrMBO)
mbo.ctrl = makeMBOControl()
mbo.ctrl = setMBOControlTermination(mbo.ctrl, iters = 5)
ctrl = makeTuneControlMBO( mbo.control = mbo.ctrl)

ベイズ最適化を実行します。

res5 = tuneParams(learner = classif.lrn, 
                     task = task_BH_train, 
                     resampling =rdesc , 
                      par.set = ps, 
                      control = ctrl, 
                      show.info = TRUE, 
                     measures = kappa)

res5
parallelStop()

$>[Tune] Started tuning learner classif.xgboost for $>parameter set:
$>
$>With control class: TuneControlMBO
$>Imputation value: 1
$>Mapping in parallel: mode = socket; level = $>mlrMBO.feval; cpus = 8; elements = 20.
$>Mapping in parallel: mode = socket; level = $>mlrMBO.feval; cpus = 8; elements = 1.
$>Mapping in parallel: mode = socket; level = $>mlrMBO.feval; cpus = 8; elements = 1.
$>Mapping in parallel: mode = socket; level = $>mlrMBO.feval; cpus = 8; elements = 1.
$>Stopped because hard maximum generation limit was $>hit.Mapping in parallel: mode = socket; level = $>mlrMBO.feval; cpus = 8; elements = 1.
$>Mapping in parallel: mode = socket; level = $>mlrMBO.feval; cpus = 8; elements = 1.
$>[Tune] Result: eta=0.108; colsample_bytree=0.957; $>min_child_weight=2; nrounds=147; max_depth=4 : $>kappa.test.mean=0.7108326
$>Tune result:
$>Op. pars: eta=0.108; colsample_bytree=0.957; $>min_child_weight=2; nrounds=147; max_depth=4
$>kappa.test.mean=0.7108326
$>Stopped parallelization. All cleaned up.

実行には少し時間がかかります。結果は、eta=0.103; colsample_bytree=0.957; min_child_weight=2; nrounds=147; max_depth=4
となり、kappa=0.71となりました。

早速、この結果を学習器にセットします。そしてセットしたパラメータを見てみます。

classif.lrn = setHyperPars(classif.lrn, par.vals = res5$x)
classif.lrn$par.vals

$>$nrounds
$>[1] 147
$>$verbose
$>[1] 0
$>$booster
$>[1] "gbtree"
$>$objective
$>[1] "multi:softmax"
$>$eval_metric
$>[1] "mlogloss"
$>$eta
$>[1] 0.1082631
$>$colsample_bytree
$>[1] 0.9569314
$>$min_child_weight
$>[1] 2
$>$max_depth
$>[1] 4

それでは、fitさせて、テストデータで評価します。

##5.訓練データによるモデルの作成とテストデータによる評価

fit5 = train(classif.lrn, task_BH_train)
testTask =subsetTask(task_BostonHousing ,subset =test.set )
pred = predict(fit5, testTask)
performance(pred,measures = kappa)
calculateConfusionMatrix(pred,sums = T)

$>    kappa 
$>0.7141204 
$>        0  1  2 -err.- -n-
$>0      60 11  0     11  71
$>1      10 50  2     12  62
$>2       0  3 16      3  19
$>-err.- 10 14  2     26  NA
$>-n-    70 64 18     NA 354

第1回のkappa0.69よりは少し向上して0.71、クラス３の誤分類率は一番良かった第2回の2/19よりは少し悪い3/19の誤分類率となりました。
偽陽性率は一番良かった第1回の1/16には及びませんが2/18となったので、クラス３に関してはバランスのとれた結果となりました。

今回のデータではXgboostは全体的な評価として、ランダムフォレストとほぼ同等な感じでしょうか。
もしかするともう少しパラメーターチューニングをすれば性能が上がるのかもしれません。

##6.autoxgboostを使ってみる
autoxgboostは、mlr とmlrMBOを使ってxgboostのモデルをベイズ最適化で自動的にチューニングするためのラッパーだそうです。
これを使って上記と比較してみたいと思います。

'CRAN'に入っていないようなので、githubからinstallします。

devtools::install_github("ja-thomas/autoxgboost")

インストールが終わりましたら、MBOコントロールを設定します。

library(autoxgboost)
set.seed(123)
parallelStartSocket(detectCores())
ctrl = makeMBOControl()
ctrl = setMBOControlTermination(ctrl, iters = 3L)

ハイパーパラメーターの自動チューニングを実行します。

res.auto = autoxgboost(task_BH_train, design.size = 15L, control = ctrl, measure = kappa, nthread = 4, tune.threshold = FALSE)
parallelStop()

$>　省略

自動チューニングされた結果を表示します。

res.auto

$>Autoxgboost tuning result
$>Recommended parameters:
$>              eta: 0.199
$>            gamma: 0.038
$>        max_depth: 14
$> colsample_bytree: 0.984
$>colsample_bylevel: 0.716
$>           lambda: 1.246
$>            alpha: 0.003
$>        subsample: 0.561
$>          nrounds: 9
$>
$>Preprocessing pipeline:
$>(fixfactors >> dummyencode >> dropconst)(fixfactors.drop.unused.levels = TRUE, fixfactors.fix.factors.prediction = $>TRUE, dummyencode.reference.cat = FALSE, dummyencode.infixdot = TRUE, dropconst.rel.tol = 1e-08, dropconst.abs.tol = 1e-08, dropconst.ignore.na = FALSE)
$>
$>With tuning result: kappa = 0.787

Recommended parametersは、訓練データで kappa = 0.79と非常に良い結果が出ています。
早速、テストデータで評価します。

testTask =subsetTask(task_BostonHousing ,subset =test.set )
pred = predict(res.auto$final.model, testTask)
performance(pred,measures = kappa)
calculateConfusionMatrix(pred,sums = T)

$>    kappa 
$>0.6233785 
$>        0  1  2 -err.- -n-
$>0      55 16  0     16  71
$>1      12 49  1     13  62
$>2       0  5 14      5  19
$>-err.- 12 21  1     34  NA
$>-n-    67 70 15     NA 354

テストデータではkappa0.62と訓練データを使った時よりも悪くなりました。
モデルは過学習している傾向が見られます。また、クラス３の誤分類率も5/19と悪くなりました。
残念ながら自動チューニングは、あまり精度は高くでませんでした。

##7.参考
mlrMBO: Quick introduction
autoxgboost を使ってみる