Kaggleの練習問題（Regression）を解いてKagglerになる

　Kaggleの練習問題の1つである、House Pricesに取り組んでみます。Regressionの練習問題はこれ１つですので、がっつり取り組んで他の（お金の絡む）コンペのための準備をしたいですね笑　

　使用言語はPythonです。基本的に、自分のKernelと内容は同じです。Qiitaはいいのでこっちにupvoteしてくれると嬉しいです。

EDA (Exploratory Data Analysis)

　まずは、データがどういう構造なのか探っていきます。Dataの項目を見ますと、ご丁寧に各列の変数の説明が載っています。例えば、今回の予測すべきターゲットは、

SalePrice - the property's sale price in dollars. This is the target variable that you're trying to predict.

ですね。他、家に関する様々な項目に関するデータから、その家の価格を推定するモデルを作るというのが、今回の課題です。いきなり70近い項目（列）が与えられるので、初めてだと面食らってしまいますが、順番に片付けていきましょう。

1) データタイプを調べる

　まずは、与えられたデータがざっくり数値なのか、それとも文章なのか、ざっくりデータタイプを把握しましょう。pandasでcsvデータを読み込み、df.dtypesで各列のデータタイプを取得します。

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)

# for visualization
import matplotlib.pyplot as plt
import seaborn as sns

# import data
train = pd.read_csv('../input/train.csv')
test = pd.read_csv('../input/test.csv')

print(train.dtypes)

　int64、float64はざっくり言うと数値なので、そのままモデルに投げることができます。ただ、このデータはobjectがとても多いですね。objectは基本的には文字情報で、このままだとモデルに与えることはできません。そこで登場するのが、Label Encoder です。

2) Label Encoder を使う

# string label to categorical values
from sklearn.preprocessing import LabelEncoder

for i in range(train.shape[1]):
    if train.iloc[:,i].dtypes == object:
        lbl = LabelEncoder()
        lbl.fit(list(train.iloc[:,i].values) + list(test.iloc[:,i].values))
        train.iloc[:,i] = lbl.transform(list(train.iloc[:,i].values))
        test.iloc[:,i] = lbl.transform(list(test.iloc[:,i].values))

　sklearnに入っているラベルエンコーダーを使うことにより、オブジェクトをカテゴリーの数値に変換することができます。例えば、

　は、ラベルエンコーダーにより、

　と変換されます。これで扱いやすくなりましたね。

3) NaNに対処する

　データを集めるのは大変です。多くの場合、完全にデータを取ることはできずに欠損値NaNが生じてしまいます。これにどう対処するかですが、まずはどれだけデータにNaNがあるのか見てみましょう。

# search for missing data
import missingno as msno
msno.matrix(df=train, figsize=(20,14), color=(0.5,0,0))

　missingnoという素晴らしいライブラリの力により、マトリックスのどこにNaNがあるのか一瞬で視覚化できます。白い場所が、データが欠損している場所です。これらのNaNをどうするかですが、今回はNaNがtraining, testデータともに多いものはマトリックスから落とし、それ以外はその列の中央値 (median)で埋めることにしました。

# keep ID for submission
train_ID = train['Id']
test_ID = test['Id']

# split data for training
y_train = train['SalePrice']
X_train = train.drop(['Id','SalePrice'], axis=1)
X_test = test.drop('Id', axis=1)

# dealing with missing data
Xmat = pd.concat([X_train, X_test])
Xmat = Xmat.drop(['LotFrontage','MasVnrArea','GarageYrBlt'], axis=1)
Xmat = Xmat.fillna(Xmat.median())

4) 新しい列 (feature) を作る

　これでマトリックスはきれいになりましたが、これまでは与えられたデータのクリーニングをしたにすぎません。良い予測をするためには、データの性質を表すFeatureが不可欠です。既存のFeature（列）を組み合わせて、新しいFeatureを作ります。

　今回の場合、'...SF'という名のFeatureが複数登場するので、それらをまとめたFeatureを作りましょう。

5) ターゲットが正規分布に従っているかチェックする

　Regressionの場合、ターゲットの値がどちらかに偏らず正規分布に従っていることが重要です。今回のターゲット、SalePriceは正規分布に従っているでしょうか？

　……従って……ない……笑

　このままだとモデルのパフォーマンスに影響するので、logを使って疑似的に正規分布に持っていきます。

　はい、これでターゲットが（だいたい）正規分布になりました。

6) どのFeatureが大切なのかチェックする

　現在、70をこえるFeatureがあります。何も考えずに全てをモデルに投げるのも手ですが、逆に、モデルにどのFeatureが予測に大切か教えてもらいましょう。今回は、ランダムフォレスト (Regressor)を使います。

# feature importance using random forest
from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor(n_estimators=80, max_features='auto')
rf.fit(X_train, y_train)
print('Training done using Random Forest')

ranking = np.argsort(-rf.feature_importances_)
f, ax = plt.subplots(figsize=(11, 9))
sns.barplot(x=rf.feature_importances_[ranking], y=X_train.columns.values[ranking], orient='h')
ax.set_xlabel("feature importance")
plt.tight_layout()
plt.show()

　sklearnのランダムフォレストのfeature_importance_を見ることにより、どのFeatureがどれくらい大切か、概観を得ることができます。

　……って、上2つ ('OverallQual' & 'TotalSF')がぶっちぎりじゃなイカ！ｗｗ

　つまり、70以上あるFeatureの中で、SalePriceの予測に大切なのはほんの数個しかない、ということがわかります。

　まぁでも、とりあえず上30個のFeatureを取ってきて、それをモデルに与え、訓練することにします。そして、上2つのFeatureを掛け合わせたものを、'Interaction'として新しいFeatureを作ります。

# use the top 30 features only
X_train = X_train.iloc[:,ranking[:30]]
X_test = X_test.iloc[:,ranking[:30]]

# interaction between the top 2
X_train["Interaction"] = X_train["TotalSF"]*X_train["OverallQual"]
X_test["Interaction"] = X_test["TotalSF"]*X_test["OverallQual"]

7) 各Featureとターゲットの関係を調べる

　現在、Interactionを入れてもせいぜい31個しかFeatureがないので、vs SalePriceをすべてプロットしてみましょう。

# relation to the target
fig = plt.figure(figsize=(12,7))
for i in np.arange(30):
    ax = fig.add_subplot(5,6,i+1)
    sns.regplot(x=X_train.iloc[:,i], y=y_train)

plt.tight_layout()
plt.show()

とりあえずInteractionは抜いて、各FeatureをX軸、SalePriceをY軸に散布図を描きました。だいたい線形の関係性が読み取れますが、いくつか外れ値 (outlier)があるようです。それらを除き、よりデータをきれいにしましょう。

8) 外れ値を除く

　マニュアルで、'TotalSF'と'GrLivArea'に見られた外れ値をマトリックスから除きます。

# outlier deletion
Xmat = X_train
Xmat['SalePrice'] = y_train
Xmat = Xmat.drop(Xmat[(Xmat['TotalSF']>5) & (Xmat['SalePrice']<12.5)].index)
Xmat = Xmat.drop(Xmat[(Xmat['GrLivArea']>5) & (Xmat['SalePrice']<13)].index)

# recover
y_train = Xmat['SalePrice']
X_train = Xmat.drop(['SalePrice'], axis=1)

　以上で、EDAパートは終了です。とうとう、モデリングの段階へと移行します。

Stacking-Emsemble modelを作る

　例えば、Deep Learningだけを使って予測値を求めてもいいのですが、多くの場合、アンサンブル (Emsemble)といって、複数のモデルの予測結果を平均した方がより良い予測に繋がる場合が多いです。今回は、機械学習モデルの最終兵器と名高いXGBoost、最近はみんな知っているNeural Network、昔も今もずっと活躍Support Vector Regressorの3つのモデルを別々に訓練し、結果を平均して提出したいと思います。

　ただ、もしかしたらXGBoostが他のモデルより圧倒的に良かったりするかもしれません。その場合、XGBoostの結果をより重んじたいですよね。つまり、

y_{prediction} = w_1 * XGB + w_2 * NN + w_3 * SVR

　のときの$w_1$ を他よりも大きくしたい。

「あれっ……これってもう一つのRegressionなんじゃないの？」

　と最初に考えた人はすごいですよね。これがスタッキング (Stacking)と呼ばれる手法で、複数のモデルの加重平均 (weighted average)によって全体の予測値を決定します。つまり、各モデルの予測値を新しい、第２のマトリックスと捉え、再びRegressionの手法を用いて最適な$w$ の値を決めていくわけです。

　ざっくり、各モデルによるトレーニングは以下のようになります。いずれのモデルも複数のハイパーパラメーター (hyperparameter) が存在する複雑なモデルですので、グリッドサーチ (Grid Search)を用いて、しらみつぶしに最適なモデルパラメーターを決定しようと試みています。

XGBoost

# XGBoost
import xgboost as xgb
from sklearn.model_selection import GridSearchCV

print("Parameter optimization")
xgb_model = xgb.XGBRegressor()
reg_xgb = GridSearchCV(xgb_model,
                   {'max_depth': [2,4,6],
                    'n_estimators': [50,100,200]}, verbose=1)
reg_xgb.fit(X_train, y_train)

Neural Network

from keras.models import Sequential
from keras.layers import Dense
from keras.wrappers.scikit_learn import KerasRegressor

def create_model(optimizer='adam'):
    model = Sequential()
    model.add(Dense(X_train.shape[1], input_dim=X_train.shape[1], kernel_initializer='normal', activation='relu'))
    model.add(Dense(16, kernel_initializer='normal', activation='relu'))
    model.add(Dense(1, kernel_initializer='normal'))

    model.compile(loss='mean_squared_error', optimizer=optimizer)
    return model

model = KerasRegressor(build_fn=create_model, verbose=0)
# define the grid search parameters
optimizer = ['SGD','Adam']
batch_size = [10, 30, 50]
epochs = [10, 50, 100]
param_grid = dict(optimizer=optimizer, batch_size=batch_size, epochs=epochs)
reg_dl = GridSearchCV(estimator=model, param_grid=param_grid, n_jobs=-1)
reg_dl.fit(X_train, y_train)

Support Vector Regressor

# SVR
from sklearn.svm import SVR

reg_svr = GridSearchCV(SVR(kernel='rbf', gamma=0.1), cv=5,
                   param_grid={"C": [1e0, 1e1, 1e2, 1e3],
                               "gamma": np.logspace(-2, 2, 5)})
reg_svr.fit(X_train, y_train)

　これらの予測値を列にまとめ、新しい、第２の訓練マトリックスを作ります。

# second feature matrix
X_train2 = pd.DataFrame( {'XGB': reg_xgb.predict(X_train),
     'DL': reg_dl.predict(X_train).ravel(),
     'SVR': reg_svr.predict(X_train),
    })

　上5行だけ見てみると、どのモデルも比較的近い値を予測していますね。

　では、最も単純な線形モデルを使って、各モデルの重みを決定し、最終的な全体の予測値を計算しましょう。最後に、ターゲットをlogにしていたので、expを使って元のスケールに戻すことを忘れないでください。

# second-feature modeling using linear regression
from sklearn import linear_model

reg = linear_model.LinearRegression()
reg.fit(X_train2, y_train)

# prediction using the test set
X_test2 = pd.DataFrame( {'XGB': reg_xgb.predict(X_test),
     'DL': reg_dl.predict(X_test).ravel(),
     'SVR': reg_svr.predict(X_test),
    })

# Don't forget to convert the prediction back to non-log scale
y_pred = np.exp(reg.predict(X_test2))

# submission
submission = pd.DataFrame({
    "Id": test_ID,
    "SalePrice": y_pred
})
submission.to_csv('houseprice.csv', index=False)

　これでめでたく、提出用のCSVが作れました！

終わりに

　Kaggleでよく使うテクニックを用いて、Regressionの課題に取り組んでみました。これで立派なKagglerの仲間入りです（多分）。グリッドサーチを使ったとはいえ、各モデルのチューニングはまだ十分とは言えないので、それ次第でもっと上に行けると思います。アドバイスがあればよろしくお願いします^^