はじめに
はじめまして．株式会社音圧爆上げくんにプロKagglerとして所属していますAshmeと申します．
業務の一環としてKaggleの様々なコンペティションに参加し，そこで得られた知見などを記事にして投稿しております．よろしくお願いいたします．
今回は今月中開催されているTabular Playground Series - Mar 2022で用いられていた特徴量エンジニアリングについて紹介します．前回の記事で私がベースラインで使用した特徴量については紹介しましたが，今回はKaggleのノートブックで見た特徴量エンジニアリングについて紹介します．
また特徴量エンジニアリングだけでなく後処理が非常に効果的だったので，こちらについても紹介します．内容としては単なるクリッピングだったのですが，これを行うか行わないかだけでMAEが0.1ほど違いました．
コンペティションの説明についても以前の記事にありますが，念の為こちらにも載せておきます．こちらのコンペティションは時系列データの予測をする回帰問題になっています．ただ時間的に変化する特徴量はターゲットの変数のみとなっています．
モデル改善のためには？
突然ですが皆さんはモデル改善といえば何を思い浮かべるでしょうか？使用するモデルを変える，はいパラーパラメータを調整するなどもあると思います．しかし，モデルの改善の中でもこれらは効果的とは言えないのではないでしょうか．より良い改善の方法としてはEDAによりデータの特徴を捉え，特徴量エンジニアリングによってより効果的な特徴量を作成することです．
ですが今回はせっかくのコンペティションであり，様々な方のコードを参考にすることができるため他の参加者の方のコードを参考に特徴量を新たに追加していこうと思います．
今回のノートブックはこちらにありますので全文が気になる方はこちらをご一読ください．また参考にさせていただいたものはノートブックにも記載しておりますが，記事の最後に参考として載せておきます．
追加した特徴量
前回のベースラインでは以下のような特徴量を作成しました．
- 日時データを三角関数を用いて変換
- ラグ特徴量
- 集約特徴量
これらに加えて今回は以下の特徴量を追加しました．
- 午前の時間帯かを判別するフラグ
- 休日（土日）かどうかを判別するフラグ
- その日付が1年のうちの何日目か，何週目に存在するか
またこれらに加えて参考にしたノートブックにあった時間特徴量の変換について紹介していきます．
午前の時間帯か判別するフラグ
時系列データでは一点の時間だけでなく，時間帯によって特徴があらわれることもありえます．例えばスーパーの各時間の売上を予測するタスクを考えてみます．そうすると，開店直後，閉店直前などはお客さんの入りが悪く売上はあまり伸びませんが，お昼の時間帯や夕方の時間帯など買い物に来る人が増える時間帯，いわゆるピークが存在することが考えられると思います．
こちらの特徴量はそのような考え方と同様なものだと思います．実装としては，pd.to_datetimeでdatetime型に変換した後，時間を取得し，その時間を用いて午前中（今回は6時〜12時）の間を1，それ以外を0とした特徴量を作成します．
df["time"] = pd.to_datetime(df["time"])
df["hour"] = df["time"].dt.hour
df["am"] = ((df["hour"] < 12) & (df["hour"] > 6)).astype("int8")
1行目，2行目は前回の時間特徴量の作成と同様です．3行目で6時〜12時までのものをTrue，それ以外がFalseに変換し，それをint型にすることでTrueを1，Falseを0に変換するようになっています．
土日かどうかを判別するフラグ
前述の例に沿って，次は曜日によって売上が変化するかどうかを考えてみます．平日と土日，どちらのほうが売上は伸びやすいでしょうか？アルバイトをしたことがある人はイメージが容易だと思いますが，土日のほうが忙しいですし，シフトとしても出てほしいと言われやすいと思います．このように曜日によっても売上が大きく変化することが考えられます．
これを利用して平日と土日を判別するフラグを作成します．こちらもpd.to_datetimeでdatetime型に変換した後，曜日情報を取得し，土曜日(5)と日曜日(6)のデータには1，それ以外のものは0とした特徴量を作成します．
df["weekday"] = df["time"].dt.weekday
df["Saturday"] = (df["weekday"] == 5).astype("int8")
df["Sunday"] = (df["weekday"] == 6).astype("int8")
こちらも1行目で曜日情報を取得してweekday列を作成しています．そしてこの列の値が5のものはSaturday列を1，6のものはSunday列を1とすることで平日，土曜日，日曜日をそれぞれ判別できるようになっています．
日付にさらなる情報を付与する
こちらはおそらくデータに時系列性を明確に持たせるためのものではないかと思います．
1つ目はある日付が1年のうちの何日目かという情報です．例えば1月1日は年の初めであり，1年全体のなかでも1日目です．例えば他に3月1日は1年全体のなかでは60日目（閏年の場合は61日目）というような情報を特徴量として作成します．
2つ目はある日付が1年のうちの何週目にあるかという情報です．先ほどと同様に考えれば，1月1日は1週目であり，3月1日は9週目，もしくは10週目にあります．これを特徴量として作成します．
こちらも実装は難しいものではなく，pd.to_datetimeによってdatetime型に変換した後，何日目かはdayofyear，何周目かはweekによって作成することができます．
df["day"] = df["time"].dt.dayofyear
df["week"] = df["time"].dt.week
それ以外の工夫
時間特徴量の変換
私のベースラインでは時間特徴量に周期性を持たせるために，三角関数を用いて変換しました．ただ参考にさせていただいたものにまた異なった変換があったため，こちらについて紹介します．
ノートブックでは'moment'という列名で作成されていました．言葉で説明するのが少し難しいので，まずはコードを見ていきましょう．
df["time"] = (df["time"].dt.hour-12)*3 + df["time"].dt.minute/20
まず第1項を見てみましょう．時間から12を引くことで午前のデータは負の値，午後のデータは正の値になります．次に3をかけているのは，今回のデータでは観測するのが1時間に3回だからでしょう．
第2項は分を20で割っています．これも同様で，観測が1時間に3回なので，60分を3で割った20で割っています．
この2つを足すことで，1日分の時間データを-36から35までの値に変換することができます．また1回目の観測と2回目の観測の間は1となっており，単調増加のためうまく時系列性と周期性を捉えることができていると考えられます．
図にするとこれらの特性を持っていることがわかるのではないでしょうか．これを図にしたものを以下に示します．
これは3日分の先程の変換した特徴量をプロットしたものです．確かに1日で-36から35になり，次の日になると-36に戻るというように順番と周期性を確認することができます．
それ以外にも参考ノートブックで使用されていたものとして，集約特徴量や，ラグ特徴量がありましたが，集約する際に参照する列が異なるだけだったため記事では割愛します．ノートブックにすべて載せていますので，気になる方はそちらをご覧ください．
Feature Importances
最後に学習したLightGBMのFeature Importancesを可視化します．以下にコードと結果を載せます．
# get feature importances
feature_importances = model.feature_importances_
# make indices to sort
indices = np.argsort(feature_importances)
# sort features importances
feature_importances = feature_importances[indices]
# sort feature names
sort_columns = np.array(columns)[indices]
# plot feature importances
plt.figure(figsize=(12, 12))
plt.barh(sort_columns, feature_importances)
plt.show()
ここで紹介した特徴量はあまり上の方にありませんが，参考に作成した特徴量(median_cong, rolling_7_std)などが上位に来ています．やはりラグ特徴量や集約特徴量などは重要度としても高いことがわかります．これらは非常に有効ですが，作成するときにはリークに注意を払うようにしましょう．
後処理
本当はスコアを上げるために特徴量エンジニアリングに力を入れようと思ったのですが，参考にしたノートブックに非常に効果的な後処理があったため紹介します．
最初にも書いたとおり，内容としてはクリッピングをしています．参考ノートブックそのままの実装ですが，テストデータが含まれる月（9月）の訓練データに存在するターゲットの15パーセンタイル値と70パーセンタイル値によってクリッピングするというものです．
sep = trans_train[trans_train['month'] >= 9]
lower = sep.groupby(['time', 'x', 'y', 'direction'])["congestion"].quantile(0.15).values[2340:]
upper = sep.groupby(['time', 'x', 'y', 'direction'])["congestion"].quantile(0.7).values[2340:]
clip_submission = submission.copy()
clip_submission["congestion"] = submission["congestion"].clip(lower, upper)
最後にスライスを用いているのは時間をテストデータのものに合わせるためです．1行目で訓練データの中の9月のデータのみを抽出しています．
2行目，3行目で時間と座標+方向（前回記事のものだとregionにあたる）ごとに15パーセンタイル値と70パーセンタイル値を取得しています．
最後にモデルで予測した値をclipを用いて15パーセンタイル値と70パーセンタイル値の間になるようにクリッピングしています．
私の場合はこれをするかしないかだけでLeaderBoardのスコアが0.1ほど差が出たので，今回のコンペティションでは非常に有用な後処理と言えるのではないでしょうか．
祝日データの消去
考え方としては土日にフラグを立てることと似ているのではないでしょうか．参考にしたノートブックの1つでは祝日のデータを削除していました．おそらく，今回のデータでは1年もないデータであるため祝日は外れ値として考えるようにしたのだと思います．
こちらも今回は参考ノートブックのものをほぼそのままで，祝日かどうかをフラグとして立てるという形で実装しました．調べたところPythonのライブラリにholidaysというライブラリがあるそうです．様々な国の祝日データをこちらのライブラリで参照できるようです．
train["time"] = pd.to_datetime(train["time"])
train["official_holiday"] = train["time"].dt.date.astype(str).str.contains('1991-05-27|1991-07-04|1991-09-02').astype('int')
test["official_holiday"] = 0
おわりに
今回はTabular Playground Series - Mar 2022で用いられていた特徴量エンジニアリングや後処理について紹介しました．
ドメイン情報がなくても時系列データは様々な特徴量が作成できるなと感じたのと同時に，どこかで見たことがあるものだったのでベースラインといえどもっと考えながらコードを書いてくべきだと身にしみました．後処理に関しても非常に効果的なものだったので，今後のコンペティションでも参考にしていきたいと思います．
最後に人材募集となりますが，株式会社音圧爆上げくんではプロKagglerを募集しています．
少しでも興味の有る方はぜひ以下のリンクをご覧の上ご応募ください．
Wantedlyリンク
参考
[1] https://www.kaggle.com/code/kotrying/tps-2022-03-lgbm
[2] https://www.kaggle.com/code/ifashion/tps-mar2022-single-lgbm-lb-4-91-run-time-91s