時系列予測モデルの実践論４　－昔のKaggleを事例に－

サトヤキ・ロイ　池野心平　鈴木雄太

GBTなどの機械学習を用いたグローバルモデルでの時系列予測

なぜグローバルモデルが必要なのか

　皆さんこんにちは。NRIの時系列予測チームのデータサイエンティスト、Satyakiです。今回も時系列モデルについて、機械学習(ML)モデルを中心にお話します。前回のブログでは、個々の時系列データを繰り返し学習・予測し、その結果をまとめて全体の予測を行う方法を説明しました。個別の時系列ごとに学習を行うことには多くの利点があります。例えば、他のデータソースから影響を受けにくいこと、デバッグのしやすさ、また柔軟性が高い、などが挙げられます。さらに、前回のブログで紹介したように、個々の時系列の過去期間すべてを訓練データとして扱えることもメリットと言えます。
しかし、グローバルモデルにも多くのメリットがあります。

最先端の機械学習やディープラーニングのアルゴリズムを使用することが可能
地域、店舗、SKUなどの階層的なデータに含まれる特徴を、複数の時系列を統合した販売履歴全体のデータから表現できる
外生変数を特徴量として使用し、正則化などの技術によりオーバーフィッティングを低減することが可能。ディープラーニングモデル等においても、様々なデータ構造をモデルの特徴量として使用できる
個々に作られた数千のモデルではなく、単一のモデルを管理すればよい

　また近年では予測性能でも、グローバルモデルが従来の統計モデルに勝るケースも見られます。例えば、2020年に行われたM5コンペティションでは上位5位までのチーム　全てがLightGBMを使用していました。私たちの最終的なモデルでも、アンサンブルモデルの主要な構成要素としてLightGBMを使用しています。もちろんグローバルモデルにも問題はありますが（後の頁で説明）、予測性能と時間短縮の点で最も有効な予測手法と言えると思います。

前処理について

　モデリングの前に、いくつかの重要な前処理について説明します。１．欠損値の取り扱い　２．外れ値の取り扱い　についてです。

１．欠損値の取り扱い

　本コンペでは、2012-11-02から39週間先の予測をします。学習データの最小販売日は2010-02-05ですが、すべての店舗や部門でこの期間すべての販売データがあるわけではありません。すべてのstoreとdeptのペアにおいて、訓練データの最小日から最大日（予測開始時点の1週間前）までの日付の欠損を埋めていきます。

　グラフからも読み取れる通り、いくつかの店舗や部門で大きなデータの欠損があることが読み取れます。例えば店舗34と部門78においては欠損値が142あり、学習データにおいて1つのデータポイントしかないということになります（本来は店舗とデパートの双方において143週分のデータが格納されているはずです。)。
データの欠損の理由はわからないので、次の2つの手法を用いて欠損値を補完します。

０で埋める（欠損期間は販売が0と想定）
季節性を考慮した補完と繰り越し

以下が季節性を考慮した補完の方法です。

  impute <- function(current_ts){
    current_ts = ts(current_ts, start = decimal_date(min(training_max_date)), frequency = 52)
    if(sum(!is.na(current_ts)) >= 3){
      as.numeric(na_seadec(current_ts))
    } else if(sum(!is.na(current_ts)) == 2){
      as.numeric(na_interpolation(current_ts))
    } else{
      as.numeric(na_locf(current_ts))
    }
  }

モデル予測の分散をより大きくするために両方のアプローチを使用します。

２．外れ値の取り扱い

まず、いくつかの系列にはマイナス値が入っているようです。

　これは、バックオーダー(在庫がなく入荷待ちの商品を受注すること)やキャンセルされた注文によるものかもしれませんが、やはり理由が明確ではないので、単純にゼロに置き換えます。データの正当性をチェックするだけでなく、従属変数がより正規化され、外れ値が少なくなるように調整する必要があります。これは機械学習(ML)モデルにおいて必ず必要な処理ではありませんが、経験上、前記を改善した方が性能向上に繋がります。今回、外れ値を除去こそはしませんでしたが、データの分散を減らし、より正規化することを試みました。相乗平均、対数処理、Box-Cox変換(変数を正規分布に従わせるような変換)などの一般的なデータ変換を試しましたが、最終的な予測誤差の観点から相乗平均を用いることが最も良い結果をもたらすことが分かりました。

　ご覧の通り、変換により分布がより平坦になっていることがわかります。

ハイパーパラメータの最適化

　機械学習モデルでは最適なハイパーパラメータを設定が重要です。柔軟すぎるパラメータはオーバーフィットにつながり、頑健すぎる値はモデルの性能不足につながる可能性があります。また、時系列データを扱っているため、学習-検証におけるデータセットの選択にも注意が必要です（これについては、パート１で詳しく説明しています）。一般的に、pythonのハイパーパラメータのチューニングライブラリはより発展しておりかつ柔軟なので、我々はこの課題にpython、特にscikit-optimizeパッケージを使用しました(具体的にはscikit-optimizeパッケージで広範囲なハイパーパラメータ検索ルーチンを実行し、Rによる機械学習APIでそれらのパラメータを使用)。以下はそのコードブロックですが、必要に応じて一般化することが可能です。

    param_space = [
        space.Integer(10,200, name="max_depth"),
        space.Integer(100, 2000, name="n_estimators"),
        space.Real(0.001,0.05, prior="log-uniform", name="learning_rate"),
        space.Integer(100,1000, name="max_bin"),
        space.Real(0.2,0.8, prior="uniform", name="feature_fraction"), ##Not in catboost
        space.Real(0.1,1, prior="uniform", name="subsample"),
        space.Real(1e-9, 100.0, prior="log-uniform", name="reg_lambda"),
        space.Real(1e-9, 100.0, prior="log-uniform", name="reg_alpha"), ##Not in catboost
        space.Integer(10,1000, name="num_leaves"),
        space.Categorical(['regression'], name="objective"),
    ]

    param_names = [
        "max_depth",
        "n_estimators",
        "learning_rate",
        "max_bin",
        "feature_fraction",
        "subsample",
        "reg_lambda",
        "reg_alpha",
        "num_leaves",
        "objective"
    ]

    def optimize(params, param_names, feats_df):
        params  = dict(zip(param_names, params))
        common_params = {
            'boosting_type': 'gbdt',
            'verbosity': -1,
             'objective': 'regression',
             'metric': 'rmse' 
            }
        params = {**params, **common_params} 
        rmse = []
        for fold in range(n_folds):
            xtrain_all =  feats_df[feats_df.kfold != fold].reset_index(drop=True)
            xtrain = xtrain_all.drop(['Date','Weekly_Sales',"kfold","id"], axis=1)
            ytrain = xtrain_all['Weekly_Sales']

            xtest_all = feats_df[feats_df.kfold == fold].reset_index(drop=True)
            xtest = xtest_all.drop(['Date','Weekly_Sales',"kfold","id"], axis=1)
            ytest = xtest_all['Weekly_Sales']

            dtrain = lgb.Dataset(xtrain, label=ytrain, categorical_feature=object_cols)
            model = lgb.train(params, dtrain, verbose_eval=100)
            preds = model.predict(xtest)
            fold_rmse = metrics.mean_squared_error(ytest, preds, squared=False)
            rmse.append(fold_rmse)
            return np.mean(rmse)


    def train_model(train_current_features):

    #     overdone_control = DeltaYStopper(delta=0.0001) 
        optimization_function = partial(optimize, param_names=param_names, feats_df = train_current_features)
        result = gp_minimize(
                optimization_function,
                dimensions=param_space,
                n_calls=50,
                n_random_starts=10,
                n_initial_points= 10,
                verbose=False
            )
        plot_convergence(result)
        best_params = (dict(zip(param_names, result.x)))
        print(best_params)
        return best_params

　ガウス過程(Gaussian Processes)を用いたベイズ最適化のためのインターフェイスを提供するgp_minimizeパッケージを使用しています。パッケージの動作詳細については、以下のリンクを参照してください。(https://scikit-optimize.github.io/stable/index.html)
今回のケースでは探索に数時間を要しました。何回か繰り返すことにより有効なパラメータのリストや検索範囲を絞り込むことができるはずです。

モデリングに用いたアプローチ

　モデル精度の改善を継続的に回すために、次のようなプロセスを作成しました。この方法は、何千もの販売履歴の時系列モデルを一括で作成する場合に有効です。

ステップ１では、ロバストなアンサンブルモデルを作るための基本として、バイアス/バリアンスの観点で学習能力や強みの異なった個別モデルを訓練します。これにより最終的なアンサンブル(複合的な)予測における分散を小さくすることができます。前回のブログでご覧になったように、我々は複数の統計モデルを開発しました。今回は、私たちがどのようにMLモデルを開発・進化させたかをご紹介します。
ステップ２では、モデルごとのパフォーマンスサマリーテーブルを作成し、全体精度をチェックします。このステップでは、パフォーマンスの悪いモデルや学習に時間がかかりすぎるモデルを削除することができます。次に、個々の時系列について、エントロピー、安定性、トレンド、スパイク、季節的強度などの時系列特性（特徴量）を作成します。次に、各特徴の寄与度を作成し、時系列特徴の寄与度ごとに各モデルのパフォーマンスのサマリーテーブルを作成します。このテーブルを分析することで、各時系列モデルの相対的な強さ/弱さについての示唆を得ることができます。
ステップ３では、ステップ2で見つかったモデルの脆弱性を補強できるような新たな特徴量やモデルを作成します。
最後に、テストデータによりアンサンブル(複合)モデルの精度をチェックします。

後処理

　最後に、クリスマスにおける予測に必要な調整を行いました。
　今回のデータセットでは、クリスマスのある週は、2010年は0日、2011年は1日、2012年は3日の祝前日があります。つまり、2012年から2010年にかけては3日、2012年から2011年にかけては2日の差があることになります。1週間（7日間）で、平均2.5日ですから、この値分の影響を第51週から第52週に移すことにより、モデルが考慮しなかった分を手動で補います。

検討されたモデル

LightGBM, and trend-adjusted LightGBM
Catboost
XGBoost
Random Forest
Seasonal Auto-ARIMA [Discussed in last blog post] 7.Error, Trend, Seasonal (ETS)
model [Discussed in last blog post]
TBATS [Discussed in last blog post]

RのparnsipインターフェースはほとんどのMLタスクに標準的なインターフェースを提供しているので，そちらを使用しました．以下はその例です。

    object_cols_num <- which(names(training_data) %in% object_cols)

    model_spec_lgbm_tuned <- boost_tree(
      mode = "regression",
      ## New parameters after hyperparameter tuning of sqrt transformed Weekly_Sales
      mtry           = 16,
      trees          = 1771,
      tree_depth     = 172,
      learn_rate     = 0.0265657835139171,
      sample_size    = 0.643235802797972

    ) %>% set_engine("lightgbm", objective = "rmse", max_bin=972,reg_lambda=0.000000001,
                     reg_alpha = 0.000537792755360404, num_leaves=952, boosting_type = "gbdt",categorical_feature = object_cols_num)

    wflw_fit_lgbm <- workflow() %>%
      add_model(
        spec = model_spec_lgbm_tuned
      ) %>%
      add_recipe(recipe_spec_lgb) %>% 
      fit(training_data)

モデル性能の比較

評価スコアとランキング

Model	Private Leaderboard Score	Rank	Rank after adjustment
LightGBM	2585	19	17
Catboost	2742	37	32
XGBoost	2950	112	127
Random Forest	3008	135	139

　アルゴリズムによって、スコアに差があることがわかります。最近の時系列コンペティションでは主にLightGBMが採択されていますが、今回も、XGBoostやRandom Forestは、LightGBMやCatboostのような最新のアルゴリズムと比較すると、パフォーマンスが劣っていました。

予測結果の比較

　いくつかの予測モデルを可視化して確認し整合性をチェックしてみましょう。

店舗34_部門65

店舗43_部門60

店舗20_部門59

　見ての通り、XGBoostの予測はかなりノイズが多く、分散も大きいので、このモデルは学習データに対してオーバーフィットしていると言えます。したがって、このモデルをアンサンブルモデルから除外することにします。

モデルの弱点 –ツリー系モデルにおけるトレンドの取り扱い

　ツリーベースの機械学習モデルの最大の弱点は、強いトレンドを捉えるのが難しいという点です。機械学習モデルはデータを分割して予測を行うため、テストデータの分布がトレーニングデータと大きく異なる場合があり、強いトレンドがる場合その影響を捉えられません。
　例えばこの系列では徐々に売上が減少しています。

　この傾向を捉えるために、トレンド調整済みLightGBMというモデルを追加で作成しました。設計は以下の通りです。

prophetモデルに学習データを学習させます。トレンドと転換点におけるオーバーフィットを避けるために頑健なパラメータを使用します。
トレンド成分を抽出し、メインのモデルから差し引きます。残差、つまりトレンド調整済の売上とトレンド成分の予測値を得ます。
残差に対してLightGBMモデルを学習させ、検証期間の残差成分を予測します。
Prophetから予測されたトレンド成分とLightGBMから予測された残差を加算し、最終的な予測値を得ます。

　　以下は、トレンド調整済みLightGBMによる予測値です。

　グラフから読み解けるように、予測値はかなり低めとなり、下降傾向をとらえています。これらの予測は低すぎたり高すぎたりするので単独で使用すべきではありませんが、アンサンブル化(複合化)することによって全体の精度が高まります。また、prophetと機械学習モデルのそれぞれのハイパーパラメータを一緒に最適化することは難しく、prophetのパラメータを都度修正する必要があります。

　アンサンブルの詳細については次回の記事で紹介しますが、LightGBMとTrend-adjusted LightGBMの単純平均だけで、ランク3を達成できました。

Model	Private Leaderboard Score	Rank	Rank after adjustment
LightGBM	2585	19	17
Trend-adjusted LightGBM	2724	36	34
Simple Average	2308	-	3

まとめ

　このブログでは、機械学習の観点からデータ処理からモデル構築までの全体的なモデリングプロセス、続いて、モデル改善と分析のプロセスをご説明しました。LightGBMのようなトップクラスの機械学習ライブラリを使用するだけでもトップ20にランクインできました。また、異なるモデルの長所と短所を組み合わせることで、より精度の良いアンサンブルモデルを作ることができることを紹介しました。次回の投稿では、より効果的なアンサンブルの方法についてご紹介できればと思います。