機械学習パッケージmlrを使った順序関係があるカテゴリ分類（その3）

1.はじめに

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

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

今回もその続きになります。
米国ボストン市郊外における地域別の住宅価格のデータセットを用います。
前回に引き続き、medvを3クラスのカテゴリーに分け、順序関係があるカテゴリ分類問題にします。
また、特にクラス3の分類に興味があると仮定します。
始めに必要なパッケージをインストールしておきます。
tidyverseは、今回はその中のggplotしか使っていませんので、ggplotでもOKです。

library(tidyverse)
library(ranger)
library(mlr)
library(mlbench)
library(parallel)
library(parallelMap)
data(BostonHousing, package = "mlbench")
df = BostonHousing
df$medvC = cut(df$medv,breaks =3,labels =c(1,2,3))
table(df$medvC)

$>  1   2   3 
$>215 243  48 

1.1回帰モデルを使った予測

今回は、分類タスクとしてではなく回帰タスクとして分類問題を扱います。
連続値で予測値を出力した後に、クラス間の閾値を最適化します。
クラス間の閾値は単純に四捨五入して決めるのではなく、どの値を境にクラス分けするかの区切りを求める最適化の計算を行います。

2.タスクの定義

いつものようにまずタスクを定義します。
回帰タスクなので、makeRegrTask関数を使います。

set.seed(123)
n = nrow(df)
train.set = sample(n, size = n*0.7)
test.set = setdiff(1:n, train.set)
df$medv = NULL
df$medvC=as.numeric(df$medvC)
task_BostonHousing = makeRegrTask(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: regr
$>Target: medvC
$>Observations: 354
$>Features:
$>   numerics     factors     ordered functionals 
$>         14           0           0           0 
$>Missings: FALSE
$>Has weights: FALSE
$>Has blocking: FALSE
$>Has coordinates: FALSE

3.学習器の構築

学習器はmakeLearner関数で作成します。［regr(回帰木）.ranger(ﾊﾟｯｹｰｼﾞ)を指定］
makeParamSet関数で、チューニングするハイパーパラメータの範囲を指定します。
makeTuneControlGrid関数で、グリッドサーチを指定します。
makeResampleDesc関数で交差検証法を指定します。

regr.lrn = makeLearner("regr.ranger")
ps = makeParamSet(makeIntegerParam("mtry", lower = 1, upper = 14))
ctrl = makeTuneControlGrid(resolution =  14)
rdesc = makeResampleDesc("CV",iters = 7)

3.1評価関数の作成

予測した連続値から、クラスを分ける閾値の関数(SQWK)を作成します。

SQWKfun = function(x = seq(1.5, 2.5, by = 1), pred) {
  preds = pred$data$response
  true = pred$data$truth 
  cuts = c(min(preds), x[1], x[2],  max(preds))
  preds = as.numeric(Hmisc::cut2(preds, cuts))
  err = Metrics::ScoreQuadraticWeightedKappa(preds, true, 1, 3)
  return(-err)
}

SQWK = makeMeasure(id = "SQWK", minimize = FALSE, properties = c("regr"), best = 1, worst = 0,
  fun = function(task, model, pred, feats, extra.args) {
    return(-SQWKfun(x = seq(1.5, 2.5, by = 1), pred))
  })

3.2パラメータのチューニング

評価関数にSQWKを指定して、チューニングを開始します。

library(Hmisc)
set.seed(123)
parallelStartSocket(detectCores())　　#cpuのコア数の並列処理で実行
parallelExport("SQWK", "SQWKfun")    
res4 = tuneParams(regr.lrn, task = task_BH_train, resampling = rdesc, par.set = ps, control = ctrl, measures = SQWK)
res4

$>Starting parallelization in mode=socket with cpus=8.
$>Exporting objects to slaves for mode socket: SQWK,SQWKfun
$>[Tune] Started tuning learner regr.ranger for parameter set:
$>
$>With control class: TuneControlGrid
$>Imputation value: -0
$>Exporting objects to slaves for mode socket: .mlr.slave.options
$>Mapping in parallel: mode = socket; level = mlr.tuneParams; cpus = 8; elements = 14.
$>[Tune] Result: mtry=8 : SQWK.test.mean=0.7661276
$>Tune result:
$>Op. pars: mtry=8
$>SQWK.test.mean=0.7661276

結果は、mtry=8のときにkappaが0.77となりました。
第2回の結果が0.73でしたので、少し向上しました。

3.3最適なパラメータの設定とクラス分け閾値の最適化

訓練データの正解クラスの分布を交差検証法を使って、SQWK関数を使ったクラス分けの最適な閾値を探します。

regr.lrn = setHyperPars(regr.lrn, par.vals = res4$x)
cv = crossval(regr.lrn, task_BH_train, iter = 7, measures = SQWK, show.info = TRUE)
parallelStop()
optCuts = optim(seq(1.5, 2.5, by = 1), SQWKfun, pred = cv$pred)
optCuts$par

$>Exporting objects to slaves for mode socket: .mlr.slave.options
$>Resampling: cross-validation
$>Measures:             SQWK      
$>Mapping in parallel: mode = socket; level = mlr.resample; cpus = 8; elements = 7.
$>
$>Aggregated Result: SQWK.test.mean=0.7748975
$>
$>Stopped parallelization. All cleaned up.
$>[1] 1.561646 2.552002

交差検証の結果は0.77でした。
連続値をクラスの真ん中で分けた場合は、
クラス１： -inf　-　1.5
クラス２：　1.5　-　2.5
クラス３：　2.5　-　inf
　であったのが、
クラス１： -inf　-　1.561646
クラス２：　1.561646 　-　2.552002
クラス３：　2.552002　-　inf
　に最適化されました。

4.学習の実行

それでは、準備ができたので学習を実行します。

fit4 = train(regr.lrn, task_BH_train)

5.テストデータによる評価

library(Metrics)
testTask =subsetTask(task_BostonHousing,subset =test.set )
pred = predict(fit4, testTask)
preds = as.numeric(Hmisc::cut2(pred$data$response, c(-Inf, optCuts$par, Inf)))
ScoreQuadraticWeightedKappa(preds,pred$data$truth)
table(pred$data$truth,preds)

$>[1] 0.8111621
$>   preds
$>     1  2  3
$>  1 61 10  0
$>  2 10 51  1
$>  3  0  5 14

全体の評価関数であるkappaは前回が0.69で今回が0.81と大幅に向上しました。
kappaの指標は0.6でgood、0.8ではexcellentになります。
しかし、クラス3の誤分類率は１回目の4/19から5/19、偽陽性率が1/16から1/15と少し悪くなりました。

6.参考

Prudential Life Insurance Assessment(kaggle)