mlrチュートリアルの写経（その2）

mlrチュートリアルの写経(その1)の続き．

基本的な機能（続き）

ここではモデル評価とリサンプリングを扱う．だんだんmlrの便利さが実感されてくるはず……．

モデルの評価

学習させたモデルの予測を評価するには，performance関数にPredictionオブジェクトを渡せばよい．

利用可能な評価指標

利用可能な評価指標の一覧はここで 見られるが，listMeasures関数によっても見ることができる．

よく使うのは

• auc: ROC曲線のAUC
• mmce: 平均誤分類率（mean misclassification error）
• rmse: 平方平均二乗誤差（root mean squared error）

あたりだろうか．変わったところでは，

• timetrain: 学習にかかる時間
• timepredict: 予測にかかる時間
• timeboth: 学習＋予測にかかる時間

といったものも計算できる．

# 多クラス分類に使える評価指標だけ表示
listMeasures("classif", properties = "classif.multi")

## [1] "featperc"       "mmce"           "timeboth"       "acc"           
## [5] "multiclass.auc" "ber"            "timepredict"    "timetrain"

モデル評価をやってみる

勾配ブースティングでボストン住宅価格を予測する．

n = getTaskSize(bh.task)
lrn = makeLearner("regr.gbm", n.trees = 1000)
mod = train(lrn, task = bh.task, subset = seq(1, n, 2))
pred = predict(mod, task = bh.task, subset = seq(2, n, 2))

# 予測をmseとmedse（メジアン二乗誤差）で評価（デフォルトではmseだけ出力する）
performance(pred, measures = list(mse, medse))

##       mse     medse 
## 42.773268  9.194542

上記の例でmeasuresを指定する際に，文字列ではなくてオブジェクトを渡しているが， 実は評価指標は全てMeasureクラスのオブジェクトとして実装されている（独自のMeasureを作ることも可能）．

class(mse)

## [1] "Measure"

# Measureオブジェクトの実体はfun関数
mse$fun

## function (task, model, pred, feats, extra.args) 
## {
##     measureMSE(pred$data$truth, pred$data$response)
## }
## <bytecode: 0x000000000b4d6dc0>
## <environment: namespace:mlr>

注意事項

Predictionオブジェクトだけでは計算できない評価指標もあるので注意を要する．

# 学習にかかった時間を出すにはモデルオブジェクトが必要
performance(pred, measures = timetrain, model = mod)

## timetrain 
##      0.11

# クラスタリングの評価にはTaskオブジェクトも必要なことが多い
lrn = makeLearner("cluster.kmeans", centers = 3)
mod = train(lrn, mtcars.task)
pred = predict(mod, task = mtcars.task)
performance(pred, measures = dunn, task = mtcars.task) # Dunn index

##      dunn 
## 0.2278991

2値分類の閾値調整

確率を出力する2値分類器の場合，望みの評価指標が最適な値になるように 確率の閾値を決めたいというのはよくある課題である． よくあるので，専用の便利な関数が提供されている．

# 線形判別で確率を出力させる
lrn = makeLearner("classif.lda", predict.type = "prob")
n = getTaskSize(sonar.task)
mod = train(lrn, task = sonar.task, subset = seq(1, n, by = 2))
pred = predict(mod, task = sonar.task, subset = seq(2, n, by = 2))

# 閾値の値と，指定した評価指標の値の組を作成
d = generateThreshVsPerfData(pred, measures = list(fpr, fnr, mmce))
head(d$data)

##    threshold    learner       fpr       fnr      mmce
## 1 0.00000000 prediction 1.0000000 0.0000000 0.4615385
## 2 0.01010101 prediction 0.3541667 0.1964286 0.2692308
## 3 0.02020202 prediction 0.3333333 0.2321429 0.2788462
## 4 0.03030303 prediction 0.3333333 0.2321429 0.2788462
## 5 0.04040404 prediction 0.3333333 0.2321429 0.2788462
## 6 0.05050505 prediction 0.3125000 0.2321429 0.2692308

# 閾値と各評価指標の値をプロット
plotThreshVsPerf(d)

リサンプリング

まず開発者の人によると，リサンプリング自体について知りたい人へのおすすめ文献は

• R. Simon, Resampling Strategies for Model Assessment and Selection（有料）
• B. Bischl et al., Resampling methods for meta-model validation with recommendations for evolutionary computation

らしい（特に後者の論文は主要開発者の人が筆頭著者なので，参考になるかも？）．

mlrにはモデル評価のためのリサンプリング手法が色々実装されている．

リサンプリング法はResampleDescクラスとして抽象化されている． ResampleDescオブジェクトはmakeResampleDesc関数で作成する．

提供されているリサンプリング法は以下の通り：

• クロスバリデーション (“CV”),
• Leave-one-out クロスバリデーション (“LOO”),
• くり返しクロスバリデーション (“RepCV”),
• ブートストラップ (“Bootstrap”),
• サブサンプリングクロスバリデーション (“Subsample”),
• ホールドアウト (“Holdout”)．

リサンプリングを実行するには，Task, Learner, ResampleDescオブジェクトを指定して resample関数を呼べばよい．

具体例を見よう．

# 3分割クロスバリデーションを設定
rdesc = makeResampleDesc("CV", iters = 3)
# Cox比例ハザードモデルで肺がん予測してクロスバリデーションで評価
r = resample("surv.coxph", lung.task, rdesc)

## [Resample] cross-validation iter: 1
## [Resample] cross-validation iter: 2
## [Resample] cross-validation iter: 3
## [Resample] Result: cindex.test.mean=0.611

r

## Resample Result
## Task: lung-example
## Learner: surv.coxph
## cindex.aggr: 0.61
## cindex.mean: 0.61
## cindex.sd: 0.01
## Runtime: 0.464625

resample関数はResampleResultオブジェクトを返す． ResampleResultの表示結果については説明不要だろう．

層別リサンプリング

分類問題，特に不均衡データの分類では，リサンプリングした各データにおけるクラス比率が 元のデータ全体と一致するようにしたい． 要するに正例と負例を層別してリサンプリングしたいわけだが， mlrではmakeResampleDesc関数でstratify=TRUEと指定するだけで，この層別リサンプリングができる．

また，説明変数の値で層別したい場合はstratify.cols引数に説明変数名のベクトルを指定する．

# 目的変数で層別
str(makeResampleDesc("CV", iters = 3, stratify = TRUE))

## List of 4
##  $ id      : chr "cross-validation"
##  $ iters   : int 3
##  $ predict : chr "test"
##  $ stratify: logi TRUE
##  - attr(*, "class")= chr [1:2] "CVDesc" "ResampleDesc"

# 説明変数で層別
str(makeResampleDesc("CV", iters = 3, stratify.cols = "chas"))

## List of 5
##  $ id           : chr "cross-validation"
##  $ iters        : int 3
##  $ predict      : chr "test"
##  $ stratify     : logi FALSE
##  $ stratify.cols: chr "chas"
##  - attr(*, "class")= chr [1:2] "CVDesc" "ResampleDesc"

リサンプリングで学習したモデルへのアクセス

リサンプリングの各反復中に学習させたモデルにアクセスしたい場合がある． デフォルトでは，resample関数はリサンプリング中に学習させたモデル達を捨てるようになっているので， 欲しい場合はmodels=TRUEと指定してやる必要がある．

# 判別分析で3分割クロスバリデーション
rdesc = makeResampleDesc("CV", iters = 3)
r = resample("classif.lda", iris.task, rdesc, models = TRUE)

## [Resample] cross-validation iter: 1
## [Resample] cross-validation iter: 2
## [Resample] cross-validation iter: 3
## [Resample] Result: mmce.test.mean=0.02

# クロスバリデーションで作ったモデル3つ（WrappedModelオブジェクトのリスト）
r$models

## [[1]]
## Model for learner.id=classif.lda; learner.class=classif.lda
## Trained on: task.id = iris-example; obs = 100; features = 4
## Hyperparameters: 
## 
## [[2]]
## Model for learner.id=classif.lda; learner.class=classif.lda
## Trained on: task.id = iris-example; obs = 100; features = 4
## Hyperparameters: 
## 
## [[3]]
## Model for learner.id=classif.lda; learner.class=classif.lda
## Trained on: task.id = iris-example; obs = 100; features = 4
## Hyperparameters:

さらに，リサンプリングで学習させたモデル全部はいらないが，モデルの一部の情報だけはとっておきたい場合がある． そのような場合は，extract引数にWrappedModelオブジェクトを引数とする関数を指定しておけば， リサンプリングで学習させた各モデルから必要な情報だけ引き抜いてくることができる．

# 決定木で3分割クロスバリデーション
rdesc = makeResampleDesc("CV", iters = 3)
# 各決定木の変数重要度だけ残す
r = resample("regr.rpart", bh.task, rdesc,
    extract = function(x) x$learner.model$variable.importance)

## [Resample] cross-validation iter: 1
## [Resample] cross-validation iter: 2
## [Resample] cross-validation iter: 3
## [Resample] Result: mse.test.mean=27.8

リサンプルオブジェクトの詳細

ResampleDescオブジェクトにはリサンプリング法の記述が入っているだけで， 具体的にデータがどうリサンプルされるのかは，resample関数を呼ぶまで分からない．

# ResampleDescには大したものは入っていない
rdesc = makeResampleDesc("CV", iters = 3)
str(rdesc)

## List of 4
##  $ id      : chr "cross-validation"
##  $ iters   : int 3
##  $ predict : chr "test"
##  $ stratify: logi FALSE
##  - attr(*, "class")= chr [1:2] "CVDesc" "ResampleDesc"

しかしこれだと困るときがある．例えば以下のような場合：

• 厳密に同じく分割したデータで複数のモデルを比較したい，
• （ランダムではなく）自分で訓練用データと検証用データを分割したい．

このようなことがやりたい場合，ResampleInstanceオブジェクトを作る必要がある． ResampleInstanceは実際にリサンプルされるデータ行番号を与えるオブジェクトである． ResampleInstanceオブジェクトはmakeResampleInstance関数等によって作成する．

# ResampleInstanceの作成
(rin = makeResampleInstance(rdesc, task = iris.task))

## Resample instance for 150 cases.
## Resample description: cross-validation with 3 iterations.
## Predict: test
## Stratification: FALSE

# ResampleInstanceの実体はサンプルすべき行番号のリスト 
str(rin)

## List of 5
##  $ desc      :List of 4
##   ..$ id      : chr "cross-validation"
##   ..$ iters   : int 3
##   ..$ predict : chr "test"
##   ..$ stratify: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "CVDesc" "ResampleDesc"
##  $ size      : int 150
##  $ train.inds:List of 3
##   ..$ : int [1:100] 29 134 62 87 67 108 65 128 55 14 ...
##   ..$ : int [1:100] 62 92 63 87 96 2 108 65 128 55 ...
##   ..$ : int [1:100] 29 134 92 63 96 2 67 14 84 22 ...
##  $ test.inds :List of 3
##   ..$ : int [1:50] 2 3 12 17 18 19 20 27 32 33 ...
##   ..$ : int [1:50] 1 4 6 8 9 10 11 14 22 24 ...
##   ..$ : int [1:50] 5 7 13 15 16 21 23 26 31 34 ...
##  $ group     : Factor w/ 0 levels: 
##  - attr(*, "class")= chr "ResampleInstance"

# 厳密に同じ3分割クロスバリデーションで決定木と判別分析を比較
# （今まではrdescを渡していたところがrinになっていることに注意）
r.rpart = resample("classif.rpart", iris.task, rin, show.info = FALSE)
r.lda   = resample("classif.lda",   iris.task, rin, show.info = FALSE)
c(r.rpart$aggr, r.lda$aggr)

## mmce.test.mean mmce.test.mean 
##           0.06           0.02

# makeFixedHoldoutInstance関数を使うと，訓練用データと検証用データを手で指定できる
# （この例では，150件のデータのうち先頭100件を訓練用，残りを検証用とするResampleInstanceを作成）
rin = makeFixedHoldoutInstance(train.inds = 1:100, test.inds = 101:150, size = 150)
str(rin)

## List of 5
##  $ desc      :List of 5
##   ..$ split   : num 0.667
##   ..$ id      : chr "holdout"
##   ..$ iters   : int 1
##   ..$ predict : chr "test"
##   ..$ stratify: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "HoldoutDesc" "ResampleDesc"
##  $ size      : int 150
##  $ train.inds:List of 1
##   ..$ : int [1:100] 1 2 3 4 5 6 7 8 9 10 ...
##  $ test.inds :List of 1
##   ..$ : int [1:50] 101 102 103 104 105 106 107 108 109 110 ...
##  $ group     : Factor w/ 0 levels: 
##  - attr(*, "class")= chr "ResampleInstance"

モデル評価指標値の集約方法

普通は，リサンプリングの各反復における検証データでの評価指標値を平均した値が， 最終的な評価指標値になる．しかし，評価指標値の平均ではなく，他の集約方法を設定することもできる．

実は，評価指標のMeasureオブジェクトが，リサンプリング時の集約方法をAggregationオブジェクトとして 内部に持っているので，これをsetAggregation関数で変更してやればよい（定義されているAggregationオブジェクトの一覧はヘルプ?mlr::aggregationsで見られる）．

# 誤分類率のデフォルトの集約方法は，検証用データでの平均値
mmce$aggr

## Aggregation function: test.mean

# 誤分類率の集約方法を検証用データでの中央値にしてみる
mmce2 = setAggregation(mmce, test.median)
rdesc = makeResampleDesc("CV", iters = 3)
r = resample("classif.rpart", sonar.task, rdesc, measures = list(mmce, mmce2))

## [Resample] cross-validation iter: 1
## [Resample] cross-validation iter: 2
## [Resample] cross-validation iter: 3
## [Resample] Result: mmce.test.mean=0.303,mmce.test.median=0.29

r$aggr

##   mmce.test.mean mmce.test.median 
##        0.3029676        0.2898551

もう少しおもしろい例として，B632ブートストラップを実装してみる．

# これだけだと普通のout-of-bagブートストラップ
rdesc = makeResampleDesc("Bootstrap", predict = "both", iters = 10)
# AggregationをB632に変更したMeasureを作成する
(b632.mmce = setAggregation(mmce, b632))

## Name: Mean misclassification error
## Performance measure: mmce
## Properties: classif,classif.multi,req.pred,req.truth
## Minimize: TRUE
## Best: 0; Worst: 1
## Aggregated by: b632
## Note:

# B632ブートストラップでリサンプリング！
r = resample("classif.rpart", iris.task, rdesc, measures = list(mmce, b632.mmce), show.info = FALSE)
r$aggr

## mmce.test.mean      mmce.b632 
##     0.05190563     0.04261769

便利な関数

このモデルでちょっと精度測ってみたいなー，という場合にいちいちリサンプリング法を定義するのは面倒である． そういうときのために，よく使うリサンプリング法には専用の関数が定義されている．

# 線形回帰してホールドアウトで平均二乗誤差と平均絶対誤差を算出
holdout("regr.lm", bh.task, measures = list(mse, mae))

## [Resample] holdout iter: 1
## [Resample] Result: mse.test.mean=20.9,mae.test.mean=3.47

## Resample Result
## Task: BostonHousing-example
## Learner: regr.lm
## mse.aggr: 20.92
## mse.mean: 20.92
## mse.sd: NA
## mae.aggr: 3.47
## mae.mean: 3.47
## mae.sd: NA
## Runtime: 0.019001

# 判別分析で3分割クロスバリデーションして平均誤分類率とbalanced error rateを算出
crossval("classif.lda", iris.task, iters = 3, measures = list(mmce, ber))

## [Resample] cross-validation iter: 1
## [Resample] cross-validation iter: 2
## [Resample] cross-validation iter: 3
## [Resample] Result: mmce.test.mean=0.02,ber.test.mean=0.0214

## Resample Result
## Task: iris-example
## Learner: classif.lda
## mmce.aggr: 0.02
## mmce.mean: 0.02
## mmce.sd: 0.02
## ber.aggr: 0.02
## ber.mean: 0.02
## ber.sd: 0.02
## Runtime: 0.0760038

# 上で作ったB632ブートストラップも簡単にできる
bootstrapB632("classif.lda", iris.task, measures = mmce)

## [Resample] OOB bootstrapping iter: 1
## [Resample] OOB bootstrapping iter: 2
## [Resample] OOB bootstrapping iter: 3
## [Resample] OOB bootstrapping iter: 4
## [Resample] OOB bootstrapping iter: 5
## [Resample] OOB bootstrapping iter: 6
## [Resample] OOB bootstrapping iter: 7
## [Resample] OOB bootstrapping iter: 8
## [Resample] OOB bootstrapping iter: 9
## [Resample] OOB bootstrapping iter: 10
## [Resample] OOB bootstrapping iter: 11
## [Resample] OOB bootstrapping iter: 12
## [Resample] OOB bootstrapping iter: 13
## [Resample] OOB bootstrapping iter: 14
## [Resample] OOB bootstrapping iter: 15
## [Resample] OOB bootstrapping iter: 16
## [Resample] OOB bootstrapping iter: 17
## [Resample] OOB bootstrapping iter: 18
## [Resample] OOB bootstrapping iter: 19
## [Resample] OOB bootstrapping iter: 20
## [Resample] OOB bootstrapping iter: 21
## [Resample] OOB bootstrapping iter: 22
## [Resample] OOB bootstrapping iter: 23
## [Resample] OOB bootstrapping iter: 24
## [Resample] OOB bootstrapping iter: 25
## [Resample] OOB bootstrapping iter: 26
## [Resample] OOB bootstrapping iter: 27
## [Resample] OOB bootstrapping iter: 28
## [Resample] OOB bootstrapping iter: 29
## [Resample] OOB bootstrapping iter: 30
## [Resample] Result: mmce.b632=0.0221

## Resample Result
## Task: iris-example
## Learner: classif.lda
## mmce.aggr: 0.02
## mmce.mean: 0.02
## mmce.sd: 0.02
## Runtime: 0.70104

ここまでのまとめ

• mlrではモデルの評価とリサンプリングのための統一されたインターフェイスが提供されている
• 組み込みの評価指標やリサンプリング法の種類も充実している

mlrチュートリアルの写経（その3）に続く（予定）……