0. はじめに

Kaggleの有名なコンペティション、House Pricesの挑戦記録です。
この記事のモチベーションについては挑戦記の第一部で触れているので、そちらをご覧ください。

とにかくあらゆる作業とその理由・思考を言語化しようというのがこの記事での目論見で、「Kaggleに勝つ」「ハイスコアを目指す」ことを主眼に置いてはいません。一つひとつの作業や試行を真の意味で理解するために、このコンペに関する他の記事はいっさい参照していません。

前回は、ランダムフォレスト回帰モデルを作成しました。ベースラインからRMSEが改善し0.125程度になりました。今回は、XGBoostでモデルを作成し、提出用データ作成から提出・評価まで一気に進もうと思います。

1. XGBoost実装

XGBoostを選択するメリットとしては、

ブースティングによる予測精度の向上: XGBoostは、複数の弱い学習器（通常は決定木）を逐次的に学習させるブースティングという手法を採用。各学習器は前の学習器の誤差を修正するように学習するため、アンサンブル全体として高い予測精度を実現。一方、ランダムフォレストは、複数の決定木を独立に学習させ、その予測結果を平均化または多数決で決定するバギングという手法を用いる。
過学習の抑制: XGBoostは、L1およびL2正則化という過学習を防ぐための仕組みを内蔵。これにより、モデルが複雑になりすぎるのを抑制し、未知のデータに対する汎化性能を高める。

一方でデメリットとしては、

解釈性の低下: ランダムフォレストの方が、個々の決定木が比較的独立しているため、モデル全体の挙動をある程度解釈しやすい。
複雑なパラメータチューニング: 調整すべきハイパーパラメータがランダムフォレストよりも多く、良い性能を得るのに複雑なチューニングが必要。

1.1 データ取得〜テストモデル実装

上記を踏まえて、実装に進みます。データは、ランダムフォレストで作成したものを用います。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os

from constants import *

plt.style.use('ggplot')
sns.set_palette('winter')

# 前回作成したデータを呼び出し
data = pd.read_csv(os.path.join(RAW_PATH, '0200_data.csv'), index_col=0)

まずは、ハイパーパラメータをいっさいいじらないプレーンなXGBモデルを作ります。

from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import root_mean_squared_error, r2_score

X = data.drop('log_price', axis=1)
y = data['log_price']

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2)

test = XGBRegressor(n_jobs=-1, random_state=86)

test.fit(X_train, y_train)
y_pred = test.predict(X_val)

rmse = root_mean_squared_error(y_val, y_pred)
r2 = r2_score(y_val, y_pred)

print('Test: XGBRegressor')
print(f'  RMSE: {round(rmse, 4)}')
print(f'  R2: {round(r2, 4)}')

出力
Test: XGBRegressor
RMSE: 0.1393
R2: 0.8847

ランダムフォレストより若干スコアが悪化していますね。

1.2 SHAP分析

テストモデルにおいて、どの特徴量がモデルの出力に影響を与えているかを可視化します。今回は、SHAPを利用します。SHAPについての詳しい説明は他の記事を参考にしてください。

import shap

shap.initjs()

explainer = shap.TreeExplainer(model=test)
shap_values1 = explainer.shap_values(X)
shap.summary_plot(shap_values1, X)

やはり、床面積と全体的な品質、そして立地がかなり強い影響を与えているようです。おおむね、ランダムフォレストと同様の傾向がありそうですね。

1.3 Optunaハイパーパラメータチューニング

ランダムフォレストと同様に、Optunaを用いてハイパーパラメータをチューニングしていきます。

import optuna

def objective(trial):
    
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2)
    
    params = {
        # 学習率
        'learning_rate': trial.suggest_float('learning_rate', 0.001, 1, log=True),
        # 最大深度
        'max_depth': trial.suggest_int('max_depth', 3, 30),
        # 決定木の葉の重みの下限
        'min_child_weight': trial.suggest_int('min_child_weight', 1, 10),
        # 各決定木においてランダムに抽出される標本の割合
        'subsample': trial.suggest_float('subsample', 0.1, 1),
        # 各決定木においてランダムに抽出される特徴量の割合
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.1, 1),
        # 決定木の葉の追加による損失減少の下限
        'gamma': trial.suggest_float('gamma', 0.1, 10),
        }
    
    model = XGBRegressor(
        **params,
        random_state=86,
        n_jobs=-1
        )
    
    model.fit(X_train, y_train)
    y_pred = model.predict(X_val)
    
    rmse = root_mean_squared_error(y_val, y_pred)
    
    return rmse

# 学習を実行
study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=300, n_jobs=-1)

# 最適なパラメータとスコアを取得
study.best_params, study.best_value

出力
({'learning_rate': 0.06258673844191955,
'max_depth': 30,
'min_child_weight': 6,
'subsample': 0.575309470827,
'colsample_bytree': 0.7570165472459157,
'gamma': 0.2323178887112457},
0.11674569213833685)

スコアはだいぶ良くなりましたね。パラメータの重要性を可視化してみましょう。

optuna.visualization.plot_param_importances(study)

gammaが最も強い影響を持っているようです。ガンマは過学習を防ぐための正則化項として機能するパラメータで、具体的には、あるノードを分割した結果、損失関数がgammaよりも小さくしか減少しない場合、その分割は行われないように抑制するものです。このパラメータが最も重要なのは、少し意外でした。

それでは、このパラメータを用いて学習をしていきましょう。

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=88)

best_params = study.best_params

opt_xgb = XGBRegressor(
    **best_params,
    random_state=88,
    n_jobs=-1
    )

opt_xgb.fit(X_train, y_train)
y_pred = opt_xgb.predict(X_val)

rmse = root_mean_squared_error(y_val, y_pred)
r2 = r2_score(y_val, y_pred)

print('Optimized XGB')    
print(f'  RMSE: {round(rmse, 4)}')
print(f'  R2: {round(r2, 4)}')

出力
Optimized XGB
RMSE: 0.138
R2: 0.8763

1.4 正則化パラメータのチューニング

次に、過学習を抑制するために、reg_alpha（L1正則化）とreg_lambda（L2正則化）を調整していきます。こちらは、学習曲線を見ながら調整していこうと思います。

L1正則化

from utils import plot_learning_curve

alphas = [0.001, 0.005, 0.01, 0.05, 0.1, 1]

plt.figure(figsize=(16, 15))

for i, alpha in enumerate(alphas):
    plt.subplot(3, 2, i+1)
    
    plot_learning_curve(
        XGBRegressor(**best_params, reg_alpha=alpha, random_state=88, n_jobs=-1),
        X_train, y_train,
        title=f"Learning Curve: L1 Regularized XGB Regr. (alpha={alpha})"
        )
    plt.tight_layout()
    
plt.show()

図を見ると、0.1まではあまり変化は見られず、　0.1から1で正則化の効果が見られます。
0.1から1の間で再度探索します。

alphas = [0.1, 0.2, 0.4, 0.6, 0.8, 1]

plt.figure(figsize=(16, 15))

for i, alpha in enumerate(alphas):
    plt.subplot(3, 2, i+1)
    
    plot_learning_curve(
        XGBRegressor(**best_params, reg_alpha=alpha, random_state=88, n_jobs=-1),
        X_train, y_train,
        title=f"Learning Curve: L1 Regularized XGB Regr. (alpha={alpha})"
        )
    plt.tight_layout()
plt.show()

0.6が、スコアと過学習抑制のバランスが良さそうなので、この値を採用します。

L2正則化

次に、L2正則化パラメータを探索します。

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=86)

alpha = 0.6
lambdas = [0.01, 0.1, 1, 5, 10, 100]

plt.figure(figsize=(18, 15))

for i, lamb in enumerate(lambdas):

    plt.subplot(3, 2, i+1)
    
    plot_learning_curve(
        XGBRegressor(**best_params, reg_alpha=alpha, reg_lambda=lamb, random_state=14, n_jobs=-1),
        X_train, y_train,
        title=f"Learning Curve: L2 Regularized XGB Regr. (lambda={lamb})"
    )
    plt.tight_layout()
    
plt.show()

1から5の間で、抑制の効果がよく現れていて、10以上になると精度が落ちているのがわかります。ふたたび1~5で探索します。

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=86)

alpha = 0.6
lambdas = [1, 2, 3, 4, 5]

plt.figure(figsize=(18, 15))

for i, lamb in enumerate(lambdas):

    plt.subplot(3, 2, i+1)
    
    plot_learning_curve(
        XGBRegressor(**best_params, reg_alpha=alpha, reg_lambda=lamb, random_state=14, n_jobs=-1),
        X_train, y_train,
        title=f"Learning Curve: L2 Regularized XGB Regr. (lambda={lamb})"
    )
    plt.tight_layout()
    
plt.show()

3から4でスコアの若干の悪化が見られるので、lambdaは3とすることにします。
正則化パラメータを加え、最適化されたXGBで再度学習をします。

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=85)

lmbd = 3

regl_xgb = XGBRegressor(
    **best_params,
    reg_alpha=alpha,
    reg_lambda=lmbd,
    random_state=15,
    n_jobs=-1
)

regl_xgb.fit(X_train, y_train)

y_pred = regl_xgb.predict(X_val)

rmse = root_mean_squared_error(y_val, y_pred)
r2 = r2_score(y_val, y_pred)

print('Regularized XGB')    
print(f'  RMSE: {round(rmse, 4)}')
print(f'  R2: {round(r2, 4)}')

出力
Regularized XGB
RMSE: 0.1367
R2: 0.8912

ランダムフォレストより若干RMSEが下がっており、また正則化による過学習抑制の効果もあるため、これを提出用モデルとして採用します。

パラメータ、モデルをそれぞれ保存します。

# 最適化パラメータに正則化パラメータを追加
best_params.update([('reg_alpha', alpha), ('reg_lambda', lmbd)])

# メトリクスDFを作成
xgb_metrics = {
    'model': 'XGB Regr.',
    'params': f'{best_params}',
    'RMSE': rmse,
    'R2': r2
}

xgb_met_df = pd.DataFrame(xgb_metrics, index=[0])

# 他のメトリクスを呼び出して合体
metrics = pd.read_csv(os.path.join(METRICS_PATH, '0200_metrics.csv'), index_col=0)

metrics = pd.concat([metrics, pd.DataFrame(xgb_metrics, index=[0])])

metrics.to_csv(os.path.join(METRICS_PATH, '0300_metrics.csv'))

# モデルの保存
import joblib

joblib.dump(regl_xgb, os.path.join(MODEL_PATH, '0300_xgb_regr.joblib'))

2. 提出用データ作成

それでは、学習させたXGBoostを用いて提出用データを作っていきましょう。
これまでは訓練用データのみを処理・使用してきたので、テスト用データと訓練用データを合わせたデータで学習・予測をさせていきます。

2.1 データ処理

これまでの処理を元に、特徴量を作っていきます。まずは読み込み。

0001_process.ipynb

# 提出用データ前処理用のノートブック

import numpy as np
import pandas as pd
import os

from constants import *

# テストデータと訓練データを読み込んで合体
test = pd.read_csv(os.path.join(RAW_PATH, 'test.csv'), index_col=0)
train = pd.read_csv(os.path.join(RAW_PATH, 'train.csv'), index_col=0)
all_data = pd.concat([train, test])
all_data.info()

# ダミー変数化するもの、そのまま使用するものを取得
cols = [
'MSZoning', 'OverallQual', 'OverallCond',
'ExterQual', 'ExterCond', 'BsmtQual', 'BsmtCond',
'KitchenQual', 'Utilities', 'PoolQC', 'Functional', 
'GarageCars', 'GarageQual','GarageCond', 'LotFrontage'
]

data = all_data[cols]

学習に利用した特徴量を作成。

all_data['log_price'] = np.log(all_data['SalePrice'])

data['MSSubClass'] = all_data['MSSubClass'].astype(str)

data['log_price'] = all_data['log_price']

data['total_sf'] = all_data['1stFlrSF'] + all_data['2ndFlrSF'].fillna(0) + all_data['TotalBsmtSF'].fillna(0)

data['log_total_sf'] = np.log(data['total_sf'])

data.drop(['total_sf'], axis=1, inplace=True)

# フラグも作成
data['has_2nd'] = (all_data['2ndFlrSF'] > 0).astype(int)

data['has_bsmt'] = (all_data['TotalBsmtSF'] > 0).astype(int)

data['has_garage'] = (all_data['GarageArea'] > 0).astype(int)

data['has_pool'] = (all_data['PoolArea'] > 0).astype(int)

data['has_central_air'] = (all_data['CentralAir'] == 'Y').astype(int)

data['has_fireplace'] = (all_data['Fireplaces'] > 0).astype(int)

data['house_age'] = all_data['YrSold'] - all_data['YearBuilt']

data['remod_age'] = all_data['YrSold'] - all_data['YearRemodAdd']

data['total_porch_sf'] = all_data['OpenPorchSF'] + all_data['EnclosedPorch'].fillna(0) + all_data['3SsnPorch'].fillna(0) + all_data['ScreenPorch'].fillna(0)

data['LotFrontage'].fillna(0.0, inplace=True)
data['Utilities'].fillna('None', inplace=True)
data['KitchenQual'].fillna('TA', inplace=True)
data['Functional'].fillna('None', inplace=True)
data['GarageCars'].fillna(0, inplace=True)

from sklearn.preprocessing import StandardScaler

std_scaler = StandardScaler()

train['log_price'] = np.log(train['SalePrice'])

neighbor_log_price = train.groupby('Neighborhood')['log_price'].mean().sort_values(ascending=False)

neighbor_log_price_stdsclaled = std_scaler.fit_transform(neighbor_log_price.values.reshape(-1, 1))

neighbor_log_price_dict = dict(zip(neighbor_log_price.index, neighbor_log_price_stdsclaled))

data['neighbor_weight'] = all_data['Neighborhood'].map(neighbor_log_price_dict).astype(float)

year_price = train.groupby('YearBuilt')['SalePrice'].mean()
year_price = year_price.rolling(window=5, center=True).mean()
year_price = year_price.interpolate(method='linear')
year_price = year_price.fillna(year_price.mean())
data['year_price'] = all_data['YearBuilt'].map(year_price)

data['year_price'].fillna(year_price.mean(), inplace=True)

data = pd.get_dummies(data, drop_first=True, dtype=int)

これらの特徴量をまとめ、処理済みデータとして分割して保存

test = data.iloc[1460:]
train = data.iloc[:1460]

test.to_csv(os.path.join(PROCESSED_PATH, 'test.csv'), index_label='Id')
train.to_csv(os.path.join(PROCESSED_PATH, 'train.csv'), index_label='Id')

次に、提出用の予測データを作成します。

9999_submission.ipynb

import numpy as np
import pandas as pd
import os

from constants import *
test = pd.read_csv(os.path.join(PROCESSED_PATH, 'test.csv'), index_col=0)

モデルをインポート。

import joblib
xgb = joblib.load(os.path.join(MODEL_PATH, '0300_xgb_regr.joblib'))

2.2 ミス発覚！！再学習！！

予測を実行しようとしたところで、問題が発生。XGBoostモデルは訓練データのみで学習を行ったが、ダミー変数化したカテゴリ変数のうち、テストデータにのみ存在するカテゴリのラベルがありました。
テストデータと訓練データを合体したもので再度学習させる必要があります。
先ほどのXGBoostモデルを実装したノートブックに戻り、上で保存した処理済みデータで再度学習させます。

0300_xgboost.ipynb

ref_data = pd.read_csv(os.path.join(PROCESSED_PATH, 'train.csv'), index_col=0)
X = ref_data.drop('log_price', axis=1)
y = ref_data['log_price']

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=22)

best_xgb = XGBRegressor(**best_params, random_state=22, n_jobs=-1)
best_xgb.fit(X_train, y_train)

joblib.dump(best_xgb, os.path.join(MODEL_PATH, '0301_xgb_regr.joblib'))

これで問題ないと思うので、このモデルで予測をしていきます。

2.3 提出データ出力

提出データのカラムはインデックスともとの販売価格とのことだったので、対数で予測した価格を指数変換して、CSVに保存します。

import joblib

# モデルの読み込み
xgb = joblib.load(os.path.join(MODEL_PATH, '0301_xgb_regr.joblib'))

# テストデータを作成
test.drop('log_price', axis=1, inplace=True)

# 予測
preds = xgb.predict(test)

# 指数変換
test['SalePrice'] = np.exp(preds)

# 提出に必要なカラムを取り出し、CSVとして保存
sub = test[['SalePrice']]
sub.to_csv(os.path.join(DATA_PATH, 'submission.csv'), index=True)

3. 結果発表

今回の結果は...

対数化した価格に対してのRMSEが0.14897となりました。
コンペティション的にはそこまで高いスコアではないですが、正則化をうまく利用し、過学習を抑制して学習データでの検証時のRMSEとの差が小さくできたことが、個人的には嬉しいポイントでした。

おわりに：コンペを終えて

今回、何の気になしに始めた「全言語化」だが、このコンペでの試行錯誤や失敗を言語化し続けることで、データ分析からモデル構築までの作業をかなり内面化することができた。
特に、モデルの性能を最大化しつつ過学習を抑止したいという思いで正則化パラメータをいじった結果、しっかりと未知のデータにも対応したモデルができたところが、個人的に嬉しいポイントだった。

ここまでの作業全体を通じて、特徴量選択とエンジニアリング、モデルの構築と解釈といったMLエンジニアの基礎的な知識を再確認することができたのではないかと思う。
また、テストデータと訓練データでダミー変数化により特徴量に差ができてしまうことなどは、実際に手を動かしてみることで初めて経験できるミスだったと感じる。

気が向いたら、他のコンペなどでも「全て言語化する」をしてみてもいいかもしれない。

全4回、ここまでお読みいただきありがとうございました。