Kaggle - House Prices挑戦記②：ベースラインモデル作成

0. はじめに

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

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

前回は、住宅価格に影響を与えそうなさまざまに特徴量を作っては可視化しました。
今回は、前回作った特徴量からベースラインモデルを作り、特徴量の貢献度などを測っていきたいと思います。

1. ベースラインモデル選定

今回は、Ridge回帰をベースラインモデルとしました。その理由は、

1. 多重共線性への対処：面積関連変数/フラグや古さ指標など、相関の高そうな変数が多いため、正則化モデルが適切
2. 解釈のしやすさ：線形モデルはモデルの解釈がしやすい。どの特徴量が重要かが係数から分かる
3. 計算効率：非線形モデルに比べて学習が速いから、素早く結果を確認できる
4. 過学習の防止：正則化パラメータを調整することで、過学習を防げる

それでは、早速ベースラインモデルにデータを読み込ませるため、データを加工していきます。

2. データ加工

2.1 データの取得〜分割

まずは、前回作成したデータを取得します。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

all_data = pd.read_csv('./data/raw/0000_train.csv', index_col=0)
all_data.head() # 110 columns

これら全てのデータから、前回作成した価格決定に寄与しそうな変数のみを抽出します。

features = [
    'log_total_sf', 'log_lot_area', 'house_age', 'remod_age',
    'year_price', 'GarageCars', 'MSSubClass_weight', 'MSZoning_weight', 
    'neighbor_weight', 'has_2nd', 'has_bsmt', 'has_garage', 'has_pool', 'has_central_ac',
    'OverallQual_weight', 'OverallCond_weight', 'ExterQual_weight', 'ExterCond_weight',
    'BsmtQual_weight', 'BsmtCond_weight',  'KitchenQual_weight', 'Utilities_weight',
    'GarageQual_weight', 'GarageCond_weight', 'PoolQC_weight', 'Functional_weight'
    ]

# 説明変数
data = all_data[features]

# 目的変数
target = all_data['log_price']

次に、訓練データと検証データに分割します。

X = data
y = target

from sklearn.model_selection import train_test_split

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

続いて、データを標準化していきます。係数化しているものはすでに標準化されているため、標準化しません。また、フラグ変数についても標準化をしません。解釈が難しくなりますからね。

# 標準化する変数
cols_to_scale = [
    'log_total_sf', 'log_lot_area', 'GarageCars'
    'house_age', 'remod_age', 'year_price'
    ] 

# 標準化しない係数
cols_to_scale = [
    'log_total_sf', 'log_lot_area', 'GarageCars',
    'house_age', 'remod_age', 'year_price'
    ] 

cols_not_scale = [
    'MSSubClass_weight', 'MSZoning_weight', 'neighbor_weight' ,
    'has_2nd', 'has_bsmt', 'has_garage', 'has_pool', 'has_central_ac'
    'OverallQual_weight', 'OverallCond_weight', 'ExterQual_weight', 'ExterCond_weight',
    'BsmtQual_weight', 'BsmtCond_weight',  'KitchenQual_weight', 'Utilities_weight',
    'GarageQual_weight', 'GarageCond_weight', 'PoolQC_weight', 'Functional_weight'
    ]

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()

X_train[cols_to_scale] = scaler.fit_transform(X_train[cols_to_scale])
X_val[cols_to_scale] = scaler.transform(X_val[cols_to_scale])

2.2 テストモデルの学習と評価

さて、データの準備ができたので、まずはこれらを全て読ませて、正則化パラメータ$\alpha=0$（すなわち単純な線形回帰）でテストモデルを作る。
なお、このコンペではスコアは対数変換した価格についてのRMSE（二乗平均平方根誤差）によって測定するため、log_priceをそのまま目的変数とする。

from sklearn.linear_model import Ridge

alpha = 0 # 正則化パラメータはいったん0

# インスタンスを生成
test_model = Ridge(alpha=alpha)

# 学習
test_model.fit(X_train, y_train)

# 予測を実行
y_pred = test_model.predict(X_val)

# 評価
from sklearn.metrics import root_mean_squared_error, r2_score

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

print(f"Ridge Regr. (alpha={alpha})")
print(f'  RMSE:\t{rmse:.4f}')
print(f'  R2:\t{r2:.4f}')

# Ridge Regr. (alpha=0)
#  RMSE:  0.1441
#  R2:	  0.8824

すでにRMSEが0.14、決定係数が0.88とかなり良い成績な気もしますが、今回の目標は「高いスコアを出すこと」ではなく「思考すること」「思考を言語化すること」なので、続行します。

今回のあてはまり度合いを可視化していきます。
seabornのregplotを使用し、予測と実測の散布図と回帰直線（95%信頼区間）を描画。

sns.regplot(x=y_pred, y=y_val,
            scatter_kws={'s': 20, 'alpha': 0.3},
            line_kws={'color': 'green', 'lw': 1},
            ci=95)
            
plt.title(f"Test Model: Ridge Regr. (alpha={alpha})")
plt.show()

3. 特徴量評価

続いて、このテストモデルにおいて各変数がどれくらい寄与しているかを視覚化します。

plt.figure(figsize=(12, 8))

# 回帰係数を取得
coef = test_model.coef_

# 重要度データフレームを作成
importance = pd.DataFrame({
    'feature': features,
    'importance': np.abs(coef),
    'coefficient': coef
    }).sortv_alues('importance', ascending=False)

# 係数の符号により色分け
colors = ['red' if c<0 else 'blue' for c in importance['coefficient']]

# 凡例を追加
from matplotlib.patches import Patch
legends = [
    Patch(facecolor='blue', label='positive influence'),
    Patch(facecolor='red', label='negative influence')
]

# 表示
sns.barplot(x='importance', y='feature', data=importance, palette=colors)
plt.legend(handles=legends, loc='lower right')
plt.xlabel('Importance')
plt.ylabel('Features')
plt.tight_layout()
plt.show()

数値も確認してみましょう。

print(importance_df)

3.1 貢献度からわかること

• 品質、床面積、全館空調あたりが特に重要
• 負の貢献をしている意外な変数がある（プール、地下室、ガレージの有無、年代別ロール平均など）

今後の対応

1. 前回の特徴量作成時に、BsmtQual_WeightやPoolQC_weightなどの変数の欠損値に対し、それらは設備がないことを示していることがわかったため'NONE'で補填した。これらの情報がhas_bsmtやhas_poolなどの持つ情報と被るため、フラグを削除していく。
2. 似通った変数がある（〇〇Condと〇〇Qual系）ので、それらの相関を調べ、必要に応じてまとめる。

まずは、不要なカラムを削除。

# 1.不要なカラムを削除
X_train_dropped = X_train.drop(columns=['has_pool', 'has_bsmt', 'has_garage'])
X_val_dropped = X_val.drop(columns=['has_pool', 'has_bsmt', 'has_garage'])

次に、相関を似通った変数同士（品質関連）の連関を調べる。

# 順序尺度同士の独立性をカイ二乗検定（有意水準0.01）

sim_cols_headers = ['Overall', 'Exter', 'Garage', 'Bsmt']
sim_cols_footers = ['Cond_weight', 'Qual_weight']

from scipy.stats import chi2_contingency

for i, h in enumerate(sim_cols_headers):
    
    x = h + sim_cols_footers[0]
    y = h + sim_cols_footers[1]
    
    cross_tab = pd.crosstab(X_train_dropped[x], X_train_dropped[y])
    
    chi2, p, _, _ =chi2_contingency(cross_tab)
    
    if p < 0.01:
        print(f'"{h}" correlated / p = {round(p, 4)} < 0.01')
    else:
        print(f'"{h}" not correlated / p = {round(p, 4)} > 0.01')

# 出力
# "Overall" correlated / p = 0.0 < 0.01
# "Exter" correlated / p = 0.0 < 0.01
# "Garage" correlated / p = 0.0 < 0.01
# "Bsmt" correlated / p = 0.0 < 0.01

だいぶ強い連関が見られる。いったんフラグ変数を削除したデータで、再度学習をしてから重要度を再度調べる。

test_model2 = Ridge(alpha=alpha, random_state=1)
test_model2.fit(X_train_dropped, y_train)

y_pred2 = test_model2.predict(X_val_dropped)

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

print(f"Ridge Regr. (alpha={alpha})")
print(f"RMSE: {rmse:.3f}, R2: {r2:.3f}")

# Ridge Regr. (alpha=0)
# RMSE: 0.145, R2: 0.881

plt.figure(figsize=(12, 8))

# 回帰係数を取得
coef = test_model2.coef_

# 重要度データフレームを作成
importance = pd.DataFrame({
    'feature': features,
    'importance': np.abs(coef),
    'coefficient': coef
    }).sortv_alues('importance', ascending=False)

# 係数の符号により色分け
colors = ['red' if c<0 else 'blue' for c in importance['coefficient']]

# 凡例を追加
from matplotlib.patches import Patch
legends = [
    Patch(facecolor='blue', label='positive influence'),
    Patch(facecolor='red', label='negative influence')
]

# 表示
sns.barplot(x='importance', y='feature', data=importance, palette=colors)
plt.legend(handles=legends, loc='lower right')
plt.xlabel('Importance')
plt.ylabel('Features')
plt.tight_layout()
plt.show()

ここで、重要度の高い変数順に、相関行列と偏相関行列を作成して変数同士の関係を再度見直す。

important_cols = importance_df['feature'].values

corr_mat = X_train_dropped[important_cols].corr().round(3)

plt.figure(figsize=(18, 16))
sns.heatmap(data=corr_mat, annot=True, cmap='coolwarm')
plt.title('Correlation Matrix: Ridge Regr.', fontsize=20, loc='left')
plt.savefig(os.path.join(IMAGE_PATH, '0101_corr_matrix.png'))
plt.show()

3.2 相関図・編相関図から考えられること

• 相関も偏相関も高い：それらについてはまとめることが可能か
• 相関も偏相関も低い：独立と仮定

より具体的に

• 「全体の品質」と「広さ」には本質的な関係がありそう（相関・偏相関ともにある程度ある）
• ガレージのConditionとQualityについては非常に強い相関があるため、加算してまとめる
• house_ageとyear_priceについては、より強い重要度をもつyear_priceのみで進めてみる（house_ageは他の変数とも比較的相関が高い）。

次のステップ：

• 広さと品質は立地（高級住宅街など）が説明しているのでは？ -> 独立性の仮説検定

quality = X_train_dropped['OverallQual_weight']
size = X_train_dropped['log_total_sf']
location = X_train_dropped['neighbor_weight'] 

# 通常の相関（立地の効果を考慮しない）
corr_qs = np.corrcoef(quality, size)[0, 1]
print(f"品質と広さの相関: {corr_qs:.4f}")

# 部分相関（立地の効果を制御）
pcorr_qs_given_l = pg.partial_corr(data=X_train_dropped, x='OverallQual_weight', y='log_total_sf', covar='neighbor_weight')
hdi = pcorr_qs_given_l['CI95%'].iloc[0]
print(f"立地を条件とした品質と広さの部分相関（95%信頼区間）: {hdi[0]}~{hdi[1]}")

print(f'相関減少率: {((corr_qs- hdi[1])/corr_qs*100).round(2)}~{((corr_qs - hdi[0])/corr_qs*100).round(2)}%')

• 出力
  品質と広さの相関: 0.6693
  立地を条件とした品質と広さの部分相関（95%信頼区間）: 0.45~0.53
  相関減少率: 20.82~32.77%

• 解釈：品質と広さの関係の約20.8~32.8%が立地によって説明されるが、残りの7〜8割程度は立地以外の要因による本質的な関連と考えられる。

とりあえずは、第三の変数によって完全に媒介されていなさそう。

4. ベースラインモデル作成

重複する指標と重要でない指標を削除したデータセットで、新たに学習する

X_train_dropped['Garage_weight'] = X_train_dropped['GarageCond_weight'] + X_train_dropped['GarageQual_weight']
X_val_dropped['Garage_weight'] = X_val_dropped['GarageCond_weight'] + X_val_dropped['GarageQual_weight']

X_train_dropped.drop(['GarageCond_weight', 'GarageQual_weight', 'house_age'], axis=1, inplace=True)
X_val_dropped.drop(['GarageCond_weight', 'GarageQual_weight', 'house_age'], axis=1, inplace=True)

ここで、ベストな正則化パラメータを設定するために、グリッドサーチを使用。

alphas = np.logspace(-3, 5, 81)

from sklearn.model_selection import GridSearchCV

gs = GridSearchCV(
    estimator=Ridge(),
    param_grid={'alpha': alphas},
    cv=5,
    scoring='neg_root_mean_squared_error'
)

gs.fit(X_train_dropped, y_train)

best_alpha = gs.best_params_['alpha']
best_ridge = gs.best_estimator_

念の為、アルファの値によって変数の計数がどのように変化しているかを可視化。

# 各アルファに対する係数を保存するリスト
coefs = []

# 各アルファで学習して係数を保存
for alpha in alphas:
    ridge = Ridge(alpha=alpha)
    ridge.fit(X_train_dropped, y_train)
    coefs.append(ridge.coef_)

# 結果を配列に変換
coefs = np.array(coefs)

# 係数の変化をプロット
plt.figure(figsize=(10, 6))
for i in range(X_train_dropped.shape[1]):
    plt.plot(alphas, coefs[:, i], label=f"{X_train_dropped.columns[i]}")

plt.xscale('log')
plt.xlabel('alpha')
plt.ylabel('Coef')
plt.title('Ridge regr.')
plt.legend(loc='upper left')
plt.grid(True)
plt.show()

モデルの評価を、CSVとして保存。

y_pred = best_ridge.predict(X_val_dropped)

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

metrics = {
    'model': 'Ridge (base)',
    'params': f'alpha={best_alpha}',
    'RMSE': round(rmse, 4),
    'R2': round(r2, 4)
}

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

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

print(f"RMSE: {rmse:.3f}, R2: {r2:.3f}")

# RMSE: 0.144, R2: 0.882

最後に、特徴量重要度を確認する。

plt.figure(figsize=(12, 8))

coef = best_ridge.coef_

importance = pd.DataFrame({
    'feature': features,
    'importance': np.abs(coef),
    'coefficient': coef
    }).sortv_alues('importance', ascending=False)

# 係数の符号により色分け
colors = ['red' if c<0 else 'blue' for c in importance['coefficient']]

# 凡例を追加
from matplotlib.patches import Patch
legends = [
    Patch(facecolor='blue', label='positive influence'),
    Patch(facecolor='red', label='negative influence')
]

# 表示
sns.barplot(x='importance', y='feature', data=importance, palette=colors)
plt.legend(handles=legends, loc='lower right')
plt.xlabel('Importance')
plt.ylabel('Features')
plt.tight_layout()
plt.show()

学習曲線を確認。

def plot_learning_curve(estimator, X, y, title="Learning Curve"):
    from sklearn.model_selection import learning_curve
    import matplotlib.pyplot as plt
    
    # 学習曲線の計算
    train_sizes, train_scores, val_scores = learning_curve(
        estimator, X, y,
        train_sizes=np.linspace(0.1, 1.0, 10),
        cv=5,
        scoring='neg_mean_squared_error',
        n_jobs=-1
    )
    
    train_mean = train_scores.mean(axis=1)
    train_std = train_scores.std(axis=1)
    val_mean = val_scores.mean(axis=1)
    val_std = val_scores.std(axis=1)
    
    
    # プロット作成
    plt.plot(train_sizes, train_mean, 
                marker='o', label='Training')
    plt.plot(train_sizes, val_mean, ls='--',
                marker='o', label='Validation')
    
    # 信頼区間の追加
    plt.fill_between(train_sizes, 
                     train_mean - train_std,
                     train_mean + train_std, 
                     alpha=0.1)
    plt.fill_between(train_sizes,
                     val_mean - val_std,
                     val_mean + val_std,
                     alpha=0.1)
    
    # グラフの装飾
    plt.title(title)
    plt.xlabel('Training Examples')
    plt.ylabel('Neg Mean Squared Error')
    plt.legend()
    plt.tight_layout()

plot_learning_curve(Ridge(best_alpha), X_train, y_train,'Ridge Regr.')
plt.show()

訓練初期はオーバーフィッティング気味だが、正則化により徐々にMSEが上がり、最終的には落ち着いている。
訓練と検証でのMSEの差が小さく、概ね過学習を避けたモデルができた。

ベースラインモデルは、以下になった。
alpha=10
特徴量：
['log_total_sf', 'log_lot_area', 'remod_age', 'year_price', 'GarageCars',
'MSSubClass_weight', 'MSZoning_weight', 'neighbor_weight', 'has_2nd',
'has_central_ac', 'OverallQual_weight', 'OverallCond_weight',
'ExterQual_weight', 'ExterCond_weight', 'BsmtQual_weight',
'BsmtCond_weight', 'KitchenQual_weight', 'Utilities_weight',
'PoolQC_weight', 'Functional_weight', 'Garage_weight']
スコア：RMSE - 0.145, R2 - 0.881
次回は、このベースラインモデルを基準に、ランダムフォレスト回帰を実装する。