tidyポエムうしろのtidy

いくつか作ったモデルの予測の精度を比較したい時、何をするのがよいのだろうか。バリデーション用に取っておいたデータへのあてはまりの良さを使ってモデルを比べていくことが多いと思うが、いずれかのモデルが偶然良かったからと言って、そのモデルが本当に良かったと言って良いのだろうか。
caretというパッケージがある。Rで機械学習のパラメータチューニングを、統一されたフレームワークで実行できるパッケージだ。このパッケージには先の疑問に対する1つの答えが関数として用意されている。rsample()関数である。これを使って実際にモデルを比較してみよう。まず適当にモデルを2つ作る。

library(caret)
library(mlbench)
library(dplyr)
library(tidyr)
data(BostonHousing)

BostonHousing <- na.omit(BostonHousing)

set.seed(1)
tr = trainControl(
  method = "cv",
  number = 10,
  savePredictions = T)

train_grid_rf = expand.grid(mtry = c(1:13), 
                            splitrule = "variance") 
set.seed(3)
rf_res = train(BostonHousing[, c(1:13)],  
               BostonHousing$medv,　　 
               method = "ranger",　　 
               tuneGrid = train_grid_rf, 
               trControl=tr,                
               preProc = c("center", "scale")
               )

train_grid_glm = expand.grid(alpha = 10 ^ (0:10 * -1), lambda = 10 ^ (0:10 * -1)) 

set.seed(3)
enet_res = train(data.matrix(BostonHousing[, c(1:13)]),  
                 BostonHousing$medv,　　 
                 method = "glmnet",　　 
                 tuneGrid = train_grid_glm, 
                 trControl=tr,                
                 preProc = c("center", "scale")
                 )

その後、得られた最適化済みモデルの各foldにおけるMAE, RMSE, Rsquaredなどを以下コードの通り抽出できる。判別分析であればROCなど、また別の評価指標を抽出することももちろん可能だ。

resamps <- resamples(list(rf = rf_res,
                          enet = enet_res))
summary(resamps)

## 
## Call:
## summary.resamples(object = resamps)
## 
## Models: rf, enet 
## Number of resamples: 10 
## 
## MAE 
##          Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## rf   1.608598 1.940034 2.050975 2.075349 2.217910 2.648400    0
## enet 2.958318 3.180185 3.382017 3.375692 3.598527 3.793749    0
## 
## RMSE 
##          Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## rf   1.895643 2.569456 2.938506 2.965669 3.418101 4.042882    0
## enet 3.813264 4.310452 4.859145 4.840166 5.323722 6.075280    0
## 
## Rsquared 
##           Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
## rf   0.8540058 0.8726932 0.8886106 0.8954780 0.9066470 0.9651061    0
## enet 0.5840123 0.6937885 0.7247992 0.7279925 0.7737212 0.8255980    0

図にするとこんな感じだろうか。

bwplot(resamps, layout = c(3, 1),
       scales = list(relation = "free"),
       xlim = list(c(1, 5), c(0, 6),  c(0.4, 1)))

これはこれで良いのだろうが、caret作者のMax Kuhn氏はこれをもうちょっとパワーアップしたかったようだ。そこでBayesian glmを使った tidyposterior なるパッケージを作成された。ようやくtidyポエムなのに初めてtidyって言う単語が出てきた。これはfoldごとの検証値が含まれていれば、さっき作ったrsampleで出力されたリストにも普通のdata.frameにも使えるとのことだ。

とりあえずどのように使うかをresampsリストのRMSEを例にやってみる。

library(tidyposterior)
RMSEs <- resamps$values %>%
  select(Resample, contains("RMSE"))

colnames(RMSEs) <- c("id", "rf", "enet")

RMSEs

##        id       rf     enet
## 1  Fold01 2.156541 4.483690
## 2  Fold02 3.429724 5.532127
## 3  Fold03 2.955112 3.986372
## 4  Fold04 2.490247 3.813264
## 5  Fold05 2.807083 5.162452
## 6  Fold06 4.042882 5.051743
## 7  Fold07 3.383231 4.666547
## 8  Fold08 3.574327 6.075280
## 9  Fold09 2.921900 5.377479
## 10 Fold10 1.895643 4.252705

broomがDependsになってる。

で、パッケージのperf_mod()関数に作ったここで作ったデータフレームRMSEsを入れてやるとどこかで見たような CHAINが4回走って事後分布を計算してくれる。

rmse_model <- perf_mod(RMSEs, seed = 71)

## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## 
## Gradient evaluation took 0.000103 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 1.03 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.391174 seconds (Warm-up)
##                0.250606 seconds (Sampling)
##                0.64178 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## 
## Gradient evaluation took 2.5e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.314732 seconds (Warm-up)
##                0.209008 seconds (Sampling)
##                0.52374 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## 
## Gradient evaluation took 2.4e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.349441 seconds (Warm-up)
##                0.20834 seconds (Sampling)
##                0.557781 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## 
## Gradient evaluation took 1.8e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.343246 seconds (Warm-up)
##                0.324269 seconds (Sampling)
##                0.667515 seconds (Total)

summaryの格納先はstanでは正直なんだなという気持ちがある。

rmse_model$stan

## stan_glmer
##  family:  gaussian [identity]
##  formula: statistic ~ model + (1 | id)
## ------
## 
## Estimates:
##             Median MAD_SD
## (Intercept) 3.0    0.2   
## modelenet   1.9    0.2   
## sigma       0.5    0.1   
## 
## Error terms:
##  Groups   Name        Std.Dev.
##  id       (Intercept) 0.52    
##  Residual             0.55    
## Num. levels: id 10 
## 
## Sample avg. posterior predictive 
## distribution of y (X = xbar):
##          Median MAD_SD
## mean_PPD 3.9    0.2   
## 
## ------
## For info on the priors used see help('prior_summary.stanreg').

ここから事後分布を取り出すにはtidy()関数を使う。

rmse_model_post <- tidy(rmse_model)
str(rmse_model_post)

## Classes 'posterior', 'tbl_df', 'tbl' and 'data.frame':   8000 obs. of  2 variables:
##  $ model    : chr  "rf" "rf" "rf" "rf" ...
##  $ posterior: num  2.72 2.78 2.81 2.86 2.64 ...

1000sampleing4chains2modelで8000個のデータが格納されている。

図にするとこんな感じ。

ggplot(rmse_model_post) + 
  # Add the observed data to check for consistency 
  geom_point(
    data = rmse_model_post, 
    aes(x = model, y = posterior), 
    alpha = .5, col = "blue"
  )

まあどう見てもrfのほうがrmse小さいよねという感じだが、ちゃんと比較するための関数contrast_models()も用意されている。

rf_v_enet <- contrast_models(rmse_model, "rf", "enet")
str(rf_v_enet)

## Classes 'posterior_diff' and 'data.frame':   4000 obs. of  3 variables:
##  $ difference: num  -1.47 -2.16 -1.8 -1.73 -2.03 ...
##  $ model_1   : chr  "rf" "rf" "rf" "rf" ...
##  $ model_2   : chr  "enet" "enet" "enet" "enet" ...

strで中身を視るとrfとenetの差分を作ってるっぽい。

summary(rf_v_enet, size = 1)

## # A tibble: 1 x 9
##     contrast probability      mean     lower     upper  size pract_neg
##        <chr>       <dbl>     <dbl>     <dbl>     <dbl> <dbl>     <dbl>
## 1 rf vs enet           0 -1.864987 -2.280377 -1.450151     1    0.9975
## # ... with 2 more variables: pract_equiv <dbl>, pract_pos <dbl>

平均の差、差の上限、下限とかを出してくれる。sizeは差の真の大きさの予想で、positive, negative方向にどのくらいの確率でずれてるのか、もしくは変わらないかがpract_neg, pract_equiv, pract_posに表示されているようだ。

ggplot(rf_v_enet, size = 1)

図でもこれを出すことができる。 tidyって名前がついてたからtidyなんだろうと思ってちょっと触ってみたけどtidyなのかどうかわからなかった。そもそもtidyがどんな感じか私にはよくわからない。うんうん、それもまたtidyだね。KaggleとかでモデルのLocal CV見るのには使えるのかもしれない。

tidyポエムうしろのtidy

@siero5335

2017/12/15

tidyポエム うしろのtidy

@siero5335

2017/12/15

tidyポエムうしろのtidy