TidyModels: modeltime + cross-validation
UPDATE: 2022-12-24 04:35:24
はじめに
tidymodelsパッケージの使い方をいくつかのノートに分けてまとめている。tidymodelsパッケージは、統計モデルや機械学習モデルを構築するために必要なパッケージをコレクションしているパッケージで、非常に色んなパッケージがある。ここでは、今回はmodeltimeパッケージ
+
cross-validationについてまとめていく。モデルの数理的な側面や機械学習の用語などは、このノートでは扱わない。
下記の公式ドキュメントやtidymodelsパッケージに関する書籍を参考にしている。
modeltimeパッケージの目的
modeltimeパッケージは、時系列モデルの構築を効率よく行うためのパッケージで、tidymodelsパッケージと組み合わせて時系列モデルを作成する際に使用される。時系列予測のための体系的なワークフローを提供し、時系列ブーストモデルなど新しいアルゴリズムもカバーしている点が特徴的。
前回は基本的な使い方をまとめており、時系列モデルのクロスバリデーションを扱う。
上記のサイトをなぞりながら、不明点を追記しながらまとめていく。
モデルの作成
まずは必要なパッケージを読み込んでおく。
library(tidymodels)
library(modeltime)
library(tidyverse)
library(lubridate)
library(timetk)
# クロスバリデーションには下記のパッケージが追加で必要
library(modeltime.resample)前回同様、ここでもtimetkパッケージ内のm4_monthlyデータセットのidがM750の部分を使用する。データは、1990-01-01から始まり、2015-06-01まである。レコード数は306レコードで、間隔は月単位のデータになっている。
m750 <- m4_monthly %>% filter(id == "M750")
m750 %>%
summarise(start_date = min(date), end_date = max(date), records = n())## # A tibble: 1 × 3
## start_date end_date records
## <date> <date> <int>
## 1 1990-01-01 2015-06-01 306このようなデータをここでは扱う。
m750 %>%
plot_time_series(.data = ., .date_var = date, .value = value, .interactive = FALSE)
クロスバリデーション用のフォールドはtime_series_cv関数で作成できる。今回は回の設定でクロスバリデーションのフォールドを作成する。
assess: 評価ウィンドウは2年initial: 学習ウィンドウは5年skip: リサンプルセット間のシフトは2年slice_limit: 生成するリサンプル数は3コ
resamples_tscv <- time_series_cv(
data = m750,
assess = "2 years",
initial = "5 years",
skip = "1 years",
slice_limit = 3
)
resamples_tscv %>%
mutate(
analysis = map(.x = splits, .f = function(x){analysis(x)}),
assessment = map(.x = splits, .f = function(x){assessment(x)}),
analysis_min_date = map(.x = analysis, .f = function(x){min(x$date)}),
analysis_max_date = map(.x = analysis, .f = function(x){max(x$date)}),
assessment_min_date = map(.x = assessment, .f = function(x){min(x$date)}),
assessment_max_date = map(.x = assessment, .f = function(x){max(x$date)}),
) %>%
unnest(analysis_min_date,
analysis_max_date,
assessment_min_date,
assessment_max_date
)## # A tibble: 3 × 8
## splits id analysis assessment analysis_mi…¹ analysis…² assessme…³
## <list> <chr> <list> <list> <date> <date> <date>
## 1 <split [60/24]> Slice1 <tibble> <tibble> 2008-07-01 2013-06-01 2013-07-01
## 2 <split [60/24]> Slice2 <tibble> <tibble> 2007-07-01 2012-06-01 2012-07-01
## 3 <split [60/24]> Slice3 <tibble> <tibble> 2006-07-01 2011-06-01 2011-07-01
## # … with 1 more variable: assessment_max_date <date>, and abbreviated variable
## # names ¹analysis_min_date, ²analysis_max_date, ³assessment_min_dateそのままplot_time_series_cv_plan関数に渡すことで、クロスバリデーションの状態を可視化できる。
resamples_tscv %>%
tk_time_series_cv_plan() %>%
plot_time_series_cv_plan(date, value, .interactive = FALSE)
データの準備が整ったので、モデルを設定する。前回同様のモデルを利用するので、実行結果を乗せておく。
# Model1 ARIMA
model_fit_arima <- arima_reg() %>%
set_engine(engine = "auto_arima") %>%
fit(value ~ date, data = m750)
# Model2 Prophet
model_fit_prophet <- prophet_reg(seasonality_yearly = TRUE) %>%
set_engine(engine = "prophet") %>%
fit(value ~ date, data = m750)
# Recipe For RandomForest & Prophet Boost
recipe <- recipe(value ~ date, m750) %>%
step_timeseries_signature(date) %>%
step_rm(
contains("am.pm"),
contains("hour"),
contains("minute"),
contains("second"),
contains("xts")
) %>%
step_fourier(date, period = 365, K = 5) %>%
step_dummy(all_nominal())
# Model Spec For RandomForest
model_randomforest <- rand_forest(trees = 1000, min_n = 50, mtry = 15) %>%
set_engine("ranger") %>%
set_mode("regression")
# Model Fitting For RandomForest
model_fit_randomforest <- workflow() %>%
add_model(model_randomforest) %>%
add_recipe(recipe %>% step_rm(date)) %>%
fit(m750)
# Model Spec For Prophet Boost
model_prophet_boost <- prophet_boost(seasonality_yearly = TRUE) %>%
set_engine("prophet_xgboost")
# Model Fitting For Prophet Boost
model_fit_prophet_boost <- workflow() %>%
add_model(model_prophet_boost) %>%
add_recipe(recipe) %>%
fit(m750)モデルのフィッティングが完了したらモデルテーブルに追加する(modeltime_table関数)。これは、クロスバリデーションするのに、訓練データでフィッティングする理由は、モデルテーブルには、フィッティングしてあるモデルしか登録できないため(All objects must be fitted workflow or parsnip models)。この後にクロスバリデーションを行なうことになるので形式的な作業である。
models_tbl <- modeltime_table(
model_fit_arima,
model_fit_prophet,
model_fit_prophet_boost,
model_fit_randomforest
)
models_tbl## # Modeltime Table
## # A tibble: 4 × 3
## .model_id .model .model_desc
## <int> <list> <chr>
## 1 1 <fit[+]> ARIMA(0,1,1)(0,1,1)[12]
## 2 2 <fit[+]> PROPHET
## 3 3 <workflow> PROPHET W/ XGBOOST ERRORS
## 4 4 <workflow> RANGERリサンプル予測の生成はmodeltime_fit_resamples関数で行なう。この時、内部的には、各モデルはリサンプルの各トレーニングセットで再学習される。モデルテーブルにリサンプル予測(.resample_results)の列が追加される。
resamples_fitted <- models_tbl %>%
modeltime_fit_resamples(
resamples = resamples_tscv,
control = control_resamples(verbose = FALSE)
)
resamples_fitted## # Modeltime Table
## # A tibble: 4 × 4
## .model_id .model .model_desc .resample_results
## <int> <list> <chr> <list>
## 1 1 <fit[+]> ARIMA(0,1,1)(0,1,1)[12] <rsmp[+]>
## 2 2 <fit[+]> PROPHET <rsmp[+]>
## 3 3 <workflow> PROPHET W/ XGBOOST ERRORS <rsmp[+]>
## 4 4 <workflow> RANGER <rsmp[+]>返されるデータフレームは入れ子の入れ子構造になっている。.resample_resultsの中を除くと、データフレームが返されるものの、予測値はまだ見えないので、もう1階層潜る必要がある。
resamples_fitted %>%
# ARIMA
pluck(".resample_results", 1)## # Resampling results
## # NA
## # A tibble: 3 × 5
## splits id .metrics .notes .predictions
## <list> <chr> <list> <list> <list>
## 1 <split [60/24]> Slice1 <tibble [1 × 4]> <tibble [0 × 3]> <tibble [24 × 4]>
## 2 <split [60/24]> Slice2 <tibble [1 × 4]> <tibble [0 × 3]> <tibble [24 × 4]>
## 3 <split [60/24]> Slice3 <tibble [1 × 4]> <tibble [0 × 3]> <tibble [24 × 4]>ということで、もう1階層潜って予測値を確認する。.metricsには評価指標が記録されている。
resamples_fitted %>%
# ARIMA
pluck(".resample_results", 1) %>%
# 1つ目のフォールドの予測値
pluck(".predictions", 1)## # A tibble: 24 × 4
## .pred .row value .config
## <dbl> <int> <dbl> <chr>
## 1 8977. 283 9030 Preprocessor1_Model1
## 2 9597. 284 9620 Preprocessor1_Model1
## 3 9907. 285 10050 Preprocessor1_Model1
## 4 10737. 286 10610 Preprocessor1_Model1
## 5 10887. 287 10760 Preprocessor1_Model1
## 6 10747. 288 10570 Preprocessor1_Model1
## 7 10807. 289 10730 Preprocessor1_Model1
## 8 10797. 290 10660 Preprocessor1_Model1
## 9 10977. 291 10900 Preprocessor1_Model1
## 10 10957. 292 11020 Preprocessor1_Model1
## # … with 14 more rowsplot_modeltime_resamples関数を使用することで、モデルのクロスバリデーションの精度を可視化して確認できる。評価指標ごとにファセットが区切られており、Y軸に各クロスバリデーションのフォールド、横軸が指標の精度である。
例えば、モデル4のランダムフォレストは薄い金色の点を見ればよく、どの指標でも精度がよろしくない。
resamples_fitted %>%
plot_modeltime_resamples(
.point_size = 2,
.point_alpha = 0.5,
.interactive = FALSE
)
グラフではなく、数表が欲しい場合は、modeltime_resample_accuracy関数を利用する。これを見る限り、Prophet
Boostが良い事がわかる。
resamples_fitted %>%
modeltime_resample_accuracy(summary_fns = mean) %>%
table_modeltime_accuracy(.interactive = FALSE)| Accuracy Table | |||||||||
| .model_id | .model_desc | .type | n | mae | mape | mase | smape | rmse | rsq |
|---|---|---|---|---|---|---|---|---|---|
| 1 | ARIMA(0,1,1)(0,1,1)[12] | Resamples | 3 | 296.32 | 2.84 | 1.02 | 2.91 | 337.44 | 0.90 |
| 2 | PROPHET | Resamples | 3 | 242.46 | 2.34 | 0.83 | 2.38 | 277.74 | 0.90 |
| 3 | PROPHET W/ XGBOOST ERRORS | Resamples | 3 | 239.41 | 2.32 | 0.82 | 2.35 | 277.30 | 0.89 |
| 4 | RANGER | Resamples | 3 | 494.54 | 4.79 | 1.69 | 4.83 | 555.23 | 0.55 |
この結果から、modeltime_refit関数を使って、Prophet
Boostモデルで、フルデータセットを再学習させて、将来の予測を行なう。ここでは1年先を予測する。
refit_tbl <- resamples_fitted %>%
filter(.model_id == 3) %>%
modeltime_refit(data = m750)
refit_tbl %>%
modeltime_forecast(h = "1 years", actual_data = m750) %>%
plot_modeltime_forecast(.interactive = FALSE)