FDUA第2回金融データ活用チャレンジ 1~位解法

# 1. 概要

この記事は、Signateで2024年1月18日から2月15日の期間に開催されたコンペ「第2回金融データ活用チャレンジ」における優勝チームの解法を紹介するものです。

ただし、このコンペでは、「開催最終日の夕方ごろからサーバーの混雑によって投稿用ページに誰もアクセスできず、参加者が最終評価用の提出ファイルを変えられない状態のまま終了を迎えてしまう」というトラブルが発生しました。平均F値という評価指標の特性も勘案すれば、このトラブルが最終順位に与えた影響は無視できないものであったと思います。

私たちのチームは、幸いにも暫定的に提出ファイルを選んでいたため、このトラブルからそれほど大きな被害を被らずに済み、結果だけを述べればPrivate LB、Public LBともに1位となる予測を提出することができました。上述の経緯により、私たちは100％の確信をもって「優勝」を誇れる状態にはありませんが、それでもこの最終結果について、チーム一同、大変嬉しく思っています。

以下、このコンペに対する私たちの解法と戦略について、本節でその概要を示したうえで、後続の節で前処理、モデル構築、後処理の方針をご紹介します。

**解法の概要** （角括弧内の数字は予測に含まれるラベル1の数＝陽性数）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1484754/05388c9f-92d7-5d08-105a-598f77501f74.png)

### 1-1. 前処理
- データの特徴量をもとに、日付・金額の加工やカテゴリ変数の結合などを行うことによって、合計26個のベース特徴量を作成しました。
- 各ベース特徴量について、水準ごとに集約したときの完済率が明らかに平均から乖離するような水準を取り出し、それらの重要水準を対象としてOne-Hotエンコーディングを施すことで、約300個のバイナリ変数を追加しました。

### 1-2. モデル構築
- 前処理済みのデータを用いて複数の分類モデルを学習させ、その中でLocal CVとPublic LBでの成績が比較的良好だった二つの分類モデル（CatBoostとLightGBM）の確率予測をロジット変換した値を入力変数として、ロジスティック回帰によるスタッキングモデルを作成しました。

### 1-3. 後処理
- スタッキングモデルが出力した予測順位をもとに、予測の陽性数（`MIS_Status=1`と予測する行数）が38,400から39,000の範囲で50個刻みの値になるように、提出用ファイルを複数作成しました。
- 本コンペの評価指標であった平均F値の性質を考慮し、私たちのチームでは「2つの最終提出を陽性数の多いものと少ないものに分ける」という方針を想定していました。最終日、暫定的に2つを選択しておいたうえで、Public LBの様子を見て夜に最終選択を行う予定でしたが、冒頭で述べた経緯により、暫定選択していた2つ（予測完済数が38,500個のものと38,800個のもの）を提出することになりました。

# 2. 前処理

このコンペのデータの重要な特徴のひとつとして、「どの変数にも、その変数が特定の値であるかどうかが完済率を大きく変えるような値が存在する」という点が挙げられます。`LowDoc=S`、`City=TUSTIN`、`RevLineCr=T`のようなカテゴリ変数の特定の水準だけでなく、`Approval_Date=2007-05-10`、`SBA_Appv=5,100,000`のように、日付変数や数値変数であっても、特定の水準における完済率だけが平均から特に大きく乖離している、という例が見られました。

このため、`LowDoc`や`RevLineCr`のような明らかなカテゴリ変数だけではなく、日付変数や数値変数についても、それらが本来持っている順序や連続性のことを忘れてカテゴリ変数と同じように取り扱ってみることが、予測モデルの精度を向上させるうえで重要なポイントの一つであったと考えています。

このコンペに関するいくつかの解法記事で指摘されている通り、本コンペのデータに見られたこのような特徴は、それが人工的に生成された過程に由来するものだと思われます。この特徴のゆえに、ドメイン知識に即した特徴量生成が予測モデルの精度向上につながるケースはあまり多くなさそうだ、というのが私個人の見立てであり、チームの方針としても、「仮にすべてがカテゴリ変数からなるデータだとして、予測精度を上げるためにはどういう手法が良さそうかを考える」という方針に舵を切っていきました。

### 2-1. ベース特徴量 
予測モデルの学習に用いるベース特徴量として、`City`や`FranchiseCode`のような水準数の多いカテゴリ変数を含め、すべての変数をなるべくそのままの形で残しました。ただし、2種類の日付変数は、それぞれを年・年月・年月日という3段階のレベルに変換しています。また、EDAの結果を参考にしつつ、交互作用のありそうな変数の組（`State`と`BankState`、`RevLineCr`と`LowDoc`など）については、その文字列を単純に結合した特徴量を追加しました。

追加した特徴量に対して、予測モデルの精度を確認しながら変数選択を行い、最終的に26個の変数をベース特徴量として採用することにしました。

### 2-2. One-Hot Encoding 特徴量
どの変数にも完済率を予測するうえで無視できない値が存在する一方で、多くの変数のほとんどの値は、出現数が少ないとか、完済率が平均からそれほど乖離しないという意味で、それほど重要ではありませんでした。そこで、各ベース変数に含まれる重要水準を単独で取り出し、積極的に予測モデルに取り込むことができるように、One-Hotエンコーディングによる特徴量生成を行いました。これによって、決定木系のモデルの`colsample_bytree`や`colsample_bynode`などのパラメーターに関連した分岐生成において、特定の水準による分岐が採用される確率を引き上げることをイメージしています。

各変数から一部の水準だけを抽出する際、その重要度を判別する基準として「母集団比率に関するZ検定の検定統計量」を用いました。つまり、学習用データの中におけるある水準の出現数$n$、その水準のデータにおける`MIS_Status`の平均値$p$、データ全体における`MIS_Status`の平均値$\mu$を用いて、$$|z| = |{p-\mu}|/\sqrt{{\mu*(1-\mu)}/{n}} > u$$であるかどうかによって判定しました。たとえば`City=NASHVILLE`について、学習用データでの出現数は$n=599$、学習用データにおける完済率は$p=0.969942$、また、学習用データ全体での平均完済率は$\mu=0.892689$であることから$z=5.71339$となり、採用基準が$u=6.0$の場合は不採用、$u=5.5$の場合は採用してOne-Hotエンコーディングの対象にしています。

この方法による重要水準の抽出は、たとえば、以下のようなコードで実装することができます。

```python:important_levels.py
# train : 学習用データの polars.DataFrame
def important_levels(column, u=6.0, mu=0.8926891531):
    levels_df = (
        train
        .group_by(column)
        .agg(pl.count(), pl.mean('MIS_Status'))
        .rename({'count':'n_', 'MIS_Status':'p_'})
        .with_columns(
            ((pl.col('p_') - mu).abs() / np.sqrt(mu * (1 - mu) / pl.col('n_'))).alias('z_')
        )
    .filter(pl.col('z_').gt(u))
    )
    return(list(*levels_df.select(column)))

important_levels('City', u=5.5)
```

抽出の基準値は、最終的な予測モデルの学習に使用したデータでは5.5または6.0としており、ベース特徴量の全水準のうち、約330個または約280個の水準をOne-Hotエンコーディングによる2値の特徴量に変換しました。

# 3. モデル構築

### 3-1. 分類モデルの選択
前処理後のデータをもとに、➀LightGBM、➁CatBoost、➂XGBoost、➃LogisticRegressionなどの分類モデルの構築を試みました。一部の他の参加者の方が解法ブログで記載されているように、CatBoostが特に有効であったと感じています。

- **LGBMClassifier**
比較的学習が高速であることから、前処理の変更などの実験を回すうえで大変重宝しました。シングルモデルとしては、私自身はCatBoostClassifierを超える精度のモデルをLightGBMで作ることはできませんでしたが、チームメンバーの一人は、さまざまな前処理の組み合わせを探索する手法を実装し、LGBMClassifierを用いてチーム内のシングルモデルによる全投稿のうちで最も高い平均F値（Private LBで3位〜4位相当）を出す予測モデルを構築していました。なお、詳細に検証したわけではありませんが、LightGBMでは、ベース特徴量の一部をカテゴリー化せず、数値のまま扱うことで精度がわずかに改善したように思います。

- **CatBoostClassifier**
今回のコンペにおいて、CatBoostは非常に有効で、ベース特徴量をすべてカテゴリー変数として学習させるだけでも、PublicLBでメダル圏内に入る平均F値を記録することができました。なお、今回のコンペでは、CatBoostでも学習用データにOHE特徴量を加えることで予測精度がかなり向上しましたが、本来CatBoostでは前処理段階でOHEを行うことは非推奨とされていると思います。もしかすると、OHE特徴量を加えなくても、CatBoostに実装されているカテゴリー変数の処理手法を適切に用いることで、同程度以上の精度向上を狙うことができたのかもしれません。

- **XGBClassifier（採用せず）**
XGBoostでも`enable_categorical=True`とすることでカテゴリー変数を扱うことができます。しかし、私個人の感触としては、LightGBMやCatBoostに比べて、今回のコンペでXGBoostの精度を十分に向上させることは難しかったように思います。
- **LogisticRegression（採用せず）**
ロジスティック回帰では、カテゴリー変数をそのまま入力とすることができません。そこで、OHE特徴量の数を1,000個～2,000個程度に増やすことで、なるべく学習用データの特徴を幅広く捉えるように工夫しました。デフォルトの設定では計算が収束しませんでしたが、`solver="newton-cholesky"`とし、`max_iter`を十分大きく設定することで安定的に計算が収束するようになりました。また、係数の安定のためにリッジ正則化項を加えました。CV上はLightGBMを超えるAUCや平均F値を記録した割に、なぜか、LB上では評価値が思うように伸びなかったため、最終的に採用を見送りました。

### 3-2. スタッキング
LightGBMとCatBoostによる確率予測をロジット変換した値を入力として、ロジスティック回帰によるスタッキングを行いました。これは、二つの予測モデルの確率予測$p_{lgb}$および$p_{cat}$をもとに、以下の算式で新たな確率予測を求めることに相当します。$$p={p_{lgb}^{0.37} * p_{cat}^{0.63}}/(p_{lgb}^{0.37} * p_{cat}^{0.63} + (1-p_{lgb})^{0.37} * (1-p_{cat})^{0.63})$$

このスタッキングモデルは、しきい値の設定に応じて、たとえば次の右側のプロットのような決定境界を定めます。結果を見る限り、確率予測の単純な加重平均アンサンブルを採用して、そのウェイトを調整するのと大差ないかもしれません。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1484754/92405f59-997e-c210-8fb9-ec9b4b7c649f.png)

なお、当初はロジット変換前後の確率予測を併用することも検討していました。しかし、ロジット変換前後の確率予測を入力に含めてoptunaでチューニングしたところ、`penalty="l1"`（ラッソ正則化）が選択され、ロジット変換前の確率予測の係数がゼロとなったため、ロジット変換後の確率予測だけを使うことにしました。

### 3-3. その他
**パラメーターチューニング**
各種パラメーターのチューニングには optuna を利用し、プレーンな3-Fold CVにおけるROC-AUCの結果を最大化するように探索を行いました。理論平均F値（しきい値を変えたときの平均F値の最大値）を最適化する方法も試しましたが、理論平均F値が正解ラベルの偶然的な並びによってかなり気まぐれに変動することを踏まえて、AUCを最適化する方が予測モデルの信頼度を高められると判断しました。

**バリデーション方針**
チューニング後のパラメーターを用いてプレーンな100-Fold CVを行い、各Foldで学習した予測モデルのOOFに対する確率予測をもとに、理論平均F値、ROC-AUC、Log-Lossを評価しました。また、しきい値を変えたときの平均F値と確率予測に基づく疑似的な期待平均F値（後述）の推移を可視化し、確率予測が極端に歪んでいないことを確認するようにしました。

**不均衡データへの対応**
今回のデータは、学習用データにおける正例と負例の割合がおよそ9対1であり、かなり不均衡な分布になっていました。しかし、私たちのチームでは、アンダー/オーバーサンプリングや、予測モデルの学習時における`sample_weight`の調整などの手法を採用していません。これらの手法が平均F値を向上させる方向ではあまり機能せず、また、予測ラベルを分けるしきい値の決定において、できるだけ歪みのない確率予測を利用する必要があったためです。

# 4. 後処理

今回のコンペの評価指標が「平均F値」だったことで、私たちは非常に頭を悩ませました。もし正確な確率予測を得ていたとしても、しきい値の選択次第で、平均F値はかなりの幅で変動します。実際、今回のコンペでは、最終提出でのしきい値の選択を誤るだけで、大きく順位を落とす結果になっていたはずです。

### 4-1. 予測ラベルの決定方法
私たちのチームでは、確率予測を予測ラベルに変換する際の判断基準として、以下のような方法を併用していました。
1. **最良しきい値法**
OOF（またはPublic LB）における平均F値を最大化する確率予測ベースのしきい値を採用する方法です。しきい値を変えていくとき、平均F値は正解ラベルの並びによって局所的にはギザギザに増減するため、この方法で導かれるしきい値は必ずしも安定しません。
2. **最良比率法**
OOFにおける平均F値を最大化する予測ラベルの比率を採用する方法です。学習用データと評価用データの正解ラベルの割合に乖離があるときは信頼できず、また、最良しきい値法と同じ理由から、最良の比率は必ずしも安定しません。
3. **モンテカルロ法**
評価用データに対する確率予測を用いて疑似的な正解ラベルを乱数的に繰り返し生成し、各試行において平均F値を最大化するしきい値を求め、その平均値や中央値を採用する方法です。確率予測が歪んでいる場合は信頼できません。また、試行回数を増やすほど計算時間が増加します。
4. **期待混同行列に基づく疑似平均F値法**
評価用データに対する確率予測を用いて混同行列の4要素の期待値を求め、この混合行列に基づく疑似的な平均F値を最大化するしきい値を採用する方法です。確率予測が歪んでいる場合は信頼できません。

なお、予測に含める陽性数と平均F値ないし疑似平均F値との関係は、以下のようなコードで計算し、可視化することができます。
```python:macro averaged f1
# label : データの第1行～第N行に対応する正解ラベルのリスト
# proba : データの第1行～第N行に対応する予測確率のリスト

import polars as pl
import matplotlib.pyplot as plt
import seaborn as sns

df = (
    pl.DataFrame({'label':label, 'proba':proba})
    .sort(by='proba', descending=True)
)

# 混同行列の4要素から平均F値を求める関数
def macro_f1(tp, fp, fn, tn):
    return tp / (2 * tp + fp + fn) + tn / (2 * tn + fp + fn)

# 混同行列の4要素を、すべての予測ラベルを0とする状態で初期化
TP, FP, FN, TN = 0, 0, df['label'].sum(), len(label) - df['label'].sum()
scores = list()
score  = macro_f1(TP, FP, FN, TN)
scores.append(score)

# 混同行列の4要素の期待値を、すべての予測ラベルを0とする状態で初期化
tp, fp, fn, tn = 0, 0, df['proba'].sum(), len(label) - df['proba'].sum()
pseudo_scores = list()
pseudo_score  = macro_f1(tp, fp, fn, tn)
pseudo_scores.append(pseudo_score)

for i in range(len(df)):
    label_, proba_ = df.row(i)
    if label_ == 1:
        TP += 1
        FN -= 1
    else:
        FP += 1
        TN -= 1
    score = macro_f1(TP, FP, FN, TN)
    scores.append(score)
    tp += proba_
    fn -= proba_
    fp += 1 - proba_
    tn -= 1 - proba_
    pseudo_score = macro_f1(tp, fp, fn, tn)
    pseudo_scores.append(pseudo_score)

sns.lineplot(x=list(range(len(df) + 1)), y=scores)
sns.lineplot(x=list(range(len(df) + 1)), y=pseudo_scores, linestyle='--')
```

また、投稿数に余裕がある範囲で、同じ確率予測に基づいて予測ラベルの陽性数を変えたものを投稿し、Public LB上での平均F値カーブの形状を確認するようにしました。

### 4-2. 最終提出

二つの最終提出を決めるために、私たちのチームでは「もし手元の確率予測が真だと仮定した場合に、平均F値、およびそれを最大化するしきい値は、どの程度、どのような理由で変動するか」を数値実験によって検討しました。詳細は別の記事に譲ることとし、ここでは、その結果から読み取った2つのポイントを紹介します。
- 確率予測が真の確率に近いほど、平均F値のカーブは高い水準に位置しやすい
- 確率予測が真であっても、平均F値を最大化するしきい値は正解ラベルの実現値に応じて大きく変動する（この変動幅はデータ数が少ないほど大きい）

**数値実験における試行の例** （真の確率予測に基づく平均F値の変動）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1484754/4e2f56f8-3e1c-74cc-0b64-e342b39f99af.png)

この考察を踏まえ、私たちのチームでは、上に掲げたいずれの決定方法も使用せず、「Local CVとPublic LBでの評価値が最良の確率予測に基づき、大小異なるしきい値によって2値の予測ラベルに変換したものを最終的な二つの予測とする」という方針を取ることとしました。

この方針のもとで、予測の陽性数を38,500個、38,800個とした二つの予測を暫定的に選択していたところ、陽性数を38,500個とした方の予測が、Public LBだけでなくPrivate LBでも1位となる平均F値を記録していました。最終選択の方針は、どちらかといえば優勝よりもshake down対策を主眼に置いたもので、優勝という結果はまさに「運が良かった」と言うほかありません。

# 5. 結び

私たちのチームの所属会社では、アクチュアリーを中心とするメンバーで、定期的にデータサイエンスに関する勉強会を行っています。今回のチームは、この勉強会の参加者に声をかけて結成したものです。業務が多忙な中でも最後までこのコンペに取り組んでくださったチームメンバーに、何よりもまず感謝したいと思います。

後日予定されている解法共有会では、上記の概要とともに、少し別の視点から、チーム内でPublic LBスコアが最良だったシングルモデルの概要や、コンペ中の投稿管理などの取り組みに関するTipsなどもご共有できたらと思っています。

この記事を通して何か一つでもご参考になる点があればとても嬉しいです。最後までお読みいただき、ありがとうございました。