##1.はじめに
仕事では順序関係があるカテゴリ分類のデータを扱うことが多いので、Rのmlrパッケージを使った機械学習にチャレンジしてみました。
ちなみに、仕事はIT関係では全くないので、全くの素人分析屋です。
前回は、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層化抽出法を使ったK-分割交差検証法 (Stratified K-Fold CV)
クラス３だけが少ない不均衡データですので、ハイパーパラメータのチューニングに層化抽出法を使ったK-分割交差検証法を使ってみたいと思います。
この方法は交差検証法を使うときに、foldごとに含まれるクラスの割合を等しくする方法です。
mlrではmakeResampleDesc関数の引数stratify = TRUEとすれば簡単にできます。
ではtrain data とtest dataに分けます。
7割をtrain dataとして残りをtest data にします。
前回と比較するために、同じデータセットを読み込みます。
n = nrow(df)
train.set = scan("train.txt") #前回の結果を読む
test.set = setdiff(1:n, train.set)
df$medv=NULL
df$medvC=as.factor(df$medvC)
task_BostonHousing = makeClassifTask(id = "BostonHousing", data = df, target = "medvC")
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
$> 12 1 0 0
$>Missings: FALSE
$>Has weights: FALSE
$>Has blocking: FALSE
$>Has coordinates: FALSE
$>Classes: 3
$> 1 2 3
$>144 181 29
$>Positive class: NA
##2.学習器からハイパーパラメータのチューニングまで
前回と同様にランダムフォレスト(ranger)を使います。
今回は学習器の構築からチューニングまで前回と同様なので一気にやります。
set.seed(123)
parallelStartSocket(detectCores())
classif.lrn = makeLearner("classif.ranger",importance="permutation")
ps = makeParamSet(makeIntegerParam("mtry", lower = 1, upper = 13))
ctrl = makeTuneControlGrid(resolution = 13)
rdesc = makeResampleDesc("CV",iters = 7, stratify = TRUE)
res2 = tuneParams(classif.lrn, task = task_BH_train, resampling = rdesc, par.set = ps, control = ctrl,measures = kappa)
parallelStop()
res2
$>Starting parallelization in mode=socket with cpus=8.
$>[Tune] Started tuning learner classif.ranger for parameter set:
$>With control class: TuneControlGrid
$>Imputation value: 1
$>Exporting objects to slaves for mode socket: .mlr.slave.options
$>Mapping in parallel: mode = socket; level = mlr.tuneParams; cpus = 8; elements = 13.
$>[Tune] Result: mtry=4 : kappa.test.mean=0.7302055
$>Stopped parallelization. All cleaned up.
$>Tune result:
$>Op. pars: mtry=4
$>kappa.test.mean=0.7302055
結果は、mtry=4のときにkappaが0.73となりました。
前回の結果が0.70でしたので、少し向上しました。
早速、fitさせて、テストデータで評価します。
##3.訓練データによるモデルの作成とテストデータによる評価
classif.lrn = setHyperPars(classif.lrn, par.vals = res2$x)
fit2 = train(classif.lrn, task_BH_train)
testTask =subsetTask(task_BostonHousing ,subset =test.set )
pred = predict(fit2, testTask)
performance(pred,measures = kappa)
calculateConfusionMatrix(pred,sums = T)
$> kappa
$>0.6810189
$> 1 2 3 -err.- -n-
$>1 54 17 0 17 71
$>2 7 54 1 8 62
$>3 0 4 15 4 19
$>-err.- 7 21 1 29 NA
$>-n- 61 75 16 NA 354
前回がkappa0.69で今回が0.68と、ほとんど同じです。
クラス３誤分類率や偽陽性率も同じでした。
今回のデータでは層化抽出法を使ったK-分割交差検証によるハイパーパラメータのチューニングはあまり効果は見られませんでした。
##4.重み付け
今度はさらに、クラスに重み付け(weights)を設定したいと思います。
w=c(1,1,10)[df$medvC]
wightの定義です。値はクラス1,2,3にそれぞれ1,1,10を充てた（値は適当に振り分けています）。
次に重み付けしたtrain_taskを設定します。
df$medv = NULL
df$medvC = as.factor(df$medvC)
task_BostonHousing = makeClassifTask(id = "BostonHousing", data = df, target = "medvC", weights = w)
task_BH_train = subsetTask(task_BostonHousing ,subset =train.set )
task_BH_train
$>Supervised task: BostonHousing
$>Type: classif
$>Target: medvC
$>Observations: 354
$>Features:
$> numerics factors ordered functionals
$> 12 1 0 0
$>Missings: FALSE
$>Has weights: TRUE
$>Has blocking: FALSE
$>Has coordinates: FALSE
$>Classes: 3
$> 1 2 3
$>144 181 29
$>Positive class: NA
makeClassifTask関数において引数にweights = wを指定します。
これで目的変数に重み付けをすることができます。説明にはHas weights: TRUEとなっています。
##5.学習器からハイパーパラメータのチューニングまで
前回までと同様なので一気にやります。
set.seed(123)
parallelStartSocket(detectCores())
classif.lrn = makeLearner( "classif.ranger",importance="permutation")
ps = makeParamSet( makeIntegerParam("mtry", lower = 1, upper = 13))
ctrl = makeTuneControlGrid(resolution = 13)
rdesc = makeResampleDesc("CV",iters = 7, stratify = TRUE)
res3 = tuneParams(classif.lrn, task = task_BH_train, resampling = rdesc, par.set = ps, control = ctrl,measures = kappa)
parallelStop()
res3
$>Starting parallelization in mode=socket with cpus=8.
$>[Tune] Started tuning learner classif.ranger for parameter set:
$>With control class: TuneControlGrid
$>Imputation value: 1
$>Exporting objects to slaves for mode socket: .mlr.slave.options
$>Mapping in parallel: mode = socket; level = mlr.tuneParams; cpus = 8; elements = 13.
$>[Tune] Result: mtry=5 : kappa.test.mean=0.7221793
$>Stopped parallelization. All cleaned up.
$>Tune result:
$>Op. pars: mtry=5
$>kappa.test.mean=0.7221793
結果は、mtry=5のときにkappaが0.72となりました。
早速、fitさせて、テストデータで評価します。
##6.訓練データによるモデルの作成とテストデータによる評価
classif.lrn = setHyperPars(classif.lrn, par.vals = res2$x)
fit3 = train(classif.lrn, task_BH_train)
testTask = subsetTask(task_BostonHousing ,subset =test.set )
pred = predict(fit3, testTask)
performance(pred,measures = kappa)
calculateConfusionMatrix(pred,sums = T)
$> kappa
$>0.6644114
$> 1 2 3 -err.- -n-
$>1 56 15 0 15 71
$>2 9 48 5 14 62
$>3 0 2 17 2 19
$>-err.- 9 17 5 31 NA
$>-n- 65 65 22 NA 354
weight（重み付け）した結果は、全体の評価関数であるkappaは前回が0.69で今回が0.66とほぼ同等でした。
次にクラス3の誤分類率は前回の4/19から2/19と向上しましたが、偽陽性率が前回の1/16から5/22と上昇しました。
偽陽性率と真陽性率（誤分類率）はトレードオフの関係があるので、重み付けすると誤分類率は下がりますが、誤って予測する確率が上がります。
ここではクラス３に最も関心があるという設定ですので、クラス3の予測が少しぐらい外れても、取りこぼしのないように広く予測をしたいのであれば、重み付けしたモデルを、取りこぼしがあっても予測が外れる確率を減らしたいのであれば、重み付けはしないモデルを選ぶといった感じでしょうか。
また、ランダムフォレストは変数の重要度を計算できますので、予測に用いた説明変数の重要度を見てみます。
重み付けしたモデル(fit3)の重要度を見てみます。
vmp=data.frame(getFeatureImportance(fit3)[1])
ggplot(vmp,aes(x=res.variable,y=res.importance))+geom_bar(stat = "identity")+coord_flip()
重要度の高い説明変数は、Istatと、ageとなりました。
次にPermutationベースの変数重要度を指定したので、Partial Dependence Plotを見てみます。
Partial Dependence Plotとは、以下のとおりです。
#Partial Dependence Plot
各変数の重要度がわかったら、次に行うべきは重要な変数とアウトカムの関係を見ることだと思います。 ただ、一般にブラックボックスモデルにおいてインプットとアウトカムの関係は非常に複雑で、可視化することは困難です。
そこで、複雑な関係を要約する手法の一つにPartial Dependence Plot(PDP)があります。PDPは興味のある変数以外の影響を周辺化して消してしまうことで、インプットとアウトカムの関係を単純化しようというものです。
Istat（人口に下層クラスが占める割合）の説明変数のPDPを見てみます。
lrn = makeLearner("classif.ranger", predict.type = "prob",mtry=5)
fit3.1 = train(lrn,task_BH_train)
pd = generatePartialDependenceData(fit3.1, task_BH_train,"lstat")
plotPartialDependence(pd)
PDPを実行するには、学習器のpredict.type を確率に変更する必要があるので、1行目で学習器を定義しなおしました。
人口に下層クラスが占める割合（Istat）が20%以上になると、クラス１の確率が最も高くなりました。
住宅価格には下層クラスの割合が最も重要なファクターのようです。
さて、次回もさらにモデルの検討を進めてみます。