モデルスタッキング時の{stacks}と{finetune}の不健全な関係


# はじめに

この記事は R Advent Calendar 2021の１６日目の記事です。

https://qiita.com/advent-calendar/2021/rlang


日本最大級のRコミュニティであるJapan.Rが2021年度も無事開催され、Rに関する様々な知見が共有されておりました。
かくいう私も昨年に引き続き、tidymodelsについて稚拙なライトニングトークをさせていただきました。

https://japanr.connpass.com/event/230381/

今回はtidymodelsについて抜粋してパッケージを紹介し、より皆様にtidymodelsを、そしてR言語を知っていただければ幸いです。

japan.Rでtidymodelsのお話をしたので、
その一部分を切り出して深掘りご紹介。

# 免責(言い訳)

「ブレンディングとスタッキングはちゃうんや」、というご意見が出そうですが、厳密な定義が見当たらず、kaggle等でも口語的に使い分けられていることから今回のように複数モデルの出力結果の組み合わせから新しくモデルを構築する方法を「スタッキング」と呼ばせていただきます。
ちなみにtidymodelsのstacksパッケージの公式の説明では「ブレンディング」という単語が使われているように感じます。コメント頂けるとその内容によっては当記事の表現も変更する準備がありますので、ご意見お待ちしております。

※本記事はtidymodelsの関数や機械学習のチューニングについてある程度ご存じの方に向けて書いておりますので、機械学習初心者やR初心者には不親切かもしれません。


# パラメータチューニング


機械学習のタスクに向き合っているとき、モデルの予測精度を上げようとパラメータをチューニングする機会に出会います。

tidymodelsではparsnipパッケージのアルゴリズムを指定する関数にtuneパッケージのプレースホルダを指定することで、チューニングされるべき対象を用意しておくことができます。
その後、tune_*からはじまる関数によってチューニングパラメタを変更しながら予測精度を出力し、どのようなパラメタの組み合わせが最も予測精度を上げるのかを得ることができます。


```{r}
library(pacman)
p_load(tidymodels,xgboost)

data(ames)
ames$Sale_Price <- log(ames$Sale_Price)

df <- ames

split_df <- initial_split(df,prop = 0.8,strata = Sale_Price)

ames %>% 
  glimpse()

train <- training(split_df)
test <- testing(split_df)

xgboost_recipe <- 
  recipe(formula = Sale_Price ~ ., data = train) %>% 
  step_novel(all_nominal(), -all_outcomes()) %>% 
  step_dummy(all_nominal(), -all_outcomes(), one_hot = TRUE) %>% 
  step_zv(all_predictors()) 

xgboost_spec <- 
  boost_tree(trees = tune(), min_n = tune(), 
             tree_depth = tune(), learn_rate = tune(), 
             loss_reduction = tune(), sample_size = tune()) %>% 
  set_mode("regression") %>% 
  set_engine("xgboost") 

xgboost_workflow <- 
  workflow() %>% 
  add_recipe(xgboost_recipe) %>% 
  add_model(xgboost_spec) 

splits <- train %>% vfold_cv(v=5,strata = Sale_Price)
set.seed(40509)
xgboost_tune <-
  tune_grid(xgboost_workflow, 
            resamples = splits, 
            grid = 5,
            verbose=T,
            control = control_grid(save_pred = T,save_workflow = T))
```


チューニングパラメタの中で最も性能が良いモデルを確かめるには、


```{r}
xgboost_tune %>% 
  autoplot(metric = "rmse")
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/f6ef0f69-3ddb-d0ba-9903-b9cea57a661e.png)


```{r}
xgboost_tune %>% 
  show_best(metric = "rmse")
```



# モデルスタッキング


この後、最も良いパラメタの組み合わせによって予測モデルを作成することになるのですが、**本当に最も良いパラメタだけ使えばよいのでしょうか？**

もしかしたらデータの都合上**偶然性能が良かった**可能性もあります。
交差検証はそのリスクに対応するため設計されたものですし、二重クロスバリデーションのようなテクニックを使い交差検証を更にレベルアップさせる選択もあります。
しかし、データのレコード数が少なくテクニックが使えない場合もあります。
そのような場合には、複数のモデルを組み合わせ、モデル同士の弱い領域をカバーしあうモデルを構築してあげます。

歴史的にアンサンブル学習によって予測精度が上がることが知られているので、出来上がったモデルをアンサンブルするのは自然なことでしょう。
実際にデータ分析コンペティションであるkaggle等でその性能の向上が見られています。
※特に似た理論によるアルゴリズムよりも、異なる理論によるアルゴリズムを組み合わせた方が性能が上がりやすいとされています。

このような考え方によって「モデルのスタッキング」が生まれました。

それではtidymodelsのstacksパッケージで、チューニングする際に作成される個別パラメタモデルをアンサンブルしてみましょう。
stacksではLASSOによって予測結果に重みづけを行い、予測精度のさらなる向上に貢献しない不要なモデルは無視することが出来ます。

そのためにはstacksの為に学習済みモデルを保存しておく必要があります。
tune_gridやfit_resampleによってモデルを学習させる際にcontrol引数に

```{r}
ctrl_grid <- control_stack_grid()
ctrl_res <- control_stack_resamples()
```

を入力しましょう。
上記は下記と同じ効果を持っています。

```{r}
control_grid(
  save_pred = T,
  save_workflow = T)

control_resamples(
  save_pred = T,
  save_workflow = T)
```


チューニング結果、もしくは交差検証を行ったfit_resamplesの結果をスタッキングします。
折角なのでglmnetを使った交差検証の結果も使いスタッキングに使いましょう。

```{r}
glmnet_recipe <- 
  recipe(formula = Sale_Price ~ ., data = train) %>% 
  step_novel(all_nominal(), -all_outcomes()) %>% 
  step_dummy(all_nominal(), -all_outcomes()) %>% 
  step_zv(all_predictors()) %>% 
  step_normalize(all_predictors(), -all_nominal()) 

glmnet_spec <- 
  linear_reg(penalty = 0.001) %>% 
  set_mode("regression") %>% 
  set_engine("glmnet") 

glmnet_workflow <- 
  workflow() %>% 
  add_recipe(glmnet_recipe) %>% 
  add_model(glmnet_spec) 

fit_res <- fit_resamples(glmnet_workflow,splits,
              control = control_resamples(save_pred = T,save_workflow = T))
```

スタッキングの手順としては、stacks()によって箱を用意してadd_candidatesで学習済みの結果を渡すだけです。

```{r}
library(stacks)

model_stacks <- 
  stacks() %>%
  add_candidates(xgboost_tune) %>%
  add_candidates(fit_res)

```

その後、blend_predictions()によってLASSOを当てはめます。

```{r}
model_stacks_blended <-
  model_stacks %>%
  blend_predictions()
```

スタッキングしたモデルを使いすべての訓練データを使って学習します。

```{r}
model_stacks_blende_trained <-
  model_stacks_blended %>%
  fit_members()
```

最後にtestデータを使って予測です。

```{r}
member_preds <- 
  test %>%
  select(Sale_Price) %>%
  bind_cols(predict(model_stacks_blende_trained, 
                    test, 
                    members = TRUE))

```


これでスタッキングしたモデルによる予測が出来ました。

# チューニング時間の短縮




上記のチューニングの場合、
- 補填するパラメタパターンの組み合わせが5パターン
- vfoldとして5個
合計25モデルが作成されることになります。

そして実際の学習にかかった時間は**518秒**でした。

もっと時間を短縮しなければ、最良のモデルを求める作業を限られた時間で繰り返し実行することが難しくなります。

時間短縮の方法の一つは「マルチコアでの実行」です。
これは[Optimizations and Parallel Processing](https://tune.tidymodels.org/articles/extras/optimizations.html)あたりをご参照ください。
もう一つの方法として[finetune](https://finetune.tidymodels.org/index.html)パッケージを紹介します。


tidymodelsのfinetuneパッケージを簡単に紹介すると「性能の上がらないパラメタの組み合わせを中止させる」ということになります。

実際に見ていただいた方が早いでしょう。

```{r}
library(finetune)
rece_res <- 
tune_race_anova(xgboost_workflow, 
                resamples = splits, 
                grid = 5, 
                control = control_race(save_workflow = T,save_pred = T))

rece_res %>% 
  plot_race()

```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/95760d12-e023-d361-1633-f4abbaeb4b66.png)


実行時間は**424秒**でした。

結果を見てみると

```{r}
rece_res
```


```
# Tuning results
# 5-fold cross-validation using stratification 
# A tibble: 5 x 6
  splits             id    .order .metrics           .notes           .predictions         
  <list>             <chr>  <int> <list>             <list>           <list>               
1 <split [1874/468]> Fold3      2 <tibble [10 x 10]> <tibble [0 x 1]> <tibble [2,340 x 10]>
2 <split [1875/467]> Fold4      1 <tibble [10 x 10]> <tibble [0 x 1]> <tibble [2,335 x 10]>
3 <split [1875/467]> Fold5      3 <tibble [10 x 10]> <tibble [0 x 1]> <tibble [2,335 x 10]>
4 <split [1871/471]> Fold1      5 <tibble [2 x 10]>  <tibble [0 x 1]> <tibble [471 x 10]>  
5 <split [1873/469]> Fold2      4 <tibble [2 x 10]>  <tibble [0 x 1]> <tibble [469 x 10]>  
```



.metricsが2つしか入っていないレコードがあります。
性能が上がらないモデルは検証されることなく中止されるため、すべてを試していないのです。
更にvfoldを増やした場合に時間的なメリットを感じることができるでしょう。

# 時間の短縮がモデルスタッキングに与える悪い影響

このfinetuneオブジェクトをstacksの関数に詰めることができます。

なんとなく予想がつくかと思いますが、先に結果から申し上げますと
「**途中までしか検証されていないモデルも含まれる**」ため、スタッキング後のモデルはよく確認しなければいけません。

```{r}
race_stack <-
  stacks() %>% 
  add_candidates(rece_res) %>% 
  add_candidates(fit_res)

race_stack <-
  race_stack%>% 
  blend_predictions() 

race_stack_trained<-
race_stack %>% 
  fit_members()
```


これが問題なく通ります。

バグではなく、公式にも認識されている仕様で今後修正される予定は今のところ出ていないようです。
（※2021/10　のチェック。見逃していたら教えてください。）



# finetuneについての余談

finetuneは時間短縮に使うことが可能です。
性能の悪い領域を見捨てることで最も良いパラメタを素早く見つけることが可能です。
つまり、広い範囲のチューニングのためのパラメータグリッドを用意しておいて、何度もfinetuneを使い最適なパラメタを求めることが短時間でできるようになります。

グリッドを複数回繰り返し実行することも可能ですが、さらに便利な関数もありベイズ最適化によって逐次探索する関数も実装されています。
ベイズ最適化や焼きなまし法は、上手に言語化できるレベルに達していないので以下をご参照ください。

[Optimization of model parameters via simulated annealing](https://finetune.tidymodels.org/reference/tune_sim_anneal.html)

# 最後に

finetuneとstacksの性質を理解し、正しく判断して利用しましょう。

enjoy!!