Modeltime -Feature Engineering with Recipes-
UPDATE: 2022-12-28 11:52:30
はじめに
modeltimeパッケージを使いながら時系列データの予測を行なうノートをいくつかに分けてまとめている。modeltimeパッケージの目的は、時系列モデルの構築を効率よく行うためのパッケージで、tidymodelsパッケージと組み合わせて時系列モデルを作成する際に使用される。時系列予測のための体系的なワークフローを提供し、時系列ブーストモデルなど新しいアルゴリズムもカバーしている点が特徴的。モデルの数理的な側面や機械学習の用語などは、このノートでは扱わない。
とりあえず、modeltimeパッケージについて一通り使い方はまとめてたので、総括として、グローバルモデルを利用し、クロスバリデーション、ハイパーパラメータチューニングを行い、予測を複数の系列に対して行なう。
使用データ
必要なパッケージを読み込んでおく。
library(tsibbledata)
library(tsibble)
library(timetk)
library(tidyverse)
library(tidymodels)
library(modeltime)オーストラリアの小売の売上データセット(tsibbledata::aus_retail)を利用する。データの各時系列は下記のキーで一意に識別される。
State: オーストラリアの州Industry: 業種
ここではAustralian Capital Territoryの州の2つの業種restaurants and catering services、Clothing, footwear and personal accessory retailingに注目する。
monthly_retail_tbl <- aus_retail %>%
tk_tbl() %>%
filter(State == "Australian Capital Territory") %>%
filter(Industry %in% unique(aus_retail$Industry)[c(1,4)]) %>%
filter(Month >= as.Date("2013-01-01")) %>%
mutate(Month = as.Date(Month),
Industry = if_else(Industry == "Cafes, restaurants and catering services",
"catering", "clothing")
) %>%
mutate(Industry = as_factor(Industry)) %>%
select(date = Month, id = Industry, value = Turnover)
ids <- unique(monthly_retail_tbl$id)
DT::datatable(monthly_retail_tbl)2013-01-01から2018-12-01の72ヶ月分のデータで、2つの業種はきれいなパネルデータになっているのでこのまま利用する。
monthly_retail_tbl %>%
dplyr::group_by(id) %>%
dplyr::summarise(
n.obs = n(),
start_date = min(date),
end_date = max(date),
min_value = min(value),
mean_value = mean(value),
max_value = max(value)
)## # A tibble: 2 × 7
## id n.obs start_date end_date min_value mean_value max_value
## <fct> <int> <date> <date> <dbl> <dbl> <dbl>
## 1 catering 72 2013-01-01 2018-12-01 30.7 40.0 49.1
## 2 clothing 72 2013-01-01 2018-12-01 24.4 33.8 61.3このノートの全体像は下記のイメージである。イメージなので厳密ではない。
image
データの前処理
tidymodelsパッケージのrecipeでは(たぶん)グループごとに処理を行なうことが現状できないので、そのような処理はレシピを作成する前に行なう。ラグやウインドウ計算などは業種ごとにグループ化(フィルタリング)して処理しておく。ここで予測期間のデータにも一部の特徴量を混ぜることができる。
df <- map_dfr(
.x = 1:length(ids),
.f = function(x){
monthly_retail_tbl %>%
dplyr::filter(id == ids[x]) %>%
arrange(date) %>%
future_frame(date, .length_out = "12 months", .bind_data = TRUE) %>%
mutate(id = ids[x]) %>%
tk_augment_fourier(.date_var = date, .periods = 12, .K = 1) %>%
tk_augment_lags(.value = value, .lags = 12) %>%
tk_augment_slidify(.value = value_lag12,
.f = ~ mean(.x, na.rm = TRUE),
.period = c(3, 6, 9, 12),
.partial = TRUE,
.align = "center")
}
)
DT::datatable(df)データ分割
予測期間は2019-01-01から2019-12-01までの1年とする。
future_df <- df %>%
filter(is.na(value))
future_df %>%
dplyr::group_by(id) %>%
dplyr::summarise(
n.obs = n(),
start_date = min(date),
end_date = max(date)
)## # A tibble: 2 × 4
## id n.obs start_date end_date
## <fct> <int> <date> <date>
## 1 catering 12 2019-01-01 2019-12-01
## 2 clothing 12 2019-01-01 2019-12-01学習データは、2014-01-01から2018-12-01の60ヶ月分のデータ。予測期間にラグやウインドウ計算の特徴量を入れた際に、2013年は欠損するので除外した。
# training
preped_df <- df %>%
filter(!is.na(value)) %>%
drop_na()
preped_df %>%
dplyr::group_by(id) %>%
dplyr::summarise(
n.obs = n(),
start_date = min(date),
end_date = max(date),
min_value = min(value),
mean_value = mean(value),
max_value = max(value)
)## # A tibble: 2 × 7
## id n.obs start_date end_date min_value mean_value max_value
## <fct> <int> <date> <date> <dbl> <dbl> <dbl>
## 1 catering 60 2014-01-01 2018-12-01 30.7 39.9 49
## 2 clothing 60 2014-01-01 2018-12-01 24.9 34.5 61.3学習データを可視化しておく。cateringは緩やかに上昇し、そこからは停滞しており、takeawayでは特定期間に突発的なスパイクが存在している。
plot_time_series(
.data = preped_df,
.date_var = date,
.value = value,
.color_var = id,
.smooth = TRUE,
.interactive = FALSE) +
scale_x_date(limits = c(min(preped_df$date), max((preped_df$date))),
labels = date_format("%Y/%m"),
breaks = date_breaks("3 month")) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.2, hjust=0.2))
学習データを分割する。ここではcateringを対象に学習データとテストデータを可視化しているが、takeawayでも同様である。
splits <- time_series_split(
data = preped_df,
assess = 12,
cumulative = TRUE
)
splits %>%
tk_time_series_cv_plan() %>%
filter(id == "catering") %>%
plot_time_series_cv_plan(date, value, .interactive = FALSE) +
scale_x_date(limits = c(min(preped_df$date), max((preped_df$date))),
labels = date_format("%Y/%m"),
breaks = date_breaks("3 month")) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.2, hjust=0.2))
モデルのパラメタを調整するために、学習データを学習データと評価データに分割する。ここでは4スライス分の時系列クロスバリデーションデータを作成する。
resamples_kfold <- training(splits) %>%
time_series_cv(
date_var = date,
assess = "12 months",
cumulative = TRUE,
skip = "3 months",
slice_limit = 4
)
resamples_kfold## # Time Series Cross Validation Plan
## # A tibble: 4 × 2
## splits id
## <list> <chr>
## 1 <split [72/24]> Slice1
## 2 <split [66/24]> Slice2
## 3 <split [60/24]> Slice3
## 4 <split [54/24]> Slice4時系列クロスバリデーションデータを可視化しておく。
resamples_kfold %>%
tk_time_series_cv_plan() %>%
filter(id == "catering") %>%
plot_time_series_cv_plan(date, value, .interactive = FALSE) +
scale_x_date(limits = c(min(preped_df$date), max(assessment(resamples_kfold$splits[[1]])$date)),
labels = date_format("%Y/%m"),
breaks = date_breaks("3 month")) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.2, hjust=0.2))
時系列クロスバリデーションデータは2014-01-01を起点に、学習データと評価データを調整している。
resamples_kfold %>%
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(c(analysis_min_date, analysis_max_date, assessment_min_date, assessment_max_date)) %>%
select(id, analysis_min_date, analysis_max_date, assessment_min_date, assessment_max_date)## # A tibble: 4 × 5
## id analysis_min_date analysis_max_date assessment_min_date assessment_ma…¹
## <chr> <date> <date> <date> <date>
## 1 Slice1 2014-01-01 2016-12-01 2017-01-01 2017-12-01
## 2 Slice2 2014-01-01 2016-09-01 2016-10-01 2017-09-01
## 3 Slice3 2014-01-01 2016-06-01 2016-07-01 2017-06-01
## 4 Slice4 2014-01-01 2016-03-01 2016-04-01 2017-03-01
## # … with abbreviated variable name ¹assessment_max_dateレシピ・モデル・ワークフローの作成
最初に行った特徴量作成に加え、その他のカレンダー系特徴量を追加する。また、目的変数は対数変換を行なう。
recipe_prophet_boost <- recipe(value ~ ., data = training(splits)) %>%
step_timeseries_signature(date) %>%
step_rm(matches("(.iso$)|(.xts$)|(week)|(day)|(hour)|(minute)|(second)|(am.pm)|(.num$)|(.lbl.*$)")) %>%
step_mutate(date_month = factor(date_month, ordered = TRUE)) %>%
step_dummy(all_nominal(), one_hot = TRUE) %>%
step_log(value)
recipe_prophet_boost %>% prep() %>% bake(training(splits)) %>% DT::datatable()次はモデルの設定。今回はBoost
Prophetを利用する。mtryとsample_sizeは上限を決める必要があるので、ここはグリッド作成の際に対処する。おそらく、ここでの設定方法が誤っていたため、これまでのノートではBoost
Prophetのフィッティングの際に、ワーニングが発生していたと思われる。
model_prophet_boost <- prophet_boost(
# prophet
mode = "regression",
growth = "linear",
changepoint_num = tune(),
changepoint_range = tune(),
seasonality_yearly = TRUE,
seasonality_weekly = FALSE, # Monthlyデータなので機能しない
seasonality_daily = FALSE, # Monthlyデータなので機能しない
season = "additive",
prior_scale_changepoints = tune(),
prior_scale_seasonality = tune(),
# prior_scale_holidays = NULL,
# xgboost
mtry = tune(),
trees = 1000,
tree_depth = tune(),
learn_rate = tune(),
loss_reduction = tune(),
sample_size = tune(),
stop_iter = tune(),
min_n = tune()
) %>%
set_engine("prophet_xgboost")
model_prophet_boost## PROPHET Regression Model Specification (regression)
##
## Main Arguments:
## growth = linear
## changepoint_num = tune()
## changepoint_range = tune()
## seasonality_yearly = TRUE
## seasonality_weekly = FALSE
## seasonality_daily = FALSE
## season = additive
## prior_scale_changepoints = tune()
## prior_scale_seasonality = tune()
## mtry = tune()
## trees = 1000
## min_n = tune()
## tree_depth = tune()
## learn_rate = tune()
## loss_reduction = tune()
## sample_size = tune()
## stop_iter = tune()
##
## Computational engine: prophet_xgboost作成したモデルとレシピをワークフローに登録し、これを利用してハイパーパラメータのチューニングを行なう。
workflow_prophet_boost <- workflow() %>%
add_model(model_prophet_boost) %>%
add_recipe(recipe_prophet_boost)
workflow_prophet_boost## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: prophet_boost()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 5 Recipe Steps
##
## • step_timeseries_signature()
## • step_rm()
## • step_mutate()
## • step_dummy()
## • step_log()
##
## ── Model ───────────────────────────────────────────────────────────────────────
## PROPHET Regression Model Specification (regression)
##
## Main Arguments:
## growth = linear
## changepoint_num = tune()
## changepoint_range = tune()
## seasonality_yearly = TRUE
## seasonality_weekly = FALSE
## seasonality_daily = FALSE
## season = additive
## prior_scale_changepoints = tune()
## prior_scale_seasonality = tune()
## mtry = tune()
## trees = 1000
## min_n = tune()
## tree_depth = tune()
## learn_rate = tune()
## loss_reduction = tune()
## sample_size = tune()
## stop_iter = tune()
##
## Computational engine: prophet_xgboostハイパーパラメータチューニング
ここではラテンハイパーキューブを使って、グリッドを20個作成する。
set.seed(1989)
grid <- grid_latin_hypercube(
extract_parameter_set_dials(model_prophet_boost) %>%
update(mtry = mtry(range = c(1, 20)),
sample_size = sample_prop(c(0.4, 0.8))),
size = 20
)
grid %>% DT::datatable()作成したグリッドを利用してチューニングを行なう。本当はここで並列化の設定をしたほうが良いかもしれないが、そこらへんは雰囲気でしか理解してないので、また後日、学習した際にノートにまとめる。
tune_workflow_prophet_boost <- workflow_prophet_boost %>%
tune_grid(
resamples = resamples_kfold,
grid = grid,
control = control_grid(verbose = TRUE,
save_pred = TRUE),
metrics = metric_set(rmse, mae)
)
tune_workflow_prophet_boost## # Tuning results
## # NA
## # A tibble: 4 × 5
## splits id .metrics .notes .predictions
## <list> <chr> <list> <list> <list>
## 1 <split [72/24]> Slice1 <tibble [40 × 15]> <tibble [0 × 3]> <tibble [480 × 15]>
## 2 <split [66/24]> Slice2 <tibble [40 × 15]> <tibble [0 × 3]> <tibble [480 × 15]>
## 3 <split [60/24]> Slice3 <tibble [40 × 15]> <tibble [0 × 3]> <tibble [480 × 15]>
## 4 <split [54/24]> Slice4 <tibble [40 × 15]> <tibble [0 × 3]> <tibble [480 × 15]>mae基準で最小のパラメタの組み合わせを10個取り出す。
tune_workflow_prophet_boost %>%
show_best("mae", n = 10) %>%
DT::datatable()パラメタと評価指標の関係を可視化する。この際に、評価指標が下がる傾向があれば、パラメタグリッドを再作成するなどして再度チューニングする。
tune_workflow_prophet_boost %>%
autoplot() +
geom_smooth(se = FALSE) 
mae基準で最小のパラメタ組み合わせは下記の通り。
tuned_best_params <- tune_workflow_prophet_boost %>%
select_best("mae")
tuned_best_params %>% t() ## [,1]
## changepoint_num "24"
## changepoint_range "0.6598805"
## prior_scale_changepoints "0.03299381"
## prior_scale_seasonality "17.86857"
## mtry "17"
## min_n "8"
## tree_depth "11"
## learn_rate "0.02808174"
## loss_reduction "0.0002741032"
## sample_size "0.5998895"
## stop_iter "10"
## .config "Preprocessor1_Model09"一旦、これをベストなパラメタとして、ワークフローのパラメタをファイナライズする。
finalize_workflow_prophet_boost <- workflow_prophet_boost %>%
finalize_workflow(parameters = tuned_best_params)
finalize_workflow_prophet_boost$fit$actions$model$spec## PROPHET Regression Model Specification (regression)
##
## Main Arguments:
## growth = linear
## changepoint_num = 24
## changepoint_range = 0.659880460824352
## seasonality_yearly = TRUE
## seasonality_weekly = FALSE
## seasonality_daily = FALSE
## season = additive
## prior_scale_changepoints = 0.0329938065497755
## prior_scale_seasonality = 17.8685654750194
## mtry = 17
## trees = 1000
## min_n = 8
## tree_depth = 11
## learn_rate = 0.0280817387862989
## loss_reduction = 0.000274103153914346
## sample_size = 0.599889543387108
## stop_iter = 10
##
## Computational engine: prophet_xgboostモデルテーブルに登録するために、とりあえずフィッティングしておく。
# ここのデータはなんでもよい
fit_finalize_workflow_prophet_boost <- finalize_workflow_prophet_boost %>%
fit(training(resamples_kfold$splits[[1]]))
fit_finalize_workflow_prophet_boost## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: prophet_boost()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 5 Recipe Steps
##
## • step_timeseries_signature()
## • step_rm()
## • step_mutate()
## • step_dummy()
## • step_log()
##
## ── Model ───────────────────────────────────────────────────────────────────────
## PROPHET w/ XGBoost Errors
## ---
## Model 1: PROPHET
## - growth: 'linear'
## - n.changepoints: 24
## - changepoint.range: 0.659880460824352
## - yearly.seasonality: 'TRUE'
## - weekly.seasonality: 'FALSE'
## - daily.seasonality: 'FALSE'
## - seasonality.mode: 'additive'
## - changepoint.prior.scale: 0.0329938065497755
## - seasonality.prior.scale: 17.8685654750194
## - holidays.prior.scale: 10
## - logistic_cap: NULL
## - logistic_floor: NULL
##
## ---
## Model 2: XGBoost Errors
##
## xgboost::xgb.train(params = list(eta = 0.0280817387862989, max_depth = 11L,
## gamma = 0.000274103153914346, colsample_bytree = 1, colsample_bynode = 0.708333333333333,
## min_child_weight = 8L, subsample = 0.599889543387108), data = x$data,
## nrounds = 1000, watchlist = x$watchlist, verbose = 0, early_stopping_rounds = 10L,
## objective = "reg:squarederror", nthread = 1)モデルキャリブレーション
これをモデルテーブルに登録してキャリブレーションする。
calibrattion_prophet_boost <- modeltime_table(fit_finalize_workflow_prophet_boost) %>%
modeltime_calibrate(testing(splits), id = "id") # confを計算するために必要
calibrattion_prophet_boost %>%
pluck(".calibration_data", 1) %>%
DT::datatable()これでテストデータの精度を確認できるので、テストデータでグローバルモデルの精度を確認する。
# Global
calibrattion_prophet_boost %>%
modeltime_accuracy(acc_by_id = FALSE) %>%
table_modeltime_accuracy()こちらはローカルモデルの精度。
# Local
calibrattion_prophet_boost %>%
modeltime_accuracy(acc_by_id = TRUE) %>%
table_modeltime_accuracy()テストデータでの予測値はこの通り。
calibrattion_prophet_boost %>%
modeltime_forecast(
new_data = testing(splits),
actual_data = preped_df,
conf_by_id = TRUE,
keep_data = TRUE
) %>%
DT::datatable()可視化すると、そこまで精度が良くないように思えるが、今回はパッケージの使い方のまとめなので、これで良いとする。
calibrattion_prophet_boost %>%
modeltime_forecast(
new_data = testing(splits),
actual_data = preped_df,
conf_by_id = TRUE,
keep_data = TRUE
) %>%
group_by(id) %>%
plot_modeltime_forecast(
.interactive = FALSE,
.legend_show = FALSE,
.facet_ncol = 3
) +
scale_x_date(limits = c(min(preped_df$date), max(preped_df$date)),
labels = date_format("%Y/%m"),
breaks = date_breaks("3 month")) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.2, hjust=0.2))
再学習・予測
あとはモデルをフルデータで再学習して、予測期間の予測値を計算する。
refit_prophet_boost <- calibrattion_prophet_boost %>%
modeltime_refit(data = preped_df)
refit_prophet_boost## # Modeltime Table
## # A tibble: 1 × 5
## .model_id .model .model_desc .type .calibration_data
## <int> <list> <chr> <chr> <list>
## 1 1 <workflow> PROPHET W/ XGBOOST ERRORS Test <tibble [24 × 5]>おさらいしておくと、予測期間は2019-01-01から2019-12-01までの1年間である。
bind_rows(
training(splits) %>% group_by(id) %>% summarise(start_date = min(date), end_date = max(date), n.obs = n()) %>% mutate(label = "train"),
testing(splits) %>% group_by(id) %>% summarise(start_date = min(date), end_date = max(date), n.obs = n()) %>% mutate(label = "test"),
future_df %>% group_by(id) %>% summarise(start_date = min(date), end_date = max(date), n.obs = n()) %>% mutate(label = "future")
) %>%
select(id, label, everything()) %>%
arrange(id, start_date)## # A tibble: 6 × 5
## id label start_date end_date n.obs
## <fct> <chr> <date> <date> <int>
## 1 catering train 2014-01-01 2017-12-01 48
## 2 catering test 2018-01-01 2018-12-01 12
## 3 catering future 2019-01-01 2019-12-01 12
## 4 clothing train 2014-01-01 2017-12-01 48
## 5 clothing test 2018-01-01 2018-12-01 12
## 6 clothing future 2019-01-01 2019-12-01 12再学習したモデルで予測値を計算する。
refit_prophet_boost %>%
modeltime_forecast(
new_data = future_df,
actual_data = preped_df,
conf_by_id = TRUE
) %>%
DT::datatable()計算した予測値を使って、データ全体を可視化する。
refit_prophet_boost %>%
modeltime_forecast(
new_data = future_df,
actual_data = preped_df,
conf_by_id = TRUE
) %>%
# recipeでlog変換したものを戻す
mutate(
.value = exp(.value),
.conf_lo = exp(.conf_lo),
.conf_hi = exp(.conf_hi)
) %>%
group_by(id) %>%
plot_modeltime_forecast(
.title = "Turnover 1-year forecast",
.interactive = FALSE,
.facet_ncol = 2
) +
scale_x_date(limits = c(min(preped_df$date), max(future_df$date)),
labels = date_format("%Y/%m"),
breaks = date_breaks("3 month")) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.2, hjust=0.2))