Kaggle Benzに挑戦

# はじめに
Kaggleのコンペに初挑戦した前回[こちらの記事](https://qiita.com/gucchi_JPN/private/25cdcddb73a604cbeab2)から1ヶ月ほど経ち、新たな課題に挑戦しました。今回は募集終了しているコンペであるものの、多くの参加者がいるコンペで腕試しには最適と思い、取り組みました。
これからPythonや機械学習を始めようと思っている方、学習をすでに始めており、力試しをしてみたいと思っている方にわかりやすい記事を目指します！
ご意見やご質問頂けたら嬉しいです。

# 挑戦した課題
- [Kaggle Benz](https://www.kaggle.com/competitions/mercedes-benz-greener-manufacturing)
### 内容
- メルセデスベンツは堅牢な車両評価手法を構築してきたものの、そのテストにかかる時間の最適化には多くの時間を要している
- 多くの特徴量から[次元の呪い](https://qiita.com/tn1031/items/96e7131cd41df1aebe85)に対処しつつ、テストに合格するまでの時間を予測する([回帰問題](https://qiita.com/maskot1977/items/c648e398b3b86cb68ae5))

# 実行環境
- プロセッサ：12th Gen Intel(R) Core(TM) i5-1240P
- 開発環境：Google Colaboratory
- 言語：Python
- 使用したライブラリ：Pandas, Numpy, Matplotlib, scikit-learn

# 分析するデータ
- データ数：trainデータ → 4209, testデータ → 4209
- 特徴量の数 : 376 (カテゴリ : 8, バイナリ : 368)
※ただし各特徴量がなにをあらわしているかの説明はなし
- 欠損値 : なし


# 分析の流れ
1. データの読み込み
1. データのチェック
1. データの前処理・特徴量エンジニアリング

- 今回は試行錯誤しながら行ったため特徴量エンジニアリングの途中でスコアの確認が入っています。ご了承ください

### 1. データの読み込み
#### 実行したコード
```python: kaggle_benz.ipynb
# 必要なモジュールのインストール
!pip install pycaret[full]

# 必要なモジュールのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
from sklearn.model_selection import train_test_split

warnings.filterwarnings("ignore")
```
- 各モジュールのURLはこちら
[pandas](https://pandas.pydata.org/), [numpy](https://numpy.org/ja/), [matplolib](https://matplotlib.org/), [seaborn](https://seaborn.pydata.org/), [PyCaret](https://pycaret.org/)

```python: kaggle_benz.ipynb
# データの読み込み
train = pd.read_csv('/content/train.csv')
test = pd.read_csv('/content/test.csv')

# 学習データの最初の5行を表示
train.head(5)

# 学習データの最初の5行を表示
test.head(5)
```
![train_head.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3838907/7ef33894-5df1-bdb5-eeb7-511569c5becf.png)

- 特徴量とてもおおいので一部を記載

```python: kaggle_benz.ipynb
# Id列削除前の型を確認
print("The train data size before dropping Id feature is : {} ".format(train.shape))
print("The test data size before dropping Id feature is : {} ".format(test.shape))

# Idだけ別に保持する
train_ID = train['ID']
# test_ID = test_data['ID']
test_ID = test['ID']

# 予測に不要なのでId列を削除する
train.drop("ID", axis = 1, inplace = True)
# test_data.drop("ID", axis = 1, inplace = True)
test.drop("ID", axis = 1, inplace = True)

# Id列削除後の型を確認
print("\nThe train data size after dropping Id feature is : {} ".format(train.shape))
print("The test data size after dropping Id feature is : {} ".format(test.shape))
```
- 学習に不要なId列が削除されていることを確認

```python: kaggle_benz.ipynb
# データの特徴を大まかにみる
train.describe()

#データ型などを確認
train.info()
```
- X0~X8はカテゴリデータ
- X10~X385はバイナリのデータ

 ## 2. データのチェック

 #### 実行したコード
 **ターゲット変数のチェック**
 ```python: kaggle_benz.ipynb
 # ターゲット変数の外れ値をクリップする
train["y"] = train['y'].clip(upper=155)
y = train['y']

# 特徴量とターゲット変数に分ける
X = train.drop(columns=['y'])
```
![y散布図.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3838907/a48a7e5a-187c-207a-3bf4-686a907a6b46.png)
- ターゲット変数yに外れ値が無いか確認
- 250を超える値が１つだけあるのでクリップする

**欠損値のチェック**
```python: kaggle_benz.ipynb
#train, testデータの各特徴量に欠損値がいくつあるかを確認

missing_values = train.isnull().sum()
print(missing_values)

missing_values = test.isnull().sum()
print(missing_values)
```

- 今回は欠損値ないのでそのまま

**特徴量のチェック**
- カテゴリデータ (X0, X1, X2, X3, X4, X5, X6, X8 の８つ)

```python: kaggle_benz.ipynb
# 各カテゴリデータのターゲット変数yとの相関を確認
import matplotlib.pyplot as plt

def plot_bar_chart(data, features):
    """
    指定された複数の特徴量の値ごとの割合を2×4のマトリクス状に棒グラフで表示する関数

    Args:
        data: データフレーム
        features: 特徴量の名称のリスト
    """
    fig, axes = plt.subplots(2, 4, figsize=(20, 12))
    axes = axes.flatten()

    for i, feature in enumerate(features):
        values = data[feature].value_counts()
        labels = values.index.to_list()
        percentages = (values / len(data)) * 100

        axes[i].bar(labels, percentages)
        axes[i].set_xlabel(feature)
        axes[i].set_ylabel("Percentage (%)")
        axes[i].set_title(f"Value Counts of {feature}")

    # 不要なサブプロットを非表示にする
    for j in range(len(features), len(axes)):
        fig.delaxes(axes[j])

    plt.tight_layout()
    plt.show()

features = ["X0", "X1", "X2", "X3", "X4", "X5", "X6", "X8"]
plot_bar_chart(train, features)
```
- X4以外はyと相関を持っていそうだとわかる
![カテゴリデータとyの相関.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3838907/a8c77df9-8b4f-b903-e878-142e520b21c0.png)


- バイナリデータ

```python: kaggle_benz.ipynb
# バイナリデータを抽出
train_binary = train_exp.drop(['X0', 'X1', 'X2', 'X3', 'X4', 'X5', 'X6', 'X8'], axis=1)
```

- yと相関のある特徴量があるのか調べる(相関係数0.2できってみる)
```python: kaggle_benz.ipynb
train_cat_corr = train_binary

# 相関係数行列を計算
corr_matrix = train_cat_corr.corr()

# ヒートマップを表示
plt.figure(figsize=(12, 10))
sns.heatmap(corr_matrix, cmap='coolwarm')
plt.title("Correlation Heatmap")
plt.show()

```

- 特徴量どうしで相関の高いものがある → 多重共線性の可能性
- 特徴量の数自体がとても多い → そのままだと次元の呪いの影響がありそう
![バイナリデータのヒートマップ.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3838907/a1100a98-3657-52cf-2090-1236022b7ded.png)

 ## 3. データの前処理・特徴量エンジニアリング
 
#### カテゴリデータ
- X4以外のデータに対してターゲットエンコーディングを実施する
- 例としてX0に対する処理を記載(ほかのカテゴリデータも同様)
#### 実行したコード

- X0について各カテゴリごとにyの平均値を算出、平均値が大きい順にカテゴリをソートしたものをX0_y_sortedに格納
```python: kaggle_benz.ipynb

X0_y_mean = train.groupby('X0')['y'].mean().sort_values(ascending=False)
X0_y_sorted = X0_y_mean.index.to_list()
# print(X0_y_sorted)
```

- 各カテゴリで得た平均値を正規化し、train, testのカテゴリの値とする
```python: kaggle_benz.ipynb

from sklearn.preprocessing import MinMaxScaler

scaler = MinMaxScaler()
X0_y_mean_normalized = scaler.fit_transform(X0_y_mean.values.reshape(-1, 1))

for i, category in enumerate(X0_y_sorted):
    train.loc[train['X0'] == category, 'X0'] = X0_y_mean_normalized[i][0]
    test.loc[test['X0'] == category, 'X0'] = X0_y_mean_normalized[i][0]
```

- test["X0"]にのみ含まれるカテゴリはtrain["X0"]の平均値で埋める

```python: kaggle_benz.ipynb

# Create sets of unique values for X0 in train and test data
train_X0_unique = set(train["X0"].unique())
test_X0_unique = set(test["X0"].unique())

test_only_categories = test_X0_unique - train_X0_unique

test.loc[test["X0"].isin(test_only_categories), "X0"] = np.nan

train_X0_mean = train["X0"].mean()

test["X0"].fillna(train_X0_mean, inplace=True)
```

#### カテゴリデータのみで予測してみる
- R2スコア : 0.54125 (Privateスコア)
- (この時点の)順位 : 2670 / 3775位
- やはり大量のバイナリデータをうまく扱う必要がありそう

#### バイナリデータへのアプローチ①
- バイナリデータについてtrain, testを結合
``` python: kaggle_benz.ipynb
all_binary = pd.concat([train_binary, test_binary], axis=0)
all_binary.shape
```

- 次元削減として一般的なPCAを実施
- tranデータに対して実施(fit_transform)したPCAをtestデータにも適用する
```python: kaggle_benz.ipynb
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
all_binary_scaled = scaler.fit_transform(all_binary.drop(columns=["y"]))
print(all_binary_scaled.shape)

#PCAの適用
pca = PCA(n_components=15)
train_binary_pca = pca.fit_transform(all_binary_scaled[:train_exp.shape[0]])
test_binary_pca = pca.transform(all_binary_scaled[train_exp.shape[0]:])

# PCA 結果をデータフレームに変換
pca_columns = [f"PC{i+1}" for i in range(train_binary_pca.shape[1])]
train_binary_pca_df = pd.DataFrame(train_binary_pca, columns=pca_columns)
test_binary_pca_df = pd.DataFrame(test_binary_pca, columns=pca_columns)

#結果の表示
train_binary_pca_df.head(10)
```

![PCA寄与率.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3838907/c9c9e47e-60d0-a96b-66b2-b04db5bd6ffd.png)
- 各主成分の寄与率はそこまで高くないことがわかる

#### PCAで得た主成分とカテゴリデータでアンサンブルして予測してみる
- 今回はカテゴリデータとの相関係数が高くない第5主成分までを採用

![カテゴリ_PCA_相関係数.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3838907/02cc6d1c-58ef-f8ae-c577-493631b6203e.png)

- カテゴリ + PCAのデータを作成する

```python: kaggle_benz.ipynb
## train, testで必要な主成分を取り出す
train_pca_selected = train_binary_pca_df.iloc[:, 0:5]
test_pca_selected = test_binary_pca_df.iloc[:, 0:5]

X = pd.concat([train[['X0', 'X1', 'X2', 'X3', 'X5', 'X6', 'X8']].reset_index(drop=True), train_pca_selected], axis=1)
X_test = pd.concat([test[['X0', 'X1', 'X2', 'X3', 'X5', 'X6', 'X8']].reset_index(drop=True), test_pca_selected], axis=1)
y = train['y'].reset_index(drop=True)
```
- カテゴリ + PCAのデータに対してLightGBMを使用
```python: kaggle_benz.ipynb
best_params = {
    'objective': 'regression',
    'metric': 'rmse',
    'verbosity': -1,
    'boosting_type': 'gbdt',
    'learning_rate': 0.06400413604943948,
    'num_leaves': 43,
    'min_child_samples': 100,
    'min_child_weight': 0.0027616347599050636,
    'subsample': 0.9951009359507461,
    'colsample_bytree': 0.9896019466335962,
    'reg_alpha': 0.039327097602159196,
    'reg_lambda': 3.9344354817602074,
    'max_depth': 6,
    'n_estimators': 660
}

# モデルの訓練
best_model = lgb.LGBMRegressor(**best_params)
best_model.fit(X, y)
```

- テストデータに対して予測

```python: kaggle_benz.ipynb
# X_testデータに対する予測
y_pred = final_model.predict(X_test)

# 予測結果の表示
# print("Predictions on X_test:")
# print(y_pred)
```

#### カテゴリ + PCAに対してRidge回帰を使用
- 次元削減によって13変数まで減っているので影響の小さい特徴量を消してしまうLassoよりも適当と判断

```python: kaggle_benz.ipynb
from sklearn.linear_model import LinearRegression, Lasso, Ridge
from sklearn.metrics import r2_score

model = Ridge(alpha=0.1)

# Fit the model to the training data
model.fit(X, y)
print(r2_score(model.predict(X), y))
y_pred_reg = model.predict(X_test)
```

#### LgihtGBMとRidge回帰の結果をアンサンブルし予測
```python: kaggle_benz.ipynb
y_final = 0.22 * y_pred + 0.78 * y_pred_reg
```
- R2スコア : 0.54155 (Privateスコア)
- (この時点の)順位 : 2641 / 3775位
- 大幅な改善は見られず → 他のアプローチが必要そう

#### バイナリデータへのアプローチ②
戦略
- バイナリデータから任意に10個取り出す
- スコアがそれまでよりも良いかどうかで10個の特徴量の重みを更新
- 重みが閾値を超えれば、その特徴量が重要だと判定し抽出する

#### 実行したコード

- 重みの初期値や閾値の準備
```python: kaggle_benz.ipynb
from numpy.random import choice
n_elements = train.shape[0]
feats = train.columns.values[9:]
n_feats = len(feats)
weights = np.array(np.zeros(n_feats))
weights += 1/n_feats

best_feats = []
best_score = 0
epochs = 3000 # number of rounds / feature bags
n_features = 10 # try out different values
wt_g = .2 #weight growth rate
w_threshold = .2 #weight to add feature
# kf = KFold(n_elements, n_folds=3, shuffle=True)
# kf = KFold(3, shuffle=True).get_n_splits(train.values)
kf = KFold(n_splits=3, shuffle=True)
y = train['y'].values
scores = []
top_feats = np.array([])
```

- 重要な特徴量の選別

```python: kaggle_benz.ipynb
clipped_y = y
for i, ind in enumerate(list(range(epochs))):
      sample_feat_ids = choice(a=n_feats, size=n_features,
                              replace=False, p=weights)
      sample_feats = np.append(top_feats, feats[sample_feat_ids])
      tst_P = np.array(np.zeros(n_elements))

      X = train.loc[:,sample_feats]
      for trn, tst in kf.split(X, clipped_y):
          # X = train.loc[:,sample_feats]
          trn_X, tst_X = X.iloc[trn,:], X.iloc[tst,:]
          trn_Y, tst_Y = clipped_y[trn], clipped_y[tst]
          #mod = SVR(C=50, epsilon=3, gamma=.2)
          #mod = GBR(alpha=.01, n_estimators=50, max_depth=5, min_samples_leaf=15, subsample=.5, random_state=1776)
          mod = RFR(n_estimators=100, max_depth=12, max_features=20, min_samples_leaf=4, n_jobs=2, random_state=1776)
          mod.fit(trn_X, trn_Y)
          tst_P[tst] = mod.predict(tst_X)

      # I don't want to overfit to outliers so I am clipping all y's at 155
      tst_rsq = r2(y_pred=tst_P, y_true=clipped_y)

      if ind > 29:
          scores.append(tst_rsq)
          ma_rsq = np.mean(scores[-30:])
          if ind % 25 == 0:
              print(ind, ma_rsq, tst_rsq)
          if tst_rsq > ma_rsq:
              weights[sample_feat_ids] *= (1+wt_g)
              sum_w = np.sum(weights)
              weights /= sum_w
              if tst_rsq > best_score:
                  best_score = tst_rsq
                  best_feats = sample_feats
          else:
              weights[sample_feat_ids] *= (1-wt_g)
              sum_w = np.sum(weights)
              weights /= sum_w
          mx_w = np.max(weights)

          # add feature to the top feats if weight is > threshold
          if mx_w > w_threshold:
              feat_imps = pd.Series(index = feats, data = weights)
              new_feats = feat_imps[feat_imps > w_threshold].index.values
              top_weights = feat_imps[feat_imps > w_threshold].values
              top_feats = np.append(top_feats, new_feats)
              feats = list(feats)
              weights = list(weights)
              for f,w in zip(new_feats,top_weights):
                  print(f, w)
                  feats.remove(f)
                  weights.remove(w)
              n_feats = len(weights)
              feats = np.array(feats)
              weights = np.array(weights)
              sum_w = np.sum(weights)
              weights /= sum_w

          if ind % 25 == 0:
              print(mx_w)
      else:
          scores.append(tst_rsq)
          print(ind, tst_rsq)

print(best_feats)
print(best_score)
```

![選ばれた特徴量.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3838907/7b8abe95-88b0-5236-2825-3fe495ff8257.png)

- 選ばれた特徴量を用いて予測
- 選別したデータをtrain_best, test_bestとしている

```python: kaggle_benz.ipynb
mod = RFR(n_estimators=100, max_depth=12, max_features=20, min_samples_leaf=4, n_jobs=2, random_state=1776)
mod.fit(train_best, clipped_y)

pred_rfr = mod.predict(test_best)
```

- 同様にLightGBMでも再び予測しサンサンブルで最終的な予測値とする
```python: kaggle_benz.ipynb
y_final = 0.9 * pred_lgbm + 0.1 * pred_rfr
```

#### 最終的なスコア・順位
- R2スコア : 0.54946 (Privateスコア)
- (この時点の)順位 : 1228 / 3775位
- スコアは大幅に改善

# さいごに
課題を通じて感じたこと
- データの特性によっては線形回帰などの古典的なモデルでもそこそこよい精度が出せることもある
- 次元削減 = PCAというイメージを持っていたが、今回のように各主成分の寄与率が低く、効果が小さいこともある
- 最後に実施した特徴量の抽出のような地道な手法も時には必要

Kaggleのようなコンペでよい成績を残すには教科書にかかれているような手法の理解はもちろん、さまざまな手法・アイデアをトライアルアンドエラーしていくことが大切だとわかりました。まだまだ駆け出しですが、これからも頑張りたいと思います。