Tidymodels による機械学習モデル構築

俺様コードによる tidymodels を用いた機械学習モデル構築のサンプルを記載する。
※基礎的な集計の段階は終了しているという前提で話を進めていく

使用したパッケージの一覧

• rsample
• recipes
• parsnip
• dials
• yardstick

データの分割(訓練/テスト)

データを訓練用とテスト用に分割する。
strata 引数に目的変数を指定する事で目的変数の分布を変えずにデータを分割する事が可能。

• rsample::initial_split
• rsample::analysis / rsample::assesment

library(tidyverse)
library(tidymodels)

# 元データを訓練用とテスト用に分割
# strata を指定する事でクラスの分布を保持したまま分割
lst.splitted <- rsample::initial_split(iris, prop = 0.75, strata = "Species") %>% {
  list(
    train = rsample::analysis(.),
    test  = rsample::assessment(.)
  )
}

クロスバリデーション用データの分割

クロスバリデーション用にデータを分割する。
strata に目的変数を指定しているのは前述した rsample::initial_split と同じ理由。

• rsample::vfold_cv

# クロスバリデーション用の分割を定義
df.cv <- rsample::vfold_cv(
  lst.splitted$train,
  v = 5,
  strata = "Species"
)

前処理の定義

事前に実施した基礎集計の結果やドメイン知識を元にデータの前処理を行う。
いわゆる特徴量エンジニアリングを行う箇所。

※ここで定義している処理はサンプルとして適当に行っているものであり、実施によって精度は下がっていると思われるw

• recipes::recipe
• recipes::step_rm
• recipes::step_log
• recipes::step_center
• recipes::step_scale

# 前処理レシピの定義
recipe <- recipes::recipe(Species ~ ., lst.splitted$train) %>%

  # Sepal.Width を削除
  recipes::step_rm(Sepal.Width) %>%

  # Sepal.Length を対数変換
  recipes::step_log(Sepal.Length) %>%

  # 説明変数を基準化
  recipes::step_center(all_predictors()) %>%
  recipes::step_scale(all_predictors())

モデルの定義

今回は RandomForest でモデルを定義する。
各ハイパーパラメータに parsnip::varying を指定する事で後から複数の値を設定可能。今回はグリッドサーチでパラメータを指定する(後述)。

parsnip::set_engine では指定した engine パッケージ特有のパラメータをセットする事ができる。
下記の例では ranger::ranger 関数の num.threads パラメータ(並列処理に用いるコア数)を指定している。

• parsnip::rand_forest
• parsnip::varying
• parsnip::set_engine

# モデルの定義
model <- parsnip::rand_forest(
  mode = "classification",
  mtry = parsnip::varying(),
  min_n = parsnip::varying(),
  trees = parsnip::varying()
) %>%
  parsnip::set_engine(engine = "ranger", num.threads = 4)

グリッドサーチの定義

3 x 3 x 3 = 27 パターンのハイパーパラメータの組み合わせを定義。
今回はグリッドサーチを行っているが dials::grid_random を用いてランダムサーチを行う事も可能。

• dials::grid_regular

# グリッドサーチ用の組み合わせパターンを定義
df.grid.params <- dials::grid_regular(
  dials::mtry  %>% dials::range_set(c(1, 3)),
  dials::min_n %>% dials::range_set(c(2, 6)),
  dials::trees %>% dials::range_set(c(500, 1500)),
  levels = 3
)

モデルの学習と評価

処理の構造としては下記の 2 重ループとなっている。

1. ハイパーパラメータ一覧のループ
2. クロスバリデーションによる分割のループ

上記 2 重ループの内部で下記の各ステップを実施。

1. モデルの学習
2. 学習済モデルによる予測
3. モデルの評価

上記ループが終了したら 1 段目のループ(ハイパラ一覧のループ)の内部にて下記を実施。

• CV 毎の評価スコアを平均してハイパーパラメータ毎の評価スコアを算出

上記によって得られる結果を評価スコアでソートする事で最も性能の良いハイパーパラメータを特定する。

• recipes::prep
• recipes::bake
• parsnip::fit
• yardstick::metric_set
• yardstick::accuracy
• yardstick::precision
• yardstick::recall
• yardstick::f_meas

df.results <-

  # ハイパーパラメータをモデルに適用
  # ※2020.05 時点で下記のコードは動かない
  merge(df.grid.params, model) %>%
  # 下記のコードに改修(2020.05.19)
  purrr::pmap(df.grid.params, set_args, object = model) %>%

  # ハイパーパラメータの組み合わせごとにループ
  purrr::map(function(model.applied) {

    # クロスバリデーションの分割ごとにループ
    purrr::map(df.cv$splits, model = model.applied, function(df.split, model) {

      # 前処理済データの作成
      df.train <- recipe %>%
        recipes::prep() %>%
        recipes::bake(rsample::analysis(df.split))
      df.test <- recipe %>%
        recipes::prep() %>%
        recipes::bake(rsample::assessment(df.split))

      model %>%

        # モデルの学習
        {
          model <- (.)

          parsnip::fit(model, Species ~ ., df.train)
        } %>%

        # 学習済モデルによる予測
        {
          fit <- (.)

          list(
            train = predict(fit, df.train, type = "class")[[1]],
            test  = predict(fit, df.test,  type = "class")[[1]]
          )
        } %>%

        # 評価
        {
          lst.predicted <- (.)

          # 評価指標の一覧を定義
          metrics <- yardstick::metric_set(
            yardstick::accuracy,
            yardstick::precision,
            yardstick::recall,
            yardstick::f_meas
          )

          # train データでモデルを評価
          df.result.train <- df.train %>%
            dplyr::mutate(
              predicted = lst.predicted$train
            ) %>%
            metrics(truth = Species, estimate = predicted) %>%
            dplyr::select(-.estimator) %>%
            dplyr::mutate(
              .metric = stringr::str_c("train", .metric, sep = "_")
            ) %>%
            tidyr::spread(key = .metric, value = .estimate)

          # test データでモデルを評価
          df.result.test <- df.test %>%
            dplyr::mutate(
              predicted = lst.predicted$test
            ) %>%
            metrics(truth = Species, estimate = predicted) %>%
            dplyr::select(-.estimator) %>%
            dplyr::mutate(
              .metric = stringr::str_c("test", .metric, sep = "_")
            ) %>%
            tidyr::spread(key = .metric, value = .estimate)

          dplyr::bind_cols(
            df.result.train,
            df.result.test
          )
        }
    }) %>%

      # CV 分割全体の平均値を評価スコアとする
      purrr::reduce(dplyr::bind_rows) %>%
      dplyr::summarise_all(mean)
  }) %>%

  # 評価結果とパラメータを結合
  purrr::reduce(dplyr::bind_rows) %>%
  dplyr::bind_cols(df.grid.params) %>%

  # 評価スコアの順にソート(昇順)
  dplyr::arrange(
    desc(test_accuracy)
  ) %>%

  dplyr::select(
    mtry,
    min_n,
    trees,

    train_accuracy,
    train_precision,
    train_recall,
    train_f_meas,

    test_accuracy,
    test_precision,
    test_recall,
    test_f_meas
  )

上記で得られる df.results のサンプル。

mtry	min_n	trees	train_accuracy	train_precision	train_recall	train_f_meas	test_accuracy	test_precision	test_recall	test_f_meas
2	2	500	1.000	1.000	1.000	1.000	0.946	0.949	0.946	0.946
2	4	500	0.991	0.991	0.991	0.991	0.946	0.949	0.946	0.946
3	4	500	0.987	0.987	0.987	0.987	0.946	0.949	0.946	0.946
2	6	500	0.985	0.985	0.985	0.985	0.946	0.949	0.946	0.946
3	6	500	0.982	0.983	0.982	0.982	0.946	0.949	0.946	0.946

ベストモデルの構築

最も性能の良い(=評価スコア最大)モデルを構築する。
この段階では訓練データの全体を用いて学習を行う事に注意。

# 訓練データに前処理を適用
df.train.baked <- recipe %>%
  recipes::prep() %>%
  recipes::bake(lst.splitted$train)

# 最も性能の良いハイパーパラメータを用いたモデルを構築
best_model <- update(
  model,
  mtry  = df.results[1,]$mtry,
  min_n = df.results[1,]$min_n,
  trees = df.results[1,]$trees
) %>%

  # 訓練データ全体を用いてモデルの学習を行う
  parsnip::fit(Species ~ ., df.train.baked)

テストデータによる検証

テストデータを用いた評価を行い、構築したモデルの汎化性能を検証する。

# テストデータに前処理を適用
df.test.baked <- recipe %>%
  recipes::prep() %>%
  recipes::bake(lst.splitted$test)  

# 汎化性能を検証
df.test.baked %>%

  # ベストモデルを用いて予測
  dplyr::mutate(
    predicted = predict(best_model, df.test.baked)[[1]]
  ) %>%

  # 精度(Accuracy)を算出
  yardstick::accuracy(Species, predicted)

まとめ

必要になる度に前に書いたコードを思い出したりコピペしたりでノウハウをまとめられていなかったので、今回は良い機会になった。
機械学習周りは scikit-learn 一択かなと思ってた頃もあったけど tidymodels がいい感じなのでぜひ使う人が増えてくれると嬉し。
他の人がどんな感じでやってるのかも知りたいところ。

時系列のクロスバリデーション

tidymodels に含まれる rsample パッケージの rolling_origin という関数が便利そうだったので試してみた内容を記載してみる。表題の通り時系列データに対するクロスバリデーションを行う。

この関数を知る切っ掛けになったのがこちら。ありがとうございます。

https://blog.hoxo-m.com/entry/2019/06/08/220307

rolling_origin 関数

時系列データのクロスバリデーションに関する考え方はこちらのサイト Rでのナウなデータ分割のやり方: rsampleパッケージによる交差検証 や書籍 前処理大全 に任せるとして、rolling_origin 関数の skip 引数の意味が分かりづらかったので少し試してみる事にする。

まず、今回使用するデータを準備。
"S4248SM144NCEN" というカラム名が何かは知らないけどデータの名称からして恐らく飲酒量とかその辺りだと思われるのでとりあえず "amount" という項目名に変更しておく。

library(tidyverse)
library(tidymodels)

df.drinks <- rsample::drinks %>%
  dplyr::rename(amount = S4248SM144NCEN)

#	date	amount
1	1992-01-01	3,459
2	1992-02-01	3,458
3	1992-03-01	4,002

このデータは 1992/01/01 から 2017/09/01 まで続く月次データであり、全部で 309 件のレコードが存在する。
このデータを下記の要件にて分割する事を考えてみる。

• 学習期間: 24ヶ月
• 検証期間: 12ヶ月
• スライド: 1ヶ月

上記を実現したのがこちら

df.rolling <- rsample::rolling_origin(
  df.drinks,
  initial = 24,
  assess = 12,
  skip = 0,
  cumulative = F
)

df.rolling
# Rolling origin forecast resampling 
# A tibble: 274 x 2
   splits          id      
                
 1 <split [24/12]> Slice001
 2 <split [24/12]> Slice002
 3 <split [24/12]> Slice003
 4 <split [24/12]> Slice004
 5 <split [24/12]> Slice005
 6 <split [24/12]> Slice006
 7 <split [24/12]> Slice007
 8 <split [24/12]> Slice008
 9 <split [24/12]> Slice009
10 <split [24/12]> Slice010
# … with 264 more rows

学習データの 1 番目を取り出してみた。確かに 1992/1/1 を含む 24 ヶ月間のデータとなっている。

df.rolling$splits[[1]] %>%
  rsample::analysis()

# A tibble: 24 x 2
   date       amount
         
 1 1992-01-01   3459
 2 1992-02-01   3458
 3 1992-03-01   4002
 4 1992-04-01   4564
 5 1992-05-01   4221
 6 1992-06-01   4529
 7 1992-07-01   4466
 8 1992-08-01   4137
 9 1992-09-01   4126
10 1992-10-01   4259
# … with 14 more rows

検証用データの 1 番目を取り出してみる。こちらも想定通り 1994 年の 12 ヶ月分のデータとなっている。

df.rolling$splits[[1]] %>%
  rsample::assessment()

# A tibble: 12 x 2
   date       amount
         
 1 1994-01-01   3075
 2 1994-02-01   3377
 3 1994-03-01   4443
 4 1994-04-01   4261
 5 1994-05-01   4460
 6 1994-06-01   4985
 7 1994-07-01   4324
 8 1994-08-01   4719
 9 1994-09-01   4374
10 1994-10-01   4248
11 1994-11-01   4784
12 1994-12-01   4971

2 番目の学習用/検証用のデータも取り出してみる。想定通り 1 ヶ月スライドしたデータとなっている。

df.rolling$splits[[2]] %>%
  rsample::analysis()
# A tibble: 24 x 2
   date       amount
         
 1 1992-02-01   3458
 2 1992-03-01   4002
 3 1992-04-01   4564
 4 1992-05-01   4221
 5 1992-06-01   4529
 6 1992-07-01   4466
 7 1992-08-01   4137
 8 1992-09-01   4126
 9 1992-10-01   4259
10 1992-11-01   4240
# … with 14 more rows

df.rolling$splits[[2]] %>%
  rsample::assessment()
# A tibble: 12 x 2
   date       amount
         
 1 1994-02-01   3377
 2 1994-03-01   4443
 3 1994-04-01   4261
 4 1994-05-01   4460
 5 1994-06-01   4985
 6 1994-07-01   4324
 7 1994-08-01   4719
 8 1994-09-01   4374
 9 1994-10-01   4248
10 1994-11-01   4784
11 1994-12-01   4971
12 1995-01-01   3370

ここで最初に疑問に思ったのが rolling_origin 関数の skip 引数の仕様である。例えば下記のようなデータ仕様を考えてみる。

• 学習期間: 24ヶ月
• 検証期間: 12ヶ月
• スライド: 12ヶ月

つまり 1 年おきにデータを学習/予測(ex. 年の始めに 2 年分のデータを使って直近 1 年間の予測を行う)するようなケースとなる。 この時に直感的には下記のようなコードを想定していた。

df.rolling2 <- rsample::rolling_origin(
  df.drinks,
  initial = 24,
  assess = 12,
  skip = 12,
  cumulative = F
)

12 ヶ月データをスキップするので skip 引数に 12 を指定しているが、これだと下記のように意図した通りにはならない。 2 番目の学習データが 1993/1/1 ではなく 1993/2/1 からのスタートになってしまっており、実際には 13 ヶ月のスキップとなってしまっている。

df.rolling2$splits[[2]] %>%
  rsample::analysis()
# A tibble: 24 x 2
   date       amount
         
 1 1993-02-01   3261
 2 1993-03-01   4160
 3 1993-04-01   4377
 4 1993-05-01   4307
 5 1993-06-01   4696
 6 1993-07-01   4458
 7 1993-08-01   4457
 8 1993-09-01   4364
 9 1993-10-01   4236
10 1993-11-01   4500
# … with 14 more rows

これはそもそも rolling_origin 関数における skip 引数のデフォルト値が 0 であり、この指定によって 1 ヶ月おきのスキップを意図している事による。 つまりそこから更にずらす幅を追加したい時に skip 引数を使いなさいという事なのだろうなと。

これはドキュメントをちゃんと読むとそれっぽい事が記載されていた。ちゃんと読みなさいよと

When skip = 0, the resampling data sets will increment by one position.

あとから考えてみると skip = 0 の指定でスライドしないのであれば無限ループになってしまう orz

なので、上記の仕様であれば下記のコードが正解となる。

df.rolling3 <- rsample::rolling_origin(
  df.drinks,
  initial = 24,
  assess = 12,
  skip = 11, # 11 = 12 - 1
  cumulative = F
)

df.rolling3$splits[[2]] %>%
  rsample::analysis()
# A tibble: 24 x 2
   date       amount
         
 1 1993-01-01   3031
 2 1993-02-01   3261
 3 1993-03-01   4160
 4 1993-04-01   4377
 5 1993-05-01   4307
 6 1993-06-01   4696
 7 1993-07-01   4458
 8 1993-08-01   4457
 9 1993-09-01   4364
10 1993-10-01   4236
# … with 14 more rows

最初から複数月のスライドを意図したデータを作ろうとしていたのが失敗で、スライド幅 1 のケースから順番に試していれば混乱しなかったんだろうなと。難しい

クロスバリデーションによるモデル比較

せっかく rolling_origin の使い方がだいたい分かったので、これを用いて時系列モデルの比較を行ってみる事にする。

まず比較対象となるモデルの一覧を list で作成する。実際には当該リストに含まれるのは学習済みモデルを返す関数である事に注意。

lst.models <- list(
  # ローカルレベルモデル + 確定的季節成分
  d1_ss = function(df.data) {
    library(KFAS)
    model <- SSModel(
      amount ~
        SSMtrend(degree = 1, Q = NA) +
        SSMseasonal(period = 12, Q = 0),
      data = df.data,
      H = NA
    )
    fitSSM(
      model,
      inits = c(1, 1),
      method = "BFGS"
    )
  },

  # ローカルレベルモデル + 確率的季節成分
  d1_sv = function(df.data) {
    library(KFAS)
    model <- SSModel(
      amount ~
        SSMtrend(degree = 1, Q = NA) +
        SSMseasonal(period = 12, Q = NA),
      data = df.data,
      H = NA
    )
    fitSSM(
      model,
      inits = c(1, 1, 1),
      method = "BFGS"
    )
  },

  # ローカル線形トレンドモデル + 確定的季節成分
  d2_ss = function(df.data) {
    library(KFAS)
    model <- SSModel(
      amount ~
        SSMtrend(degree = 2, Q = list(matrix(NA), matrix(NA))) +
        SSMseasonal(period = 12, Q = 0),
      data = df.data,
      H = NA
    )
    fitSSM(
      model,
      inits = c(1, 1, 1),
      method = "BFGS"
    )
  },

  # ローカル線形トレンドモデル + 確率的季節成分
  d2_sv = function(df.data) {
    library(KFAS)
    model <- SSModel(
      amount ~
        SSMtrend(degree = 2, Q = list(matrix(NA), matrix(NA))) +
        SSMseasonal(period = 12, Q = NA),
      data = df.data,
      H = NA
    )
    fitSSM(
      model,
      inits = c(1, 1, 1, 1),
      method = "BFGS"
    )
  }
)

リストの要素を指定してデータを渡せば KFAS による学習済みモデルが返却される。

lst.models$d1_ss(df.drinks)

### ↓以下結果 ###

$optim.out
$optim.out$par
[1]  9.355875 12.207159

$optim.out$value
[1] 2290.374

$optim.out$counts
function gradient 
      29       10 

$optim.out$convergence
[1] 0

$optim.out$message
NULL


$model
Call:
SSModel(formula = amount ~ SSMtrend(degree = 1, Q = NA) + SSMseasonal(period = 12, 
    Q = 0), data = df.data, H = NA)

State space model object of class SSModel

Dimensions:
[1] Number of time points: 309
[1] Number of time series: 1
[1] Number of disturbances: 2
[1] Number of states: 12
Names of the states:
 [1]  level        sea_dummy1   sea_dummy2   sea_dummy3   sea_dummy4   sea_dummy5   sea_dummy6   sea_dummy7   sea_dummy8 
[10]  sea_dummy9   sea_dummy10  sea_dummy11
Distributions of the time series:
[1]  gaussian

Object is a valid object of class SSModel.

こんな感じに予測も可能。

predict(lst.models$d1_ss(df.drinks)$model)

Time Series:
Start = 1 
End = 309 
Frequency = 1 
             fit
  [1,]  2719.484
  [2,]  3162.725
  [3,]  4213.908
  [4,]  4128.718
  [5,]  4758.920 ...

上記を用いて下記のようにモデル毎に評価結果を返す関数を作って実験を行う。
評価指標は CV データ毎に算出した RMSE の平均を用いる事にする。

# h: 予測する範囲 ex. h=1 で 1 期先予測
# skip: スキップ数
rmse_per_model_cv <- function(h, skip = 0) {
  # 時系列 CV の作成
  df.rolling <- rsample::rolling_origin(df.drinks, initial = 24, assess = h,  cumulative = F, skip = skip)

  # 各モデルの適用
  purrr::map(lst.models, function(model) {

    # クロスバリデーション
    purrr::map_dbl(df.rolling$splits, function(split) {
      df.train <- rsample::analysis(split)
      df.test  <- rsample::assessment(split)

      # モデルの学習
      fit <- model(df.train)

      # 学習済みモデルによる予測と RMSE の算出
      # columns: date, amount, predicted
      df.test %>%
        dplyr::mutate(
          predicted = predict(fit$model, n.ahead = h)[, "fit"] %>% as.numeric()
        ) %>%
        yardstick::rmse(amount, predicted) %>%
        .$.estimate
    })
  })
}

モデル毎に 24 ヶ月の学習データを用いて 12 期先までの予測を行う試行を 23 回行い、23 回分の RMSE の平均値をモデルの評価値とする。

# 予測結果の一覧を取得
rmses <- rmse_per_model_cv(h = 12, skip = 11)

# モデル毎に RMSE の平均値を算出
purrr::map_dbl(rmses, ~ mean(.))

   d1_ss    d1_sv    d2_ss    d2_sv 
561.9146 472.5309 442.3523 449.9628

今回のデータだとモデル d2_ss(ローカル線形トレンド+確定的季節成分) の予測精度が最も高いという結果になった。

Iris で k-NN Feature Extraction

前回 は k-NN Feature Extraction を自分で作成したテストデータに適用して試してみた。 そこで今回は Iris データに対して適用してみたいと思う。

関数のコードは https://github.com/you1025/knn_feature_extraction/blob/master/knn_feature_extraction.R

Iris

library(tidyverse)

df.iris <- iris %>%
  tibble::as_tibble()

head(df.iris)

#	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
1	5.1	3.5	1.4	0.2	setosa
2	4.9	3.0	1.4	0.2	setosa
3	4.7	3.2	1.3	0.2	setosa
4	4.6	3.1	1.5	0.2	setosa
5	5.0	3.6	1.4	0.2	setosa
6	5.4	3.9	1.7	0.4	setosa

散布図で可視化

df.iris %>%
  ggplot(aes(Sepal.Length, Sepal.Width)) +
    geom_point(aes(colour = Species), alpha = 2/3)

k-NN Feature Extraction の適用

df.knn_d_columns <- df.iris %>%
  add_knn_d_columns(Species)

head(df.knn_d_columns)

#	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species	setosa_1	versicolor_1	virginica_1
1	5.1	3.5	1.4	0.2	setosa	0.1000000	2.0904545	3.5916570
2	4.9	3.0	1.4	0.2	setosa	0.1414214	1.9131126	3.4799425
3	4.7	3.2	1.3	0.2	setosa	0.1414214	2.0856654	3.6083237
4	4.6	3.1	1.5	0.2	setosa	0.1414214	1.9157244	3.4205263
5	5.0	3.6	1.4	0.2	setosa	0.1414214	2.1424285	3.6166283
6	5.4	3.9	1.7	0.4	setosa	0.3316625	2.0566964	3.4263683

2 軸を変更しながら可視化

# setosa & versicolor
g1 <- df.knn_d_columns %>%
  ggplot(aes(setosa_1, versicolor_1)) +
    geom_point(aes(colour = Species), alpha = 2/3)

# versicolor & virginica
g2 <- df.knn_d_columns %>%
  ggplot(aes(versicolor_1, virginica_1)) +
    geom_point(aes(colour = Species), alpha = 2/3)

# virginica & setosa
g3 <- df.knn_d_columns %>%
  ggplot(aes(virginica_1, setosa_1)) +
    geom_point(aes(colour = Species), alpha = 2/3)

gridExtra::grid.arrange(g1, g2, g3)

どの組み合わせもいい感じに分離できそう。

3 次元プロットも試してみた

plotly::plot_ly(
  x = df.knn_d_columns$setosa_1,
  y = df.knn_d_columns$versicolor_1,
  z = df.knn_d_columns$virginica_1,
  type = "scatter3d",
  mode = "markers",
  color = df.knn_d_columns$Species,
  size = 0
)

いいね