UPDATE: 2022-12-26 18:28:37

はじめに

tidymodelsパッケージの使い方をいくつかのノートに分けてまとめている。tidymodelsパッケージは、統計モデルや機械学習モデルを構築するために必要なパッケージをコレクションしているパッケージで、非常に色んなパッケージがある。ここでは、今回はmodeltimeパッケージで複数の時系列を予測する方法についてまとめていく。モデルの数理的な側面や機械学習の用語などは、このノートでは扱わない。

下記の公式ドキュメントやtidymodelsパッケージに関する書籍を参考にしている。

modeltimeパッケージの目的

modeltimeパッケージは、時系列モデルの構築を効率よく行うためのパッケージで、tidymodelsパッケージと組み合わせて時系列モデルを作成する際に使用される。時系列予測のための体系的なワークフローを提供し、時系列ブーストモデルなど新しいアルゴリズムもカバーしている点が特徴的。

ここでは、複数の時系列データがある場合に効率よく予測モデルを構築する方法をまとめる。

上記のドキュメントやサイトをなぞりながら、不明点を追記しながらまとめていく。

グローバルモデリング

多数の時系列を予測を行う場合、方法としては2つある。1つ目は、各時系列に対してモデルを作成し、繰り返して予測を行なう方法。2つ目は、時系列分野での呼び方なのか、modeltimeパッケージ特有の言葉なのかわからないが、GlobalModelを用いる方法。これは、すべての時系列の予測を生成する1つのモデルのみが作成され、このモデルを用いて、各時系列の予測を行なう。双方、メリット・デメリットは存在するものの、スケーラビリティを考えるのであれば、後者の方法が望ましい。

例えば、精度に関しては1つ目の方法の方が優れている可能性があるが、数百の店舗がある場合に、数百のモデルを作成する必要があり、時間がかかる、多くのメモリが必要になる可能性がある。

グローバルモデリングの文脈ではパネルデータが重要になる。パネルデータは各種様々な指標に関して、複数の系列が同じ間隔で記録されているデータのフォーマット。パネルデータの文脈で先程のことを考えると、パネルが100個あれば、100個分のモデル作成を行い、100個分の予測を繰り返すことになる。グローバルモデリングでは、100個分のパネルに対して、1つのモデルを作成し、100個分の予測を行なう。

ただ、グローバルモデルは特徴量エンジニアリングがうまくできないと、精度が低くなってしまう。そのため、特徴量エンジニアリングが非常に重要になる。

モデルの作成

まずは必要なパッケージを読み込んでおく。

library(tidymodels)
library(modeltime)
library(tidyverse)
library(timetk)
library(DT)

使用するデータは、modeltimeパッケージのwalmart_sales_weeklyデータ。StoreとDeptを識別するid、週開始日date、週売上valueを利用する。

data_table <- walmart_sales_weekly %>%
    select(id, Date, Weekly_Sales) %>%
    set_names(c("id", "date", "value"))

data_table
## # A tibble: 1,001 × 3
##    id    date        value
##    <fct> <date>      <dbl>
##  1 1_1   2010-02-05 24924.
##  2 1_1   2010-02-12 46039.
##  3 1_1   2010-02-19 41596.
##  4 1_1   2010-02-26 19404.
##  5 1_1   2010-03-05 21828.
##  6 1_1   2010-03-12 21043.
##  7 1_1   2010-03-19 22137.
##  8 1_1   2010-03-26 26229.
##  9 1_1   2010-04-02 57258.
## 10 1_1   2010-04-09 42961.
## # … with 991 more rows

idのユニーク数は7個なので、ウォルマートの7店舗分の2010-02-05から2012-10-26までの週売上が143レコードづつ記録されている。

data_table %>% 
  group_by(id) %>% 
  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: 7 × 7
##   id    n.obs start_date end_date   min_value mean_value max_value
##   <fct> <int> <date>     <date>         <dbl>      <dbl>     <dbl>
## 1 1_1     143 2010-02-05 2012-10-26    14537.     22513.    57592.
## 2 1_3     143 2010-02-05 2012-10-26     6166.     13150.    51159.
## 3 1_8     143 2010-02-05 2012-10-26    31061.     35718.    42664.
## 4 1_13    143 2010-02-05 2012-10-26    32782.     38693.    44863.
## 5 1_38    143 2010-02-05 2012-10-26    55974.     79978.   127812.
## 6 1_93    143 2010-02-05 2012-10-26    55082.     71699.    97980.
## 7 1_95    143 2010-02-05 2012-10-26    93359.    120772.   148798.

plot_time_series関数を利用すると、一気に売上推移を可視化できる。これをみると、売上規模として、4つくらいに分類できそうで、どの時系列も特有の特徴を持っていることがわかる。つまり、各時系列に特有の特徴量がなければ、予測が難しそうなことがわかる。

data_table %>%
  #group_by(id) %>%
  plot_time_series(
    date, 
    value, 
    .color_var = id,
    .smooth = FALSE,
    .interactive = FALSE, 
    .facet_ncol = 3
  )

時系列を分けて可視化したいときはグループ化すればよい。

data_table %>%
  group_by(id) %>%
  plot_time_series(
    date, 
    value, 
    .color_var = id,
    .interactive = FALSE, 
    .facet_ncol = 3
  )

ドキュメントに従って、time_series_split関数を利用して、データを学習セットとテストセット(3 か月分)に分割する。重複するタイムスタンプが検出された。スライディングウィンドウを使って、重なり合った時系列を一緒に処理する、というインフォメーションが表示される。

splits <- data_table %>% 
  time_series_split(
    assess     = "3 months", 
    cumulative = TRUE
  )
df_train <- splits %>% training()
df_test <- splits %>% testing()

bind_rows(
  df_train %>% summarise(start_date = min(date), end_date = max(date), records = n()),
  df_test %>% summarise(start_date = min(date), end_date = max(date), records = n())
)
## # A tibble: 2 × 3
##   start_date end_date   records
##   <date>     <date>       <int>
## 1 2010-02-05 2012-08-03     917
## 2 2012-08-10 2012-10-26      84

グラフ化しておく。

splits %>%
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(date, value, .interactive = FALSE,
                           .title = "Time Series Validation Plan") + 
  scale_x_date(limits = c(as.Date("2010-02-01"), as.Date("2012-11-01")),
               labels = date_format("%Y/%m"),
               breaks = date_breaks("3 month")) + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.2, hjust=0.2))

モデリングのためのレシピを作成する。このレシピが行っている処理を理解するために、必要な列だけ残して可視化する。大枠として、2010-02-05のレコードを対象にすると、2010-02-05がidの分だけ積み上げられており、id_X*_*というカラムでレコードが管理されていることがわかる。

recipe(value ~ ., df_train) %>%
  step_mutate_at(id, fn = droplevels) %>%
  step_timeseries_signature(date) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE) %>%
  prep() %>%
  bake(df_train) %>%
  select(date, value, date_week, contains("id")) %>% 
  DT::datatable()

このような形で、各idごとに917/7=131レコードを積み上げているデータを利用する。

recipe(value ~ ., df_train) %>%
  step_mutate_at(id, fn = droplevels) %>%
  step_timeseries_signature(date) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE) %>%
  prep() %>%
  bake(df_train) %>%
  select(date, value, date_week, contains("id")) %>% 
  filter(id_X1_1 == 1) %>% 
  datatable()

この後で利用するモデルでは、dateカラムは不要なので、先程のレシピからstep_rm(date)で削除しておく。

recipe <- recipe(value ~ ., df_train) %>%
    step_mutate_at(id, fn = droplevels) %>%
    step_timeseries_signature(date) %>%
    step_rm(date) %>%
    step_zv(all_predictors()) %>%
    step_dummy(all_nominal_predictors(), one_hot = TRUE)

summary(prep(recipe)) %>% print(n = 50)
## # A tibble: 38 × 4
##    variable          type      role      source  
##    <chr>             <list>    <chr>     <chr>   
##  1 value             <chr [2]> outcome   original
##  2 date_index.num    <chr [2]> predictor derived 
##  3 date_year         <chr [2]> predictor derived 
##  4 date_year.iso     <chr [2]> predictor derived 
##  5 date_half         <chr [2]> predictor derived 
##  6 date_quarter      <chr [2]> predictor derived 
##  7 date_month        <chr [2]> predictor derived 
##  8 date_month.xts    <chr [2]> predictor derived 
##  9 date_day          <chr [2]> predictor derived 
## 10 date_mday         <chr [2]> predictor derived 
## 11 date_qday         <chr [2]> predictor derived 
## 12 date_yday         <chr [2]> predictor derived 
## 13 date_mweek        <chr [2]> predictor derived 
## 14 date_week         <chr [2]> predictor derived 
## 15 date_week.iso     <chr [2]> predictor derived 
## 16 date_week2        <chr [2]> predictor derived 
## 17 date_week3        <chr [2]> predictor derived 
## 18 date_week4        <chr [2]> predictor derived 
## 19 date_mday7        <chr [2]> predictor derived 
## 20 id_X1_1           <chr [2]> predictor derived 
## 21 id_X1_3           <chr [2]> predictor derived 
## 22 id_X1_8           <chr [2]> predictor derived 
## 23 id_X1_13          <chr [2]> predictor derived 
## 24 id_X1_38          <chr [2]> predictor derived 
## 25 id_X1_93          <chr [2]> predictor derived 
## 26 id_X1_95          <chr [2]> predictor derived 
## 27 date_month.lbl_01 <chr [2]> predictor derived 
## 28 date_month.lbl_02 <chr [2]> predictor derived 
## 29 date_month.lbl_03 <chr [2]> predictor derived 
## 30 date_month.lbl_04 <chr [2]> predictor derived 
## 31 date_month.lbl_05 <chr [2]> predictor derived 
## 32 date_month.lbl_06 <chr [2]> predictor derived 
## 33 date_month.lbl_07 <chr [2]> predictor derived 
## 34 date_month.lbl_08 <chr [2]> predictor derived 
## 35 date_month.lbl_09 <chr [2]> predictor derived 
## 36 date_month.lbl_10 <chr [2]> predictor derived 
## 37 date_month.lbl_11 <chr [2]> predictor derived 
## 38 date_month.lbl_12 <chr [2]> predictor derived

レシピが完成したので、モデルとワークフローを作成し、学習データでフィッティングする。パラメタチューニングなどはここでは一旦スルーする。

model <- boost_tree("regression") %>% 
  set_engine("xgboost")

workflow <- workflow() %>%
    add_model(model) %>%
    add_recipe(recipe)

workflow_fit <- workflow %>% 
  fit(df_train)

workflow_fit
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: boost_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 5 Recipe Steps
## 
## • step_mutate_at()
## • step_timeseries_signature()
## • step_rm()
## • step_zv()
## • step_dummy()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## ##### xgb.Booster
## raw: 46.5 Kb 
## call:
##   xgboost::xgb.train(params = list(eta = 0.3, max_depth = 6, gamma = 0, 
##     colsample_bytree = 1, colsample_bynode = 1, min_child_weight = 1, 
##     subsample = 1), data = x$data, nrounds = 15, watchlist = x$watchlist, 
##     verbose = 0, nthread = 1, objective = "reg:squarederror")
## params (as set within xgb.train):
##   eta = "0.3", max_depth = "6", gamma = "0", colsample_bytree = "1", colsample_bynode = "1", min_child_weight = "1", subsample = "1", nthread = "1", objective = "reg:squarederror", validate_parameters = "TRUE"
## xgb.attributes:
##   niter
## callbacks:
##   cb.evaluation.log()
## # of features: 37 
## niter: 15
## nfeatures : 37 
## evaluation_log:
##     iter training_rmse
##        1     46315.130
##        2     33003.670
## ---                   
##       14      3554.521
##       15      3423.419

ここからmodeltimeのワークフローを利用する。フィッティングが完了しているモデルをモデルテーブルに登録する。

model_table <- modeltime_table(
    workflow_fit
)
model_table
## # Modeltime Table
## # A tibble: 1 × 3
##   .model_id .model     .model_desc
##       <int> <list>     <chr>      
## 1         1 <workflow> XGBOOST

モデルテーブルに登録したあとは、キャリブレーションを行なう。id別に処理する必要があるので、引数でidを設定する。ここでのキャリブレーションは、モデルをテストデータでフィッティングし、テストデータで予測値を計算している。

calibration_table <- model_table %>%
    modeltime_calibrate(
      new_data = df_test, 
      id = "id"
    )
calibration_table
## # Modeltime Table
## # A tibble: 1 × 5
##   .model_id .model     .model_desc .type .calibration_data
##       <int> <list>     <chr>       <chr> <list>           
## 1         1 <workflow> XGBOOST     Test  <tibble [84 × 5]>

.calibration_dataカラムの中身を見ると、各idの日付ごとに予測値が算出されていることがわかる。

calibration_table %>% 
  pluck(".calibration_data", 1)
## # A tibble: 84 × 5
##    date       .actual .prediction .residuals id   
##    <date>       <dbl>       <dbl>      <dbl> <fct>
##  1 2012-08-10  16120.      19091.     -2971. 1_1  
##  2 2012-08-10  28257.      31012.     -2755. 1_3  
##  3 2012-08-10  37270.      35225.      2045. 1_8  
##  4 2012-08-10  42241.      39167.      3074. 1_13 
##  5 2012-08-10  74483.      69858.      4625. 1_38 
##  6 2012-08-10  84119.      78213.      5907. 1_93 
##  7 2012-08-10 137408.     130517.      6891. 1_95 
##  8 2012-08-17  17331.      19091.     -1760. 1_1  
##  9 2012-08-17  31906.      40312.     -8406. 1_3  
## 10 2012-08-17  36294.      35225.      1069. 1_8  
## # … with 74 more rows

idを1_1のレコードに絞ってみるとわかりやすい。

calibration_table %>% 
  pluck(".calibration_data", 1) %>% 
  filter(id == "1_1")
## # A tibble: 12 × 5
##    date       .actual .prediction .residuals id   
##    <date>       <dbl>       <dbl>      <dbl> <fct>
##  1 2012-08-10  16120.      19091.     -2971. 1_1  
##  2 2012-08-17  17331.      19091.     -1760. 1_1  
##  3 2012-08-24  16286.      18658.     -2372. 1_1  
##  4 2012-08-31  16680.      19449.     -2769. 1_1  
##  5 2012-09-07  18322.      20331.     -2009. 1_1  
##  6 2012-09-14  19616.      20332.      -715. 1_1  
##  7 2012-09-21  19252.      19899.      -648. 1_1  
##  8 2012-09-28  18948.      19109.      -161. 1_1  
##  9 2012-10-05  21904.      27505.     -5601. 1_1  
## 10 2012-10-12  22764.      20332.      2432. 1_1  
## 11 2012-10-19  24185.      23751.       434. 1_1  
## 12 2012-10-26  27391.      23751.      3639. 1_1

予測値の算出が終わったので、テストデータに対するグローバルモデルとローカルモデルの精度を確認する。modeltime_accuracy関数のacc_by_idがFALSEだとグローバルモデルの精度、TRUEだとローカルモデルの精度が計算される。

bind_rows(
     calibration_table %>% modeltime_accuracy(acc_by_id = FALSE) %>% mutate(id = "Global"),
     calibration_table %>% modeltime_accuracy(acc_by_id = TRUE)
     ) %>% 
  select(id, everything())
## # A tibble: 8 × 10
##   id     .model_id .model_desc .type   mae  mape  mase smape  rmse     rsq
##   <chr>      <int> <chr>       <chr> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
## 1 Global         1 XGBOOST     Test  3699.  8.30 0.114  8.17 5059. 0.982  
## 2 1_1            1 XGBOOST     Test  2126. 10.8  1.59  10.3  2597. 0.500  
## 3 1_3            1 XGBOOST     Test  3397. 18.7  0.571 17.6  4048. 0.912  
## 4 1_8            1 XGBOOST     Test  2087.  5.45 0.950  5.62 2284. 0.671  
## 5 1_13           1 XGBOOST     Test  1596.  3.86 0.703  3.96 1914. 0.529  
## 6 1_38           1 XGBOOST     Test  8345. 10.8  1.03  11.0  9686. 0.00787
## 7 1_93           1 XGBOOST     Test  4279.  5.37 0.420  5.57 5440. 0.754  
## 8 1_95           1 XGBOOST     Test  4060.  3.14 0.513  3.18 4875. 0.655

グローバルモデルの精度とは、idを考慮しない予測値と実測値の評価を行っており、

calibration_table %>% 
  pluck(".calibration_data", 1) %>% 
  summarise(
    mae = mean(abs(.actual - .prediction))
  )
## # A tibble: 1 × 1
##     mae
##   <dbl>
## 1 3699.

ローカルモデルの精度とは、idを考慮した予測値と実測値の評価を行っている。

calibration_table %>% 
  pluck(".calibration_data", 1) %>% 
  group_by(id) %>% 
  summarise(
    mae = mean(abs(.actual - .prediction))
  )
## # A tibble: 7 × 2
##   id      mae
##   <fct> <dbl>
## 1 1_1   2126.
## 2 1_3   3397.
## 3 1_8   2087.
## 4 1_13  1596.
## 5 1_38  8345.
## 6 1_93  4279.
## 7 1_95  4060.

グローバル精度、ローカル精度を確認する意図としては、どの時系列でモデルがうまく機能しているのか、反対に、どのモデルがうまく機能していないのかを判断できる。例えば、今回のモデルにおいては、1_38はうまく予測できていないが、1_13はうまく予測できていることがわかる。

各idの予測値をmodeltime_forecast関数で可視化しておく。calibration_tableにも予測値は入っているが、この関数に予測したいテストデータ、フルデータを渡すことで、観測値と予測値をあわせて可視化できる。

calibration_table %>%
    modeltime_forecast(
        new_data    = df_test,
        actual_data = data_table,
        conf_by_id  = TRUE
    ) %>%
    group_by(id) %>%
    plot_modeltime_forecast(
        .facet_ncol  = 3,
        .interactive = TRUE
    )

グローバルモデルがいったんうまく機能していると仮定して、modeltime_refit関数でフルデータで再学習を行なう。

refit_table <- calibration_table %>%
  modeltime_refit(data = data_table)

refit_table
## # Modeltime Table
## # A tibble: 1 × 5
##   .model_id .model     .model_desc .type .calibration_data
##       <int> <list>     <chr>       <chr> <list>           
## 1         1 <workflow> XGBOOST     Test  <tibble [84 × 5]>

この関数を利用することで、モデルテーブルのモデルの更新を行っている。見た目は全く同じなのでわかりくいが、内部ではデータが更新されていることがわかる。なぜかうまく値を取り出せないので、コピペしておく。上がcalibration_tableで、下がrefit_table。predictorsのレコード数が変わっていることがわかる。

# calibration_table
str(calibration_table)
mdl_tm_t [1 × 5] (S3: mdl_time_tbl/tbl_df/tbl/data.frame)
 $ .model           :List of 1
  ..$ :List of 4
  .. ..$ pre    :List of 3
  .. .. ..$ mold        :List of 4
  .. .. .. ..$ predictors: tibble [917 × 37] (S3: tbl_df/tbl/data.frame)

# refit_table
str(refit_table)
mdl_tm_t [1 × 5] (S3: mdl_time_tbl/tbl_df/tbl/data.frame)
 $ .model_id        : int 1
 $ .model           :List of 1
  ..$ :List of 4
  .. ..$ pre    :List of 3
  .. .. ..$ mold        :List of 4
  .. .. .. ..$ predictors: tibble [1,001 × 37] (S3: tbl_df/tbl/data.frame)

フルデータで再学習してたので、将来の日付に対する予測を行なう。まずは、将来の日付データをfuture_frame関数で作成する。

future_table <- data_table %>%
  group_by(id) %>%
  future_frame(
    .length_out = 52, # weeklyデータなので52週=1年
    .bind_data = FALSE
    ) %>% ungroup()

bind_rows(
  df_train %>% summarise(start_date = min(date), end_date = max(date), records = n()) %>% mutate(id = "train"),
  df_test %>% summarise(start_date = min(date), end_date = max(date), records = n()) %>% mutate(id = "test"),
  future_table  %>% summarise(start_date = min(date), end_date = max(date), records = n()) %>% mutate(id = "future")
) %>% 
  select(id, everything())
## # A tibble: 3 × 4
##   id     start_date end_date   records
##   <chr>  <date>     <date>       <int>
## 1 train  2010-02-05 2012-08-03     917
## 2 test   2012-08-10 2012-10-26      84
## 3 future 2012-11-02 2013-10-25     364

modeltime_forecast関数の挙動をidを限定して可視化しておく。後ろの方に予測値がついたデータフレームが作成されている。

refit_table %>%
  modeltime_forecast(
    new_data    = future_table,
    actual_data = data_table, 
    conf_by_id  = TRUE
  ) %>% 
  filter(id == "1_1") %>% 
  DT::datatable()

このデータフレームをplot_modeltime_forecast関数で可視化しおく。

refit_table %>%
  modeltime_forecast(
    new_data    = future_table,
    actual_data = data_table, 
    conf_by_id  = TRUE
  ) %>%
  group_by(id) %>%
  plot_modeltime_forecast(
    .interactive = FALSE,
    .facet_ncol  = 2
  )

これでグローバルモデルを利用したモデリングの基礎的な部分は終わり。下記はおまけ。

refit_table %>%
  modeltime_forecast(
    new_data    = future_table,
    actual_data = data_table, 
    conf_by_id  = TRUE
  ) %>% 
  ggplot(aes(.index, .value, col = id)) + 
  geom_line() + 
  geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi), alpha = 0.1) +
  geom_vline(xintercept = as.Date("2012-11-02")) + 
  theme_bw() + 
  scale_x_date(limits = c(as.Date("2010-02-01"), as.Date("2013-10-25    ")),
               labels = date_format("%Y/%m"),
               breaks = date_breaks("3 month")) + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.2, hjust=0.2))