kaggleのHousePriceに挑戦してみた

#はじめに
2018年夏の記事です
私は初挑戦のコンペとしてkaggleのHousePricesに挑戦しました。
HousePricesはアメリカ、アイオワ州のエイムズにおける居住宅の販売価格を79個の特徴量から予測する回帰問題です。

#データの観察
##データの大きさの確認

#データの大きさを調べる
print(train.shape)
print(test.shape)

(1460,81)
(1459,80)

各データには、特徴量79個とId,SalePrice(trainのみ)のcolumnsがあり、1460個ほどのデータが揃っていることがわかる。

##特徴量
特徴量が膨大であるので、中でも重要なものと意味が分かりにくいものを下にあげます。

##目的関数の観察
目的関数の"SalePrice"をグラフ化してみる。

グラフからわかるように、住宅の販売価格が上がるにつれて住宅の個数は少なくなっている。
このままでは分析に影響を及ぼし来るので、対数をとって正規化する。

また、SalePriceの各統計量は以下のようになった。

count      1460.000000
mean     180921.195890
std       79442.502883
min       34900.000000
25%      129975.000000
50%      163000.000000
75%      214000.000000
max      755000.000000
Name: SalePrice, dtype: float64

##データの種類の確認
以下の作業はtrain,testデータ共に行うが、下では簡略のため、trainデータでの作業を示す。
object型の変数はアルゴリズムにかける前に変換する必要があるため、train,testデータともに、種類の確認を行います

#データの種類の確認
print(train.dtypes.head())

Id               int64
MSSubClass       int64
MSZoning        object
LotFrontage    float64
LotArea          int64
dtype: object

##欠損値の個数を確認

#trainデータの欠損値の確認
print(train.isnull().sum().head()[train.isnull().sum()>0])

LotFrontage     259
MasVnrType        8
MasVnrArea        8
BsmtQual         37
BsmtCond         37
BsmtExposure     38
BsmtFinType1     37
BsmtFinType2     38
Electrical        1
GarageType       81
GarageYrBlt      81
GarageFinish     81
GarageQual       81
GarageCond       81
dtype: int64

上記より、train,testデータ共に欠損値の多い、"Alley","PoolQC","Fence","MiscFeature","FireplaceQu"の項目を削除します。

#データの前処理

##欠損値の補完
地下室に関する項目(Bsmt関係)や車庫に関する項目(Garage関係)は欠損値の個数が同等の値となっているので、建物に、地下室やガレージがないものと考え"None"で補完します。
また、欠損値の少ない項目については、その項目の最頻値で補完します。

#欠損値の補完
#BasementとGarageの欠損値は存在しない場合が考えられるためNoneでReplace
train["BsmtQual"].fillna('None', inplace=True)
test["BsmtQual"].fillna('None', inplace=True)
train["BsmtCond"].fillna('None', inplace=True)
test["BsmtCond"].fillna('None', inplace=True)
train["BsmtExposure"].fillna('None', inplace=True)
test["BsmtExposure"].fillna('None', inplace=True)
train["BsmtFinType1"].fillna('None', inplace=True)
test["BsmtFinType1"].fillna('None', inplace=True)
train["BsmtFinSF1"].fillna('None', inplace=True)
test["BsmtFinSF1"].fillna('None', inplace=True)
train["BsmtFinType2"].fillna('None', inplace=True)
test["BsmtFinType2"].fillna('None', inplace=True)
train["BsmtFinSF2"].fillna('None', inplace=True)
test["BsmtFinSF2"].fillna('None', inplace=True)
train["BsmtUnfSF"].fillna('None', inplace=True)
test["BsmtUnfSF"].fillna('None', inplace=True)
train["TotalBsmtSF"].fillna('None', inplace=True)
test["TotalBsmtSF"].fillna(ans["TotalBsmtSF"].mean(), inplace=True)
train["GarageType"].fillna('None', inplace=True)
test["GarageType"].fillna('None', inplace=True)
#GarageYrBltの欠損値を平均値補完
train["GarageYrBlt"].fillna(train["GarageYrBlt"].mean(), inplace=True)
test["GarageYrBlt"].fillna(test["GarageYrBlt"].mean(), inplace=True)
train["GarageFinish"].fillna('None', inplace=True)
test["GarageFinish"].fillna('None', inplace=True)
train["GarageCars"].fillna('None', inplace=True)
test["GarageCars"].fillna(ans["GarageCars"].mean(), inplace=True)
train["GarageArea"].fillna('None', inplace=True)
test["GarageArea"].fillna(ans["GarageArea"].mean(), inplace=True)
train["GarageQual"].fillna('None', inplace=True)
test["GarageQual"].fillna('None', inplace=True)
train["GarageCond"].fillna('None', inplace=True)
test["GarageCond"].fillna('None', inplace=True)
train["MasVnrType"].fillna("None",inplace=True)
test["MasVnrType"].fillna("None",inplace=True)
train["MasVnrArea"].fillna(0,inplace=True)
test["MasVnrArea"].fillna(0,inplace=True)
test["Exterior1st"].fillna('None',inplace=True)
test["Exterior2nd"].fillna('None',inplace=True)
test["BsmtFullBath"].fillna('None',inplace=True)
test["BsmtHalfBath"].fillna('None',inplace=True)
test["KitchenQual"].fillna('None',inplace=True)
test["Functional"].fillna('None',inplace=True)
test["SaleType"].fillna('None',inplace=True)
#Electricalは最頻要素で補完
train["Electrical"].fillna('SBrkr',inplace=True)
test["MSZoning"].fillna('RL',inplace=True)
#Neighborhoodでグループし平均をとって補完
f = lambda x: x.fillna(x.mean())
train["LotFrontage"] = train.groupby("Neighborhood")["LotFrontage"].transform(f)
test["LotFrontage"] = test.groupby("Neighborhood")["LotFrontage"].transform(f)

改めて欠損値の確認をします。

#欠損値の確認
print(train.isnull().any().head())
print("\n")
print(test.isnull().any().head())

Id             False
MSSubClass     False
MSZoning       False
LotFrontage    False
LotArea        False
dtype: bool


Id             False
MSSubClass     False
MSZoning       False
LotFrontage    False
LotArea        False
dtype: bool

.head()を外して表示すると、欠損値が全て補完されたことがわかる。

##Feature Engineering
"1stFlrSF"と"2ndFlrSF"などのまとめられそうな特徴量はまとめます。

#FeatureEngineering
#総床面積を算出する
train["TotalFlrSF"]=train["1stFlrSF"]+train["2ndFlrSF"]
test["TotalFlrSF"]=test["1stFlrSF"]+test["2ndFlrSF"]
#総ベランダ面積を算出する
train["TotalHousePorchSF"]=train["EnclosedPorch"]+train["OpenPorchSF"]+train["WoodDeckSF"]+train["3SsnPorch"]+train["ScreenPorch"]
test["TotalHousePorchSF"]=test["EnclosedPorch"]+test["OpenPorchSF"]+test["WoodDeckSF"]+test["3SsnPorch"]+test["ScreenPorch"]                               
#建設年から築年数への変換
train["Age"]=2018-train["YearBuilt"]
test["Age"]=2018-test["YearBuilt"]
#リフォームされた年から、されてからの期間への変換
train["YrSinceRemod"]=2018-train["YearRemodAdd"]
test["YrSinceRemod"]=2018-test["YearRemodAdd"]                             
#売年から売れてからの年への変換
train["YrSinceSold"]= 2018-train["YrSold"]
test["YrSinceSold"]= 2018-test["YrSold"]
#garageに関しても同じ処理を施す
train["GarageSinceYrBlt"]=2018-train["GarageYrBlt"]
test["GarageSinceYrBlt"]=2018-test["GarageYrBlt"]

##LabelEncoderを使ってmapping
先述したように、object型の説明変数に関して、以下のone-hotエンコーディングでmappingします。

#RabelEncoderでmapping
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))

(参考:https://qiita.com/katsu1110/items/a1c3185fec39e5629bcb)

##各項目の重要度を調べる
さて、欠損値の補完と、object変数の変換ができたので、RandomForestRegressorを用いて各columnsの重要度を算出します。
まず、データを以下のように成形します

X_data=train.drop(["Id","SalePrice"],axis=1)
y_data=data["SalePrice"]

これらを用いてRandomForestRegressorにかけてみます。

#各項目の重要度を算出する
from sklearn.ensemble import RandomForestRegressor
rfr=RandomForestRegressor(random_state=0,n_estimators=100)
rfr.fit(X_data,y_data)

実装できたら、重要度を降順に並べ替えて、グラフを用いて重要度を可視化します。

ranking=np.argsort(-rfr.feature_importances_)
ax=plt.subplots(figsize=(10,12))
sns.barplot(x=rfr.feature_importances_[ranking],y=X_data.columns[ranking])
plt.show()

(参考:https://qiita.com/katsu1110/items/a1c3185fec39e5629bcb)

上のグラフからは、"OverallQual"が圧倒的に重要度が高いとわかります。
このグラフから重要度が高いcolumnsを上から15個取り出し,"X_pickup"とします。

##重要度の高い各columnsについて相関関係を調べる

#important_columnsとSalePriceとの相関をみる
import seaborn as sns
fig=plt.figure(figsize=(12,9))
for i in np.arange(14):
    ax=fig.add_subplot(5,3,i+1)
    sns.regplot(x=X_pickup.iloc[:,i],y=y_data)
    
plt.tight_layout()


plt.show()

(参考:https://qiita.com/katsu1110/items/a1c3185fec39e5629bcb)

以上から、相関関係と重要度が共に高い4つの説明変数"OverallQual","TotalFlrSF","TotalBsmtSF","GarageCars"と自分が不動産の価格決定の際に重要そうな"Age"を分析に用いる変数とします。

##外れ値の除去
上のregplotより、明らかな外れ値を除去します。

data_2=pd.concat([X_pickup,y_data],axis=1)
data_2=data_2.drop(data_2[(data_2["TotalFlrSF"]>4000)&(data_2["SalePrice"]<13)].index)
data_2=data_2.drop(data_2[(data_2["TotalBsmtSF"]>5000)&(data_2["SalePrice"]<14)].index)
data_2=data_2.drop(data_2[(data_2["GrLivArea"]>4000)&(data_2["SalePrice"]<13)].index)
data_2=data_2.drop(data_2[(data_2["GarageArea"]>1200)&(data_2["SalePrice"]<12)].index)                          
data_2=data_2.drop(data_2[(data_2["Age"]>120)&(data_2["SalePrice"]>12.5)].index)

#機械学習のアルゴリズムでデータを分析する
データの前処理が終わったので、いよいよデータをアルゴリズムにかけます。
今回は、SalePriceを予測する回帰分析であるため、一般的な重回帰分析と、Gradient BoostingとRandom Forestsを組み合わせたアンサンブル学習であるXGboostを用いました。

##1)重回帰分析にかける
多重共線性を加味して、説明変数は以下のように限定しました。

#特に重要な項目を抽出
X_1_pickup=data_2.loc[:,["OverallQual","TotalFlrSF","TotalBsmtSF","GarageArea","Age"]]
y_1_data=data_2.loc[:,["SalePrice"]]

データをtestデータとtrainデータに分割します（testデータの割合は3割）

from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test=train_test_split(X_1_pickup,y_1_data,test_size=0.3,random_state=0)

データの分割までできたので、いよいよtrainデータを分析にかけます。

from sklearn.linear_model import LinearRegression
lr=LinearRegression()
lr.fit(X_train,y_train)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
In [114]:

adjusted R^2
train score: 0.850568
test score : 0.840148

train,testデータともに調整済み決定係数を算出すると、85%程の精度で予測できていることがわかりました。

##2) XGboostにかける
次に、XGboostにデータをかけてみます。この際多くのパラメータの調整が必要であり、最適なパラメータはグリッドサーチを用いて求めました。

import xgboost as xgb
from sklearn.model_selection import GridSearchCV

xgb_model = xgb.XGBRegressor()
xgb_gs = GridSearchCV(xgb_model,
                   {'max_depth': [1,2,3],
                    'n_estimators': [80,90,100]})

xgb_gs.fit(X_train, y_train)

こちらに関しても調整済み決定係数を出力すると以下のようになりました。

adjusted R^2
train score: 0.883398
test score : 0.843075

若干過学習を起こしいますが、testデータの精度も85%近く出ています。

今回はこの二つの手法で得られた予測値の平均値を用いてkaggleに提出しました。

y_pred=np.exp(lr.predict(ans_data).reshape(-1))
y_2_pred=np.exp(reg_xgb.predict(ans_data).reshape(-1))
y_stack=(y_pred)*(1/2)+(y_2_pred)*(1/2)

#結果/今後の展望
最終的なスコアは、0.15741 (2947位/4187人中) でした。
今回は、決定係数による精度が85%ほどでした。
今回のデータ分析では、LotFrontage(道から家屋までの距離)などの不動産の価格決定ではあまり必要ないと思っていた特徴量が重要であることが分かりました。

また、今回は、特徴量を５つ抽出して分析しましたが、別の特徴量を用いて分析して見るとまた別の結果が得られると思うので、今後試して見たいです。

アルゴリズムに、RandomForestを選択して分析してみましたが決定係数が test:96%,train:82%とかなり過学習してしまったのでスタッキングの対象からは外しました。今後時間があるときにランダムフォレストでの過学習を防ぐ施策を試して見たいです。
今回は、データ分析の基礎の流れを実際に体感できる良い経験になりました。今後さらに色々な分析に挑戦していきます。