水烏 - Kaggleに挑み続けるブログ

とあるインフラ企業のデータ分析担当がKaggleに挑んだその軌跡

不動産価格を機械学習で予測するKaggleに挑戦する [発展編2:上位3%]

移転しました。

前回のは色々複雑すぎる...

不動産価格を機械学習で予測するKaggleのHouse Pricesですが、前回のスコアは0.12を切り、順位も上位20%に入ることができました。
引き換えに色々な試行錯誤、処理、機械学習をテストしたため、chain of thoughts, モデルが複雑になってしまいました。。

今回はできるだけシンプルなデータ処理と予測モデリングでより高い順位を目指したいと思います。
シンプルな予測モデルといえば、やはり線形回帰です。今回はラッソ回帰を中心に考えてみたいと思います。

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

from scipy.stats import skew
from sklearn.linear_model import Lasso
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

import lightgbm as lgb

今回は上記のライブラリ以外使いません。途中から新しくライブラリをインポートすることも行いません。
これで果たして上位に行けるのか、見てみましょう。今回は、できるだけシンプルにやります。

#データの読み込み
df_train = pd.read_csv('train.csv')
df_test = pd.read_csv('test.csv')

#外れ値の処理
df_train = df_train[~((df_train['GrLivArea'] > 4000) & (df_train['SalePrice'] < 300000))]

#Id及びSaleprice以外の全データを結合(左端=Id, 右端はSalePrice)
df_all = pd.concat((df_train.loc[:,'MSSubClass':'SaleCondition'], df_test.loc[:,'MSSubClass':'SaleCondition']))

#Feature Engineering - 前回の中で断トツでFeature Importanceが高かったTotalHousePorchSFを採用
df_all['TotalHousePorchSF'] = df_all['EnclosedPorch']+df_all['OpenPorchSF']+df_all['WoodDeckSF']+df_all['3SsnPorch']+df_all['ScreenPorch']
df_train['TotalHousePorchSF'] = df_train['EnclosedPorch']+df_train['OpenPorchSF']+df_train['WoodDeckSF']+df_train['3SsnPorch']+df_train['ScreenPorch']

外れ値も、前回わかっていた広大なエリアに対して破格の値段がついている(恐らくは何もない土地)の2点だけを除きました。

回帰分析の基礎に立ち返る

線形回帰分析を行う上で、残差の正規性や多重共線性の有無、説明変数ができるだけ正規分布している、線形の関係性があるなど気にすることがあります。
L1正則化項の存在で回帰係数が不安定になりづらいとはいえ、多重共線性の有無、説明変数・目的変数の分布の正規化ぐらいはやっておこうと思います。
尚、不動産価格に対して各種特徴量が概ね線形を保っているのはHouse Pricesについては自明なので、特に触れません。

多重共線性の確認

多重共線性を確認する方法はVIFや相関係数行列を目視するなどいくつかあります。
今回は相関係数行列を適当に眺めて相関が高い奴(0.8+)を省きます。カテゴリカル変数の共線性の評価はかなり複雑になるのでやりません...

#相関係数行列の可視化
fig, ax = plt.subplots(1, 1, figsize=(30, 30))
sns.heatmap(df_train.corr(), vmax=1, vmin=-1, center=0, annot=True, ax=ax)

f:id:prml2synth:20180610104857p:plain

相関係数が絶対値で0.8を超えているものはほぼ確実に多重共線性を引き起こすので、特徴量から削除します。4ペアあります
落とす際は、目的変数であるSalePriceに対しての相関係数が低いものから順に落としていきます。
この場合落とすのは1stFlrSF, TotRmsAbdGrd, GarageArea, GarageYrBltの4つを落とします。

df_all.drop(['1stFlrSF','GarageArea','TotRmsAbvGrd', 'GarageYrBlt'], axis=1, inplace=True)

説明変数・目的変数の対数変換

次に、SalePriceのヒストグラムを見てみます。

df_train["SalePrice"].hist(bins=30)

f:id:prml2synth:20180610104402p:plain

かなり歪度が高いですね。特にこのように高いほうに向かって裾野が広い分布は、対数変換をかけてやるとうまい具合に正規分布してくれます。
ちゃんとQ-Q Plotやシャピロウィルク検定などで正規性の確認をしても良いと思います。

df_train["SalePrice"] = np.log1p(df_train["SalePrice"])

#対数変換後の分布を確認
df_train["SalePrice"].hist(bins=30)

f:id:prml2synth:20180610104415p:plain

それっぽくなりました。尚、対数変換をかける場合は最後の予測の際にEXPを取ることを忘れないようにする必要があります。
同じことをすべてのノンカテゴリカルの特徴量に対して行います。scipyのskewで歪度が計算できるので、
今回は0.6より歪度が大きいものはすべて対数変換をかけることにします。やや雑ですが、シンプルにしておきたいので。

#カテゴリカルでない特徴量
non_categoricals = df_all.dtypes[df_all.dtypes != "object"].index

skewed_feats = df_train[non_categoricals].apply(lambda x: skew(x.dropna()))
skewed_feats = skewed_feats[skewed_feats > 0.6].index

#歪度0.7より大きい特徴量を対数変換
df_all[skewed_feats] = np.log1p(df_all[skewed_feats])

あとはOne Hot Encodingと、欠損値の補完です。最近気づいたんですがHouse Pricesについては欠損値は
全部平均で補完しても精度に影響ないようなので、今回はシンプルにすべて平均値で補完します。

#One Hot Encoding
df_all = pd.get_dummies(df_all)

#欠損値を平均値で補完
df_all = df_all.fillna(df_all.mean())

今回は、バイアス・バリアンスの観点からラッソ回帰と勾配ブースティングの二つのスタッキングを行います。

#学習データ、テストデータに分割
X = df_all[:df_train.shape[0]]
X_for_test = df_all[df_train.shape[0]:]
y = df_train.SalePrice

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1111)
reg = Lasso(alpha=0.0004)
reg.fit(X_train, y_train)
Lasso(alpha=0.0004, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)
y_pred = reg.predict(X_test)
print("ラッソ回帰でのRMSE:",np.sqrt(mean_squared_error(y_pred, y_test)))
ラッソ回帰でのRMSE: 0.11467346011593593
lgb_train = lgb.Dataset(X_train,y_train)
params = {'task': 'train','boosting_type': 'gbdt','objective': 'regression','metric': {'l2'},'num_leaves': 256,
          'learning_rate': 0.01,'num_iterations':2000,'feature_fraction': 0.4,'bagging_fraction': 0.7,'bagging_freq': 5}

gbm = lgb.train(params, lgb_train, num_boost_round=1500)
 
y_test_pred_lgb = gbm.predict(np.array(X_test), num_iteration=gbm.best_iteration)

#RMSE
print("LightGBMでのRMSE:",np.sqrt(mean_squared_error(y_test, y_test_pred_lgb)))
LightGBMでのRMSE: 0.11716578498521343
print("LightGBM+Lassoのスタッキング時のRMSE:",np.sqrt(mean_squared_error(y_test, (y_pred*0.7+y_test_pred_lgb*0.3))))
LightGBM+Lassoのスタッキング時のRMSE: 0.11062571022201569

いくつかテストした結果これがもっとも良さそうなパラメーターでした。
グリッドサーチやCVも熱心には行っていませんが、これでスコアを見てみたいとおもいます。

#全データで学習
reg.fit(X, y)
lgb_train_full = lgb.Dataset(X,y)
gbm = lgb.train(params, lgb_train_full, num_boost_round=1500)
#ラッソ・LightGBMの予測及びスタッキング
pred = np.expm1(reg.predict(X_for_test))
pred2 = np.expm1(gbm.predict(X_for_test))
pred3 = (pred*0.7+pred2*0.3)
solution = pd.DataFrame({"id":df_test.Id, "SalePrice":pred3})
solution.to_csv("test_lasso.csv", index = False)

f:id:prml2synth:20180610104514p:plain

スコアは0.11482、順位は165位/5486位で上位3%にまで食い込みました。
そこまで複雑でない前処理・Feature Engineering、ラッソ回帰と勾配ブースティングのみ使ってシンプルなスタッキングを行うだけでもここまで上がれることを理解頂けたかと思います。
第3回となったHouse Pricesですが、この辺で一度終わりにして次の題材に進みたいと思います。

もし上記を元にスコアを更に上げたい!という場合はやはりIsolationForestなど多変量での外れ値処理、
深いFeature Engineering及びモデルの複雑化(アルゴリズムの追加と複層のスタッキング)などが必須になってくるでしょう。
NetFlixの例を考えても、実務的には上記ぐらいが運用の限界かなと思います。