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

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

不動産価格を機械学習で予測するKaggleに挑戦する [ベンチマーク編]

移転しました。

第一回ということで、Kaggleの不動産価格予測コンペ - House Prices: Advanced Regression Techniquesに挑戦してみたいと思います。
House Prices: Advanced Regression Techniques | Kaggle
タイタニックが分類問題のチュートリアルならば、House Pricesは回帰問題のチュートリアルという立ち位置でしょう。

個人的には分類問題の方が好きですが、早速本題に入ります。

私の場合はデータ分析の流れとして、まず最低限の処理(欠損値及びカテゴリカルデータの処理)だけを行ってベンチマークの精度とし、その後掘り下げるというアプローチをよくとります。
今回はベンチマーク編ということで、最低限の処理のみを行ってどのくらいのスコアになるかを確認したいと思います。

1. データ分析の目的

まず今回のデータ分析の目的を確認します。
家を買うときに何LDKなのか、何坪なのかといったデータはもちろん不動産価格に影響するでしょう。

しかし今回はそれら以外にも多くの要因が不動産価格に影響することを確認します。
具体的には79の変数を元に、不動産価格を予測することになります。データはアメリカ合衆国アイオワ州エイムズの不動産です。

#今回使うライブラリ
import numpy as np
import pandas as pd

from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import Lasso
from sklearn import svm
from sklearn.model_selection import GridSearchCV, train_test_split

2. データセットについて

SalePriceが予測したい不動産価格で、それら以外に1階の広さや地下の天井の高さ、基礎の種類からプールの有無まで多数の特徴量が含まれています。
全部が全部不動産に影響することはないでしょうが、常識的には部屋の広さ、高さ、築年数などは大きく影響するはずです。
特徴量としては面積のように連続値のものもあれば、品質のようなカテゴリカル変数も存在しています。

ベンチマーク編では特徴量同士のデータ可視化や外れ値の処理などは行いません。
79もの特徴量があるため、いくつか重要と考えられる変数だけ抽出したとしてもそこには主観が入ってしまいます。
例えばSalePriceと1回の広さ1stFlrSFをプロットして、外れている点を削除するというのも単変量しか見ておらず、特徴量同士の相互作用を考慮できていません。

データはtrain.csvとtest.csvに分かれており、またKaggleにアップロードする際の形式としてsample_submission.csvが用意されています。
データの説明(Data_Description.txt)も用意されており、必ず読む必要があります。79変数あるのでここでは省きますが、欠損値の補完ですぐに必要になります。

#学習データとテストデータの読み込み
df_train = pd.read_csv("train.csv")
df_test = pd.read_csv("test.csv")

#予測に使わないIdを保持・データフレームから削除
df_train_index = df_train["Id"]
df_test_index = df_test["Id"]
df_train.drop(["Id"], axis=1, inplace=True)
df_test.drop(["Id"], axis=1, inplace=True)
#目的変数であるSalePriceを別に取っておく
y_train = df_train["SalePrice"].values
df_train.drop(["SalePrice"], axis=1, inplace=True)

#学習用データとテストデータを一度統合する
df_all = pd.concat((df_train,df_test)).reset_index(drop=True)
#データ数
print(df_train.shape)
print(df_test.shape)
(1460, 79)
(1459, 79)

データの数は学習用データ・テストデータともにほぼ同数。学習用データが1500レコード弱で比較的扱いやすいですね。

3. 欠損値の補完

データセットの中身を見てみると、結構データが存在しない(NaN)レコードが多いことに気が付きます。これは本当に欠損しているのか、あるいは単純に"該当なし"という意味なのか紐解きながら補完する必要があります。
いずれかのデータが欠損しているレコードを削除してしまうというのは一つの手ですが、可能な限り手元にあるデータを生かすためまずは補完します。

#欠損値の個数確認 / 2919データ
df_all.isnull().sum()[df_all.isnull().sum() != 0].sort_values(ascending=False)
PoolQC          2909
MiscFeature     2814
Alley           2721
Fence           2348
FireplaceQu     1420
LotFrontage      486
GarageFinish     159
GarageYrBlt      159
GarageQual       159
GarageCond       159
GarageType       157
BsmtExposure      82
BsmtCond          82
BsmtQual          81
BsmtFinType2      80
BsmtFinType1      79
MasVnrType        24
MasVnrArea        23
MSZoning           4
BsmtFullBath       2
BsmtHalfBath       2
Utilities          2
Functional         2
Exterior2nd        1
Exterior1st        1
SaleType           1
BsmtFinSF1         1
BsmtFinSF2         1
BsmtUnfSF          1
Electrical         1
KitchenQual        1
GarageCars         1
GarageArea         1
TotalBsmtSF        1
dtype: int64

こうしてみると結構な特徴量で欠損値があるようですが、例えばPoolQCはData Descriptionを読むとPoolQCがNAならプールなしという意味のようです。
ほかにもガレージ系でそれぞれ159個ずつ仲良く欠損していたりするので、値がないというより、"NA"のつもりが欠損値になっているのでしょう。
それでは早速上記リストとData Descriptionを見ながら欠損値を補完していきましょう。

#欠損値の補完

#以下はNaN = NAかNoneの特徴量リスト。よって欠損値をそれぞれNAとNoneで補完する。
df_all["PoolQC"].fillna('NA', inplace=True)
df_all["MiscFeature"].fillna('None', inplace=True)
df_all["Alley"].fillna('NA', inplace=True)
df_all["Fence"].fillna('NA', inplace=True)
df_all["FireplaceQu"].fillna('NA', inplace=True)
df_all["GarageQual"].fillna('NA', inplace=True)
df_all["GarageFinish"].fillna('NA', inplace=True)
df_all["GarageCond"].fillna('NA', inplace=True)
df_all["GarageType"].fillna('NA', inplace=True)
df_all["BsmtCond"].fillna('NA', inplace=True)
df_all["BsmtExposure"].fillna('NA', inplace=True)
df_all["BsmtQual"].fillna('NA', inplace=True)
df_all["BsmtFinType2"].fillna('NA', inplace=True)
df_all["BsmtFinType1"].fillna('NA', inplace=True)
df_all["MasVnrType"].fillna('None', inplace=True)

#以下はNaN = 0の特徴量リスト。例えば地下なら、地下がないんだから0。みたいな。
df_all["GarageYrBlt"].fillna(0, inplace=True) # ガレージ築年数を0にするのも不思議な気はしますが、そもそもガレージがないので他に妥当な数字が思いつかず。
df_all["MasVnrArea"].fillna(0, inplace=True)
df_all["BsmtHalfBath"].fillna(0, inplace=True)
df_all["BsmtFullBath"].fillna(0, inplace=True)
df_all["TotalBsmtSF"].fillna(0, inplace=True)
df_all["BsmtUnfSF"].fillna(0, inplace=True)
df_all["BsmtFinSF2"].fillna(0, inplace=True)
df_all["BsmtFinSF1"].fillna(0, inplace=True)
df_all["GarageArea"].fillna(0, inplace=True)
df_all["GarageCars"].fillna(0, inplace=True)

#欠損レコード数が少なく、大半が一つの値をとっているためあまりに予測の役に立たなさそうな特徴量は単純に最頻値を代入
df_all["MSZoning"].fillna('RL', inplace=True)
df_all["Functional"].fillna('Typ', inplace=True)
df_all["Utilities"].fillna("AllPub", inplace=True)
df_all['SaleType'] = df_all['SaleType'].fillna(df_all['SaleType'].mode()[0])
df_all['Exterior2nd'] = df_all['Exterior2nd'].fillna(df_all['Exterior2nd'].mode()[0])
df_all['Exterior1st'] = df_all['Exterior1st'].fillna(df_all['Exterior1st'].mode()[0])
df_all['KitchenQual'] = df_all['KitchenQual'].fillna(df_all['KitchenQual'].mode()[0])
df_all['Electrical'] = df_all['Electrical'].fillna(df_all['Electrical'].mode()[0])

#LotFrontage - Linear feet of street connected to property
#これは補完方法が明らかかつ簡単で、近くのStreet名=Neighborhoodでグループし平均を取れば良い精度で補完できそう。
f = lambda x: x.fillna(x.mean())
df_all["LotFrontage"] = df_all.groupby("Neighborhood")["LotFrontage"].transform(f)
#欠損値がすべて補完されているか確認
df_all.isnull().sum()[df_all.isnull().sum() != 0].sort_values(ascending=False)
Series([], dtype: int64)

欠損レコードが消えていますね。

かなり地道な作業になりますが、特に欠損レコード数の多いものに対してどう補完するかで他のデータサイエンティストと差別化できます。
平均値が良いのか最頻値が良いのかなどはトライ&エラーしてみるのが一番簡単です。
もちろん離散値の場合に平均値を使うのは感覚的に変だとか、そういった判断をもとに決めるのも良いと思います

4.前処理

すぐに機械学習したいところですが、カテゴリー変数の処理を行います。

カテゴリー変数の扱い方はLabel EncodingやOne Hot Encodingなど様々な方法が存在します。
今回は大したレコード数ではないことと、特徴量の数も1000すらいかないと思われるのでOne Hot Encodingでダミー変数化します。

#One Hot Encoding
df_all = pd.get_dummies(df_all)
#特徴量数の確認
df_all.shape
(2919, 302)

学習データ・テストデータ合わせて3000レコード、300特徴量です。ややスパースなデータになっていそうですが、ベンチマーク編ではこれ以上の処理は行いません。

5. 機械学習で不動産価格を予測する

楽しいパートがやってきました。みんな大好き機械学習です。
ベンチマーク編なので、とりあえずラッソ回帰、ランダムフォレスト、SVRでやろうと思います。

チューニングはとりあえずざっくり辺りをつけてグリッドサーチします。

#学習用データとテストデータを切り分ける
ntrain = df_train.shape[0]
train = df_all[:ntrain]
test = df_all[ntrain:]

y = y_train
X = train.loc[:, train.columns != 'SalePrice']

#30%でチューニングを行う
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1234)
#モデルの呼び出し
lasso = Lasso()
rf = RandomForestRegressor()
svr = svm.SVR()
#グリッドサーチ用パラメータの設定
lasso_parameters = {'alpha':[0.1, 0.5, 1]}
rf_parameters= {'n_estimators':[100, 500, 2000], 'max_depth':[3, 5, 10]}
svr_parameters = {'C':[1e-1, 1e+1, 1e+3], 'epsilon':[0.05, 0.1, 0.3]} 

#グリッドサーチ
lasso_gs = GridSearchCV(lasso, lasso_parameters)
lasso_gs.fit(X_train,y_train)

rf_gs = GridSearchCV(rf, rf_parameters)
rf_gs.fit(X_train,y_train)

svr_gs = GridSearchCV(svr, svr_parameters)
svr_gs.fit(X_train,y_train)
GridSearchCV(cv=None, error_score='raise',
       estimator=SVR(C=1.0, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma='auto',
  kernel='rbf', max_iter=-1, shrinking=True, tol=0.001, verbose=False),
       fit_params=None, iid=True, n_jobs=1,
       param_grid={'C': [0.1, 10.0, 1000.0], 'epsilon': [0.05, 0.1, 0.3]},
       pre_dispatch='2*n_jobs', refit=True, return_train_score='warn',
       scoring=None, verbose=0)

それでは、予測してみたいと思います。

#ラッソ回帰
y_pred = lasso_gs.predict(X_test)
print("ラッソ回帰でのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred)))

#ランダムフォレスト
y_pred2 = rf_gs.predict(X_test)
print("ランダムフォレストでのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred2)))

#SVR
y_pred3 = svr_gs.predict(X_test)
print("SVRでのRMSE:",np.sqrt(mean_squared_error(y_test, y_pred3)))
ラッソ回帰でのRMSE: 38914.187417202746
ランダムフォレストでのRMSE: 24883.658315015313
SVRでのRMSE: 72294.23251731148

安定のランダムフォレストですね。今回のデータセットではやはり正規化などの処理を行わないとSVRは辛そうです。
では、実際にどのような値になっているか確認します。

df_visualize = pd.concat((pd.DataFrame(y_test), np.round(pd.DataFrame(y_pred))), axis=1)
df_visualize.columns = ['実際の値','予測値']
df_visualize.head(10)
実際の値 予測値
0 205000 203899.0
1 345000 355549.0
2 173900 183026.0
3 93500 45687.0
4 265900 223558.0
5 212000 228800.0
6 221000 257087.0
7 102000 108516.0
8 290000 328641.0
9 140000 156622.0

早速4つ目のデータは大外ししていますね。ほぼ半分。
これではスコアもあまり期待できなさそうですが、まずは確認したいと思います。
切り分けておいたテストデータで提出用のSalePrice表を作ります。

y_pred_final = rf_gs.predict(test)

次に提出用データのフォーマットを確認し、合わせます。

pd.read_csv('sample_submission.csv').head()
Id SalePrice
0 1461 169277.052498
1 1462 187758.393989
2 1463 183583.683570
3 1464 179317.477511
4 1465 150730.079977

予測値と、最初の方で分けておいたIdを合体させる。

submission = pd.concat((df_test_index, pd.DataFrame(y_pred_final)), axis=1)
submission.columns = ['Id', 'SalePrice']
#確認
submission.head()
Id SalePrice
0 1461 129225.399612
1 1462 159871.544036
2 1463 173464.563263
3 1464 183416.153801
4 1465 208079.458478

それでは、csvに吐き出してKaggleにアップロードします。
本来もっともワクワクする瞬間ですが、今回ベンチマーク編ではあまり期待できないのがわかっている為流れ作業的に...

submission.to_csv("submission.csv",sep=',',index=False)

では、一体スコアは何点だったでしょうか?
f:id:prml2synth:20180601225945p:plain


結果は2690位で、2018/6/1時点で参加者は5458チームのため、ちょうど真ん中あたりですね。

ベンチマーク編では予想通りスコアは残念な結果となってしまいましたが、次の発展編でデータ探査、精度評価及びチューニングを追求し、
どこまで順位を伸ばせるのか試してみたいと思います。

発展編1(上位20%)書きました。
www.mizugarasu.jp
発展編2(上位3%)やってみました。
www.mizugarasu.jp