【Kaggle】Tabular Playground Series - Mar 2022のBaseline

はじめに

はじめまして．株式会社音圧爆上げくんにプロKagglerとして所属していますAshmeと申します．

業務の一環としてKaggleの様々なコンペティションに参加し，そこで得られた知見などを記事にして投稿しております．よろしくお願いいたします．

今回は今月中開催されているTabular Playground Series - Mar 2022のベースラインとして作成したモデルについて説明していきます．モデルとしてLightGBMを用いて，時系列データに対する基本的な特徴量エンジニアリング，簡単なアンサンブルについても説明していきます．今回のコードについてはこちらにあります．気になる方はぜひご一読ください．

こちらのコンペティションは時系列データの予測をする回帰問題になっています．通常のテーブルデータでは用いないような時系列データに対する特徴量エンジニアリングが必要になります．

またTabular Playground Seriesは毎月開催されており，通常のコンペティションとは異なり，テーブルデータが対象となっており難易度もそこまで高くないため，コンペティションに慣れていない方でも楽しむことができると思います．

特徴量エンジニアリング

データを見るとわかりますが，今回のデータは時系列データですが，時間で変わるものは観測した日時と，ターゲットのみになっています．またカテゴリ変数が多く，数値変数はターゲットのみとなっています．

日時データの変換

まず日時データに関する変換について説明します．日時データを利用するときはPandasを用いてdatetime型に変換すると月，日などを取得するのが簡単になるため最初にdatetime型に変換します．

これはpd.to_datetime(pdはpandasの略称）を用いることで簡単に変換することができます．

    df['time'] = pd.to_datetime(df['time'])
    df['month'] = df['time'].apply(lambda x: x.month)
    df['day'] = df['time'].apply(lambda x: x.day)
    df['hour'] = df['time'].apply(lambda x: x.hour)
    df['minute'] = df['time'].apply(lambda x: x.minute)

一番上の行でpd.to_datetimeを用いることで時間データをdatetime型に変換しています．それより下の行では月，日，時間，分をそれぞれ特徴量として抽出しています．

これはdatetime型に変換した後にmonth, day, hour, minuteを取得することで簡単にできます．またyearを用いることで年も取得することができますが，今回のデータはすべて同年のデータのため作成していません．

LightGBMなどの決定木を用いたアルゴリズムでは数値の大小関係のみを予測に用いるためこのままでも特に問題ありません．ただこれだけでは周期性を持っていません．

月は1月から12月までしかなく，12月の次は1月になります．時間，分などに関しても同様です．

このような周期性を持たせるための変換として三角関数を用いたものがあります．コサイン，サインは角度によって-1から1を取り，360度回転するとサインは0，コサインは1に戻るようになっています．これを利用して特徴量が周期性を持つように変換してみます．

    # それぞれ三角関数を用いて周期性をもたせる
    # 期間内で4月，6月，9月は30日しかないため少し処理を変えている
    month_list = [4, 6, 9]
    df[['day_sin', 'day_cos']] = 0
    for month in df['month'].unique():
        if month in month_list:
            df.loc[df['month'] == month, 'day_sin'] =\
                df.query("month == @month")['day'].apply(lambda x: np.sin((2*np.pi*x) / 30))
            df.loc[df['month'] == month, 'day_cos'] =\
                df.query("month == @month")['day'].apply(lambda x: np.cos((2*np.pi*x) / 30))
        else:
            df.loc[df['month'] == month, 'day_sin'] =\
                df.query("month == @month")['day'].apply(lambda x: np.sin((2*np.pi*x) / 31))
            df.loc[df['month'] == month, 'day_cos'] =\
                df.query("month == @month")['day'].apply(lambda x: np.cos((2*np.pi*x) / 31))

    df['month_sin'] = np.sin(2*np.pi*df['month'] / 12)
    df['month_cos'] = np.cos(2*np.pi*df['month'] / 12)
    df['hour_sin'] = np.sin(2*np.pi*df['hour'] / 24)
    df['hour_cos'] = np.cos(2*np.pi*df['hour'] / 24)
    df['minute_sin'] = np.sin(2*np.pi*df['minute'] / 60)
    df['minute_cos'] = np.cos(2*np.pi*df['minute'] / 60)

まず変数に2πをかけ，変数の最大値で割ります．これによりすべての値が0から2πの間に変換されます．また大小関係はそのまま保っています．

次にそれをサイン，コサインを用いて変換することで角度として近いものは近い値になり，最小値と最大値も近い値になります．これにより周期的に近い特徴量が値としても近いものになります．

例えば1月と12月はただの数値として1や12とだけでは値として遠いものになります．しかし，実際には12月の次は1月と周期性を持っており近いものになっています．上記のような変換を用いることでこのような特徴量を持たせることができると考えられます．

三角関数を用いて変換する手法は私もあまり見たことがありませんが，周期性をもたせることで効果を発揮するデータに対しては有効な可能性があるため，知っておいて損はないのではないでしょうか．

ラグ特徴量

時系列データでよく用いられる特徴量エンジニアリングとしてラグ特徴量があります．これは対象の時刻よりも前の特徴量を入力とするものです，同時刻の情報だけでなく，前時刻の特徴量を用いることができるため時系列的に関連する予測をするときに効果的であると考えられます．

ラグ特徴量はpd.DataFrameに対してshiftメソッドを用いることで簡単に実装することができます．

    region_df["lag_1"] = region_df["congestion"].shift(1)

ここでのラグ特徴量は観測された座標（x, y）と方向（direction）ごとにターゲットをずらしています．

例として出したものは特徴量を1つずらすことで1つ前の時刻のターゲットを入力としています．shiftに対する引数の分だけずらすことができるため，関連のある時刻のデータを入力にすることができます．

今回のデータに対しては20分毎にデータを観測しているため，20分前，40分前，1時間前，2時間前といった直近のデータに加えて1日前のデータをラグ特徴量として作成しています．

またただずらすだけでなく，ある期間のデータの統計量を特徴量として作ることもできます．こちらはrollingメソッドを用いることで実装することができます．

    region_df[f"lag_avg{i}"] = region_df["congestion"].rolling(window=i, min_periods=1).mean().shift(1)

この場合はある時刻のi個前から1つ前のデータまでの平均値を計算しています．そのままだと対象の時刻を含んだi個のデータの平均値を特徴量としてしまいリークしてしまいます．そのため１つずらすために最後にshift(1)を入れています．

これと同様の方法で標準偏差も計算しています．

集約特徴量

集約特徴量はある条件ごとにデータを分割し，それらの平均などを取ったものです．今回は先程も利用したregionと，曜日ごとに大まかな傾向があると考えられるためregion+weekdayごとのターゲットの平均と中央値，標準偏差を作成しました．

まずはregionとregion+weekdayごとのターゲットの平均や中央値，標準偏差を計算する部分です．

    region_mean_map = train_df.groupby("region")["congestion"].mean()
    region_median_map = train_df.groupby("region")["congestion"].median()
    region_std_map = train_df.groupby("region")["congestion"].std()

    region_weekday_mean_map = train_df.groupby("region_weekday")["congestion"].mean()
    region_weekday_median_map = train_df.groupby("region_weekday")["congestion"].median()
    region_weekday_std_map = train_df.groupby("region_weekday")["congestion"].std()

コード内のregion_weekdayはregion+weekdayによって作成した列です．groupbyを用いることでregionやregion_weekdayごとに平均，中央値，標準偏差を簡単に計算することができます．

次に計算した値で新しい列を作成します．

    train_df["region_target_mean"] = train_df["region"].map(region_mean_map)
    train_df["region_target_median"] = train_df["region"].map(region_median_map)
    train_df["region_target_std"] = train_df["region"].map(region_std_map)

    train_df["region_weekday_target_mean"] =\
        train_df["region_weekday"].map(region_weekday_mean_map)
    train_df["region_weekday_target_median"] =\
        train_df["region_weekday"].map(region_weekday_median_map)
    train_df["region_weekday_target_std"] =\
        train_df["region_weekday"].map(region_weekday_std_map)

mapメソッドを用いることで各region，region_weekdayに対応する平均値，中央値，標準偏差に置き換えたものを新しい列として作成しています．

エクスパンディング特徴量

エクスパンディング特徴量とはある時点よりも前のすべての期間の集計料を特徴量としたものです．私も時系列データの特徴量エンジニアリングで調査し，今回初めて聞きましたので利用してみました．

今回はある時点よりも前のターゲットすべてのターゲットの中央値を特徴量として作成しています．コード中のmedianをmeanに変えるだけで平均値も作成することができます．

    region_df.loc[i, "expand_congestion"] = region_df.loc[:i-1, "congestion"].median()

私のコードでは愚直に計算するように実装していますが，pandasにはexpanding()というメソッドもあるようです．より簡潔な実装はこちらにありましたので参考にしていただければと思います．

One-hot encoding

最後に時系列データのみに対するものではないですが，カテゴリ変数をone-hot encodingしています．実装は様々な方法が有りますが，個人的にはpd.get_dummies関数を用いるのが簡単だともいます．実装は以下の通りです．

    tmp_train = pd.get_dummies(trans_train[column], prefix=column, drop_first=False)
    trans_train = pd.concat((trans_train, tmp_train), axis=1)

ここではone-hot encodingを使用しましたが，決定木がもとのモデルではlabel encodingでも問題ないと聞いたことも有ります．

私としては線型回帰やニューラルネットにも適用できるone-hot encodingのほうが好ましいためこちらをよく使用しています．

特徴量エンジニアリングは以上になります．実装ではそれぞれの特徴量の作成を関数にしていき，最終的にすべてをまとめたものをpreprocessing関数としています．また作成した新しい列をすべてfeature_columnsとしてまとめてモデルに渡す時に使用しています．

基本的にはすべてを1つの関数で実装するのではなく，機能ごとに関数，クラスの場合はメソッドを作成していくことでデバッグや修正が簡単にできるのでこのようにするのが通常だと思います．

CVの測定

submissionを作成する前に訓練データでどれほどの性能があるか測定します．今回は訓練データを3つ(train, validation, test)に分割し，trainとvalidationを用いてハイパーパラメータの調整，testを用いることでモデルの性能を評価しています．

まずは先程作成したpreprocessing関数を用いることで訓練データに新しい特徴量を作成します．その後，今回はテストデータが訓練データよりあとの12時間分のデータ（2340サンプル）だったため，同様の条件にするためにvalidationデータをtrainデータのあとの2340サンプル，testデータをvalidationデータのあとの2340サンプルとして分割しています．

train_df, val_df, test_df =\
    trans_train.iloc[:trans_train.shape[0]-4680, :],\
    trans_train.iloc[trans_train.shape[0]-4680:trans_train.shape[0]-2340, :],\
    trans_train.iloc[trans_train.shape[0]-2340:, :]

また数値型の特徴量に対する前処理（標準化など）は今回使用していません．決定木は主に特徴量の大小関係でノードが分岐するため，標準化などによってスケールを合わせても合わせなくてもほとんど変わらないためです．

その後ハイパーパラメータの調整をoptunaを用いて行っています．optunaはベイズ最適化によってハイパーパラメータを最適化してくれるライブラリです．ほとんどのモデルで利用できるため知っておくと良いと思います．

optunaを用いてtrainデータでモデルを訓練し，validationデータでモデルを評価して，validationデータに対する誤差（今回のコンペティションではMAE）が最小になるようにハイパーパラメータを調整しています．

def objective(trial):
    # 最適化するパラメータの指定
    params = {
        "num_leaves": trial.suggest_int("num_leaves", 10, 500),  # 2のmax_depth乗が良いらしい
        "max_depth": trial.suggest_int("max_depth", 1, 20),
        "learning_rate": trial.suggest_loguniform("learning_rate", 1e-4, 1e-1),
        "n_estimators": trial.suggest_int("n_estimators", 100, 1000),
        "min_child_samples": trial.suggest_int("min_child_samples", 10, 100),
        "reg_alpha": trial.suggest_loguniform("reg_alpha", 1e-3, 1e-1),
        "reg_lambda": trial.suggest_loguniform("reg_lambda", 1e-3, 1e-1),
        "metric": "mae"
    }

    # モデルの訓練
    feature_columns = dummy_columns + num_columns + use_columns
    train_X, val_X = train_df[feature_columns], val_df[feature_columns]
    train_y, val_y = train_df[target], val_df[target]

    model = lgb.LGBMRegressor(**params)
    model.fit(train_X, train_y,
              eval_set=(val_X, val_y),
              early_stopping_rounds=10, verbose=0)
    pred = model.predict(val_X)

    # 評価
    score = mean_absolute_error(val_y, pred)
    return score


study = optuna.create_study()
study.optimize(objective, 100)
print(study.best_params)
print(study.best_value)

その後，optunaで最適化されたハイパーパラメータを用いて最適化に使用していないtestデータを用いて汎化性能を測定しています．

# モデルの訓練
feature_columns = dummy_columns + num_columns + use_columns
train_X, val_X = train_df[feature_columns], val_df[feature_columns]
train_y, val_y = train_df[target], val_df[target]

test_X, test_y = test_df[feature_columns], test_df[target]

model = lgb.LGBMRegressor(**study.best_params)
model.fit(train_X, train_y,
          eval_set=(val_X, val_y),
          early_stopping_rounds=10, verbose=10,
          eval_metric="mae")
pred = model.predict(test_X)

score = mean_absolute_error(test_y, pred)
print(f"Validation Score: {score}")

記事を書いている途中に気づきましたが，testデータに対する性能を測定するときはtrainデータとvalidationデータを結合してそれらを用いてモデルを学習したほうが良かったかなと思います．後ほどNotebookを修正しようと思います．

Submissionの作成

これらによってモデルの汎化性能を知ることができたため，ついにsubmissionを作成します．

まずはCVの測定でも行ったのと同様にtrainデータとtestデータ（ここではKaggleでのtrain.csvとtest.csvを指す）に対してpreprocessing関数を用いて特徴量を作成します．

train_df = pd.read_csv(os.path.join(DIR, "train.csv"))
test_df = pd.read_csv(os.path.join(DIR, "test.csv"))

train_df, test_df = preprocessing(train_df, test_df)

for column in cat_columns:
    tmp_train = pd.get_dummies(train_df[column], prefix=column, drop_first=False)
    tmp_test = pd.get_dummies(test_df[column], prefix=column, drop_first=False)

    train_df = pd.concat((train_df, tmp_train), axis=1)
    test_df = pd.concat((test_df, tmp_test), axis=1)

    dummy_columns.extend(tmp_train.columns.values.tolist())

for column in train_df.columns:
    if column not in test_df.columns:
        test_df[column] = 0

print(f"train data shape: {train_df.shape}")
print(f"test data shape: {test_df.shape}")

その後モデルを学習するのですが，今回はLightGBMのモデルを8つ作成し平均値や中央値を予測とするような簡単なアンサンブル学習を用います．まずはモデルがそれぞれ異なる学習をするようにrandom_stateを指定するのですが，先に乱数を作成しています．ここでは乱数を作成していますが，再現性を持たせるために事前8つ数字を指定すべきだったと思います．

random_state_list = np.random.randint(0, 1e+4, size=8)

次にそれぞれのモデルを学習します．学習したモデルはmodel_dictに格納し，予測を行う時に呼び出せるようにしています．

model_dict = dict()

train_X = train_df[feature_columns]
train_y = train_df[target]

for i in range(8):
    random_state = random_state_list[i]
    model = lgb.LGBMRegressor(**study.best_params, random_state=random_state)
    model.fit(train_X, train_y)

    model_dict[f"model{i}"] = model

これでモデルの学習が終わりました．ついに予測をしていくのですが，このままでは少し問題が有ります．

今回はターゲットのラグ特徴量を作成しています．ラグ特徴量では1つ前の時刻のターゲットも使用していますが，testデータではこれを作成することができていません．なのでtestデータすべてを一気にモデルに渡して予測するということができないようになっています．

そこで予測を時系列順に行っていきます．特徴量エンジニアリングでも作成したregion列を再度使用し，regionごとのDataFrameを作成します．そしてテストデータの1番目から予測を行っていき，その予測結果を用いてラグ特徴量などを作成することでそれ以降のtestデータでもラグ特徴量などを使用できるようになっています．これを実装したのが以下になります．

submission_dict = dict()

region_list = test_df["region"].unique()
for i in range(8):
    submission = pd.DataFrame(columns=["row_id", "congestion"])
    model = model_dict[f"model{i}"]
    for region in region_list:
        # regionごとにDataFrameを抽出
        train_region_df = train_df.query("region == @region")
        test_region_df = test_df.query("region == @region")
        train_size = train_region_df.shape[0]

        region_df = pd.concat((train_region_df, test_region_df)).reset_index(drop=True)

        # 時系列順に予測，予測結果を用いてラグ特徴量を計算
        for j in range(test_region_df.shape[0]):
            target_id = train_size + j

            region_df.loc[target_id, "lag_1"] = region_df.loc[target_id-1, "congestion"]
            region_df.loc[target_id, "expand_congestion"] = region_df.loc[:target_id-1, "congestion"].median()
            for k in [2, 3, 6, 72]:
                region_df.loc[target_id, f"lag_{k}"] = region_df.loc[target_id-k, "congestion"]
                region_df.loc[target_id, f"lag_avg{k}"] = region_df.loc[target_id-k:target_id-1, "congestion"].mean()
                region_df.loc[target_id, f"lag_std{k}"] = region_df.loc[target_id-k:target_id-1, "congestion"].std()

            # 予測
            pred = model.predict(region_df.loc[target_id, feature_columns].values.reshape(1, -1))[0]

            # 予測した値を提出ファイルに格納
            submission = submission.append({"row_id": [region_df.loc[target_id, "row_id"]][0],
                                            "congestion": pred}, ignore_index=True)
            # ラグ特徴量の計算のためにDataFrameにも格納
            region_df.loc[target_id, target] = pred

    submission = submission.sort_values("row_id").reset_index(drop=True)
    submission["row_id"] = submission["row_id"].astype(int)
    submission_dict[f"submission{i}"] = submission

submission_dictはそれぞれのモデルの予測結果を格納します．コード中のregion_dfがregionごとに抽出したデータであり，時系列順に並べてあります．

事前にregion_dfのうちのtrainデータの数をtrain_sizeとして保持しておき，それよりあとのtestデータに対して予測をするようにしています．

実装では予測するインデックスを取得→そのインデックスを用いてラグ特徴量，集約特徴量，エクスパンディング特徴量を作成→予測→結果の格納というようになっています．

これで8つのモデルによる予測結果ができたので，あとは平均値や中央値を取ってアンサンブルするだけとなります．

またKaggleのNotebookを見ていると，元のターゲットが整数だったため，予測結果を整数にするといったこともされていたためこれも取り入れてみました．Notebookでは中央値を予測としたものにしか適用していませんが，同様に平均値を予測としたものに対しても適用することができます．

これでsubmissionを作成するところまで終わりました．今回のコンペティションはコードコンペティション（KaggleのNotebookを提出する形式のもの）ではないため，ここまでで作成したcsvファイルを提出することで完了となります．

今回は結果としてはLeaderBoardの真ん中よりほんのちょっと上に位置するくらいになっていますが，時系列順に予測するためより長い時間に対する予測だと予測による誤差が大きくなっていく可能性が考えられます．また予測→特徴量の生成といった処理になるためデータ量が多い場合はより時間がかかる可能性が有ります．

おわりに

今回は毎月開催されているTabular Playground Seriesで作成したベースラインモデルについて紹介しました．こちらのコンペティションはテーブルデータを用いたコンペティションであるため，Kaggleの参加経験が少ない方の練習にもなると思いますのでぜひ参加してみてはいかがでしょうか．

最後に人材募集となりますが，株式会社音圧爆上げくんではプロKagglerを募集しています．
少しでも興味の有る方はぜひ以下のリンクをご覧の上ご応募ください．
Wantedlyリンク