半歩ずつ進める機械学習 ~Kaggleに挑戦 HousePrices③~

##前回まで
前回の投稿ではKaggleのHousePricesコンペに参加して
とりあえず最低限の労力でモデルを作り、4700人中4000位というnoobらしいスコアを叩きだしました
実際にKaggleに挑戦して改めてよく分かったのが前処理などの難しさ・大変さです
前回はすでに数値のデータだけを用いて回帰を行いましたが
今回はオブジェクトデータの値型データへの変換と欠損値の補完を行い
とりあえず与えられたデータを全て使ってモデルを作り、スコアを確認してみます。
次回以降は、さらに外れ値の除去や特徴量エンジニアリングを行っていこうと思います。

##データの準備

# This Python 3 environment comes with many helpful analytics libraries installed
# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python
# For example, here's several helpful packages to load in 

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

# Input data files are available in the "../input/" directory.
# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory

import os
print(os.listdir("../input"))

from sklearn.preprocessing import StandardScaler
import pandas_profiling as pdp
# Any results you write to the current directory are saved as output.

train = pd.read_csv("../input/train.csv")
test = pd.read_csv("../input/test.csv")
y = train.SalePrice
train = train.drop("SalePrice",axis = 1)
allData = pd.concat([train,test])
allData = allData.drop("Id",axis = 1)

ここまでは前回と同じです、訓練データと正解値となる教師データ
それと実際に予測を行うテストデータをロードします。

欠損値は訓練データ・テストデータ両方に存在するので、一気に処理するため
一旦訓練データとテストデータを結合します
また、全体のデータ分布等がどうなっているかの確認をしてみます。

allData = pd.concat([train,test])

allData.info()
#pdp.ProfileReport(allData)

データ内容

<class 'pandas.core.frame.DataFrame'>
Int64Index: 2919 entries, 0 to 1458
Data columns (total 80 columns):
MSSubClass       2919 non-null int64
MSZoning         2915 non-null object
LotFrontage      2433 non-null float64
LotArea          2919 non-null int64
Street           2919 non-null object
Alley            198 non-null object
LotShape         2919 non-null object
LandContour      2919 non-null object
Utilities        2917 non-null object
LotConfig        2919 non-null object
LandSlope        2919 non-null object
Neighborhood     2919 non-null object
Condition1       2919 non-null object
Condition2       2919 non-null object
BldgType         2919 non-null object
HouseStyle       2919 non-null object
OverallQual      2919 non-null int64
OverallCond      2919 non-null int64
YearBuilt        2919 non-null int64
YearRemodAdd     2919 non-null int64
RoofStyle        2919 non-null object
RoofMatl         2919 non-null object
Exterior1st      2918 non-null object
Exterior2nd      2918 non-null object
MasVnrType       2895 non-null object
MasVnrArea       2896 non-null float64
ExterQual        2919 non-null object
ExterCond        2919 non-null object
Foundation       2919 non-null object
BsmtQual         2838 non-null object
BsmtCond         2837 non-null object
BsmtExposure     2837 non-null object
BsmtFinType1     2840 non-null object
BsmtFinSF1       2918 non-null float64
BsmtFinType2     2839 non-null object
BsmtFinSF2       2918 non-null float64
BsmtUnfSF        2918 non-null float64
TotalBsmtSF      2918 non-null float64
Heating          2919 non-null object
HeatingQC        2919 non-null object
CentralAir       2919 non-null object
Electrical       2918 non-null object
1stFlrSF         2919 non-null int64
2ndFlrSF         2919 non-null int64
LowQualFinSF     2919 non-null int64
GrLivArea        2919 non-null int64
BsmtFullBath     2917 non-null float64
BsmtHalfBath     2917 non-null float64
FullBath         2919 non-null int64
HalfBath         2919 non-null int64
BedroomAbvGr     2919 non-null int64
KitchenAbvGr     2919 non-null int64
KitchenQual      2918 non-null object
TotRmsAbvGrd     2919 non-null int64
Functional       2917 non-null object
Fireplaces       2919 non-null int64
FireplaceQu      1499 non-null object
GarageType       2762 non-null object
GarageYrBlt      2760 non-null float64
GarageFinish     2760 non-null object
GarageCars       2918 non-null float64
GarageArea       2918 non-null float64
GarageQual       2760 non-null object
GarageCond       2760 non-null object
PavedDrive       2919 non-null object
WoodDeckSF       2919 non-null int64
OpenPorchSF      2919 non-null int64
EnclosedPorch    2919 non-null int64
3SsnPorch        2919 non-null int64
ScreenPorch      2919 non-null int64
PoolArea         2919 non-null int64
PoolQC           10 non-null object
Fence            571 non-null object
MiscFeature      105 non-null object
MiscVal          2919 non-null int64
MoSold           2919 non-null int64
YrSold           2919 non-null int64
SaleType         2918 non-null object
SaleCondition    2919 non-null object
IsTrainData      2919 non-null int64
dtypes: float64(11), int64(26), object(43)
memory usage: 1.8+ MB

コメントアウトしていますが、pandas_profiling.ProfileReportを使用すれば
各特徴量の分布等をヒストグラムで確認する事が出来るので
データの分布を確認するのに非常に便利でした。

とりあえず、データ件数が全部で2919件あり、欠損値を含む特徴量が複数ある事が分かりました。
まずは欠損値を埋めていきましょう

※欠損の数があまりに多い場合、無理にデータを埋めるより
特徴量自体を削除してしまうケースもあります。(結局信頼出来ないので）
どれくらい欠損していたら削除するか、といった統計的な定石は見つける事が出来ませんでした
特徴量の重要度等を考慮して、手探りで行うものだと思います。

##欠損したデータか、欠損というデータか
欠損値は、入力作業者のミス等により純粋に欠損している物と
欠損している事自体が何かを表現したデータであるケースがあります。
具体的に言うと、PoolQCという特徴量は2919件中2909件欠損していますが
Kaggleから与えられているdata_description.txtを確認するとPoolQCというデータの説明は以下の様になっています。

PoolQC: Pool quality

   Ex   Excellent
   Gd   Good
   TA   Average/Typical
   Fa   Fair
   NA   No Pool

NA = No Poolという事で、欠損値はプールが無い事を表現しているようです。
プールなんて高級物件にしかなさそうですし、2919件中2909件にプールがないと言われれば、そうだよねって感じです。
こういったケースではallData["PoolQC"].fillna("None")とかで単純に適当なカテゴリとして変換します。
PoolQCと同様に欠損自体が何かを表現している特徴量は他にもあるので、それらも同様に処理します。

##データが純粋に欠損している場合
次に純粋に欠損しているケースを考えます。
このケースでも大まかに二つの対応が考えられます。

まず第一に単純にデータの平均等で置き換える方法です
数値データであれば、欠損していないデータの平均や中央値を
カテゴリカルデータであれば最頻出のデータで置き換えるといった方法が考えられます。

第二の方法は他のデータや、ドメイン知識から推測するという方法です。
例えばLotFrontageというデータは500件近い欠損がある数値データで、データの内容は物件に面している通りの直線距離です。
数値データなので、単純に平均しても良いのですが、ある公開カーネルでなるほどな～と思った埋め方がありました。
それは、道路の長さには地域制があり、また同じ地域の物件は同じ通りに面している可能性があるので
地域毎のLotFrontageの平均で埋める、という方法でした。
また、他にも不動産業界に勤めている人や、その地域に住んでいる人は、その分野の知識（ドメイン知識）に基づいて
単純な平均よりも現実に近い推測を行う事が出来るかもしれません。

##とりあえず埋める
何はともあれ埋めていきます。
コードが冗長的かもしれませんが、Python初心者なので御容赦下さい。

#地域毎のLotFrontage平均で埋める。
allData.LotFrontage[allData.LotFrontage.isnull()] = allData.groupby("Neighborhood").LotFrontage.mean()[allData.loc[allData.LotFrontage.isnull(),"Neighborhood"]]

#Basement
allData.BsmtQual[allData["BsmtFinType1"].isnull()] = "None"
allData.BsmtCond[allData["BsmtFinType1"].isnull()] = "None"
allData.BsmtExposure = allData.BsmtExposure.fillna("None")
allData.BsmtFinType1[allData.BsmtFinType1.isnull()] = "None"
allData.BsmtFinType2[allData.BsmtFinType2.isnull()] = "None"

#最頻出のカテゴリで埋める。
allData.Electrical = allData.Electrical.fillna("SBrkr")

allData.Alley[allData.Alley.isnull()] = "None"

allData.FireplaceQu[allData.Fireplaces == 0] = "None"

allData.PoolQC[allData.PoolArea == 0] = "None"
allData.PoolQC = allData.PoolQC.fillna("TA")

allData.Fence[allData.Fence.isnull()] = "None"

allData.SaleType = allData.SaleType.fillna("WD")

allData.KitchenQual = allData.KitchenQual.fillna("TA")

allData.MiscFeature[allData.MiscVal == 0] = "None"
allData.MiscFeature = allData.MiscFeature.fillna("Gar2")

allData.Functional = allData.Functional.fillna("Typ")

allData.Exterior1st = allData.Exterior1st.fillna(modeVal["Exterior1st"])
allData.Exterior2nd = allData.Exterior2nd.fillna(modeVal["Exterior2nd"])

#Garage
allData.GarageFinish = allData.GarageFinish.fillna("None")
allData.GarageType[allData.GarageArea == 0] = "None"
allData.GarageQual[allData.GarageArea == 0] = "None"
allData.GarageCond[allData.GarageArea == 0] = "None"
allData.GarageQual = allData.GarageQual.fillna("TA")
allData.GarageCond = allData.GarageCond.fillna("TA")
allData.loc[allData.GarageType == "None","GarageYrBlt"] = 0
allData.GarageYrBlt = allData.GarageYrBlt.fillna(allData[allData.GarageType != "None"].GarageYrBlt.mean())
allData.GarageCars = allData.GarageCars.fillna(allData[allData.GarageType != "None"].GarageCars.median())
allData.GarageArea = allData.GarageArea.fillna(allData[allData.GarageType != "None"].GarageArea.mean())


allData.MasVnrType = allData.MasVnrType.fillna("None")
allData.MasVnrArea = allData.MasVnrArea.fillna(0)

allData.Utilities = allData.Utilities.fillna("AllPub")
allData.MSZoning = allData.MSZoning.fillna("RL")

allData.BsmtFullBath = allData.BsmtFullBath.fillna(0)
allData.BsmtHalfBath = allData.BsmtHalfBath.fillna(0)

allData.BsmtCond = allData.BsmtCond.fillna("None")
allData.BsmtQual = allData.BsmtQual.fillna("None")
allData.BsmtFinSF1 = allData.BsmtFinSF1.fillna(0)
allData.BsmtFinSF2 = allData.BsmtFinSF2.fillna(0)
allData.BsmtUnfSF = allData.BsmtUnfSF.fillna(0)
allData.TotalBsmtSF = (allData.BsmtFinSF1+allData.BsmtFinSF2+allData.BsmtUnfSF)

このデータセットは
【暖炉の有無】【暖炉の品質】
【ガレージの有無】【ガレージの品質】
といった感じで設備の有無と、その品質や規模の説明を行う特徴量が多かったです。
なので、例えば【暖炉の品質】が欠損している場合
同じデータの【暖炉の有無】を確認して、暖炉が無かった場合はNoneで埋める
暖炉が存在した場合は最頻出で埋める、という具合に他のデータから推測出来る場合はそのようにして
そうでない場合は平均や最頻出のデータで埋めました。

欠損値の補完を行う前に、特徴量の重要性を確認して
重要性の低い場合は予め削除しておけば、こういった作業が軽減されるかと思いますが
今回は、ひとまず愚直に全ての欠損値を埋めました。

##カテゴリカルデータのエンコード
このデータには多数のカテゴリカルデータが存在するので、それらをOne-hot表現に変換します。
allDataからobject型のデータを引き抜いて処理を行いますが
この時、MSZoningだけは値型ですが、名義尺度なのでobject型データ同様one-hot表現に変換します。

#オブジェクト型の特徴量のリストを作成
objList = allData.select_dtypes(include = object).columns.values.tolist()
objList.append("MSZoning")

import category_encoders as ce
#引数にカラム名を与える事で、該当するカラム名のデータをone-hot表現にエンコードしてくれる。
ce_bin = ce.BinaryEncoder(cols = objList,handle_unknown = "impute")
trans_allData = ce_bin.fit_transform(allData)

print(trans_allData.shape)
#(2919, 205)

元々79個あった特徴量が205個に増えました。
データによってはone-hotエンコードは特徴量を爆発的に増やす可能性があるので注意が必要です。

##結合したデータを分割する
これでようやくデータが全て値型になりました。
モデル作りの前に、データの標準化を行います。
またallDataは訓練データとテストデータが一緒になっているので分割します。

sc = StandardScaler()
allData_z = sc.fit_transform(trans_allData)

train_z = allData_z[:1460,:]
test_z = allData_z[1460:,:]

##教師データを対数変換する
公開カーネルを見て知った方法です。
このデータの教師データの分布は以下のようになっています。

ヒストグラム

QQプロット

QQプロットの赤線に一致していればいるほど、正規分布に近い分布であるという事なんですが
これを見ると正規分布に近いとは言い難いようです。

線形回帰はもとになるデータが正規分布である事を前提としています。
なので、回帰モデルを作る前に、データの分布を正規分布に近づけたいと思います。

また、このコンペでのモデルの評価基準は**予測値を対数変換したRMSE（平均二乗平方誤差）**です

Metric
Submissions are evaluated on Root-Mean-Squared-Error (RMSE) between the logarithm of the predicted value
and the logarithm of the observed sales price. (Taking logs means that errors in predicting expensive houses and cheap houses will affect the result equally.)

メトリック
提出は、予測値の対数と観察された販売価格の対数の間の二乗平均平方根誤差（RMSE）で評価されます。 （ログをとることは、高価な家屋と安い家屋を予測する際の誤りが結果に等しく影響することを意味します。）

という事で、予測値を対数変換してからRMSEで評価するなら、モデルも対数変換した教師データを使ってクロスバリデーションを行い
最もRMSEが小さくなるモデルにしたいところです。

※100%評価尺度と同じ尺度でモデルを最適化するべきというワケではありません。
他のコンペで評価尺度と異なる損失関数を最小化したモデルが上位に来たケースもあるようです。

とりあえず教師データを対数変換して標準化してみましょう。

log_y = np.log(y)
train_yz = sc.fit_transform(log_y.values.reshape(-1,1)) 

ここで、対数変換した教師データの分布を再度確認してみます。

ヒストグラム

QQプロット

かなり正規分布に近づいたように見えます、これは回帰分析を行う上でも都合が良いです。

##予測モデルを作る
ようやく予測モデルを作ります。
特徴量が多いので、多重共線性が発生しないように、Ridgeでモデルを作って交差検証を行います。
検証のスコアはRMSEにしたいのですが、scikit-learnのcross_validateではRMSEが無いっぽいので自前で作ります。

from sklearn.linear_model import Ridge
from sklearn.model_selection import KFold,cross_validate
from sklearn.metrics import make_scorer
import pprint

def RMSE(y_true,y_pred):
    #標準化された正解ラベルと予測データが渡されるので元に戻す
    y_true_in = sc.inverse_transform(y_true)
    y_pred_in = sc.inverse_transform(y_pred)
    return np.sqrt(((y_true_in - y_pred_in) ** 2).mean())


rg = Ridge()
rg.fit(train_z,train_yz)

kf = KFold(n_splits = 5,shuffle = True,random_state = 1)
splitter = kf.split(train_z,train_yz)

pprint.pprint(cross_validate(rg,train_z,train_yz,cv = splitter,scoring = make_scorer(RMSE)))

{'fit_time': array([0.01935649, 0.02468395, 0.0262053 , 0.02337122, 0.02046204]),
 'score_time': array([0.00635767, 0.00594354, 0.00853395, 0.00335193, 0.00622702]),
 'test_score': array([0.14089677, 0.23849176, 0.12243584, 0.18486678, 0.11186544]),
 'train_score': array([0.10871249, 0.10305512, 0.11247659, 0.10262155, 0.11475878])}

前回のスコアが0.21394だったので、良くなってる･･･るかな？
実際にサブミットして確認してみます。

交差検証では0.23とか低いスコアも出てましたが、サブミットしてみるとかなり良い結果が出ました。
公開カーネル様のお蔭なので、自分の実力とは言い難いですが、ここから更に

• 外れ値の除去
• 特徴量エンジニアリング
• ハイパーパラメータのチューニング
• 複数モデルによるアンサンブル

等を行う事で、もっとスコアが高くなるかも！と思うとワクワクしますね

##次回へ
次回は外れ値の除去と特徴量エンジニアリングの手法について学びたいと思います