はじめに
大変お世話になったページ
Kaggleの住宅価格予測問題(House Prices)を解いてみる
Kaggle House Prices -Qiita
こんにちは、新卒2年目、データエンジニア兼データサイエンティスト(の見習い)をしている者です。
現在、データ分析及び機械学習の経験を積みたくて、新人研修で「Titanic: Machine Learning from Disaster」に挑戦して以来、約1年ぶりに「Kaggle」に挑戦しています。
まず手始めに回帰問題版Titanicとも言われる、「House Prices: Advanced Regression Techniques」にトライしてみることにしました。
こんな人向けの記事です
- KaggleでTitanicはやってみたけど、そのあと何したらいいかわかんないよ!
- 分類問題はやったことあるけど、回帰問題は初めてやるよ!
- House Pricesに今まさにトライしてみるところだよ!
こんな人には少し説明不足な記事かも
- Kaggle自体をはじめてやるよ!
(Kaggle自体の使い方は割愛しています。) - 表データに対する機械学習予測を初めてやるよ!
(わりと丁寧には書いたつもりですが、欠損値など超基礎的な内容は割愛しています。)
Kaggleの攻略法はパクることと見つけたり
人は誰しもはじめは初心者です。
何事も人の型を真似するところから始まります。
一方でKaggleの場合はなぜか「自分なりのアイデアで分析をして結果を出すぞ!」と意気込んで挫折する人が多いイメージです。(僕ももともとはそのうちの一人です。)
しかし、それはとても愚かなことだと気付きました。
数学の勉強を始めようというときに、先人が幾千年もかけて解き明かした公式や解き方が教科書に載っているにも関わらず、「いや、俺は自力で数学の謎を解き明かすから!」と、教科書を使わずに学び始めるくらい愚かなことです。
先人の知恵や経験は、自分が近道を辿るための「ツール」です。
幸いなことにKaggleには「Notebook」や「Discussion」というKagglerの巣窟とも呼べるような場があり、パクリ元の宝庫になっています。
また、今回の「House Prices」のような初心者向けの課題であれば、日本語で解説をしたページもいくらか出てきます。
単純にスコアをあげるための「攻略法」というだけではなく、先人の型を真似て吸収して、同じ土俵に立つための「攻略法」としてパクリをガンガン許容していくべきだと思います。
もちろん他人のNotebookをまるまるパクって提出して、「やったー、スコアが上がった!」というだけではダメで、何がスコアを上げる要因になったか学ぶ必要があるということですね。
House Pricesの概要
さて、精一杯パクリを正当化したところで本題に入ります。
「House Prices: Advanced Regression Techniques」はアメリカ合州国 アイオワ州 エイムズ市の戸建て住宅の価格を79個の変数から予測する、テーブルデータの回帰予測問題です。
学習データ : 1460件
予測データ : 1459件
スコアの評価方法 : RMSE(二乗平均平方根誤差) すごく雑にいうと値が小さいほどハイスコア
参考:いろいろな誤差の意味 -具体例で示す数学
予備情報
アイオワ州 - Wikipedia
めっちゃ余談ですが、僕の大好きなSlipknotのほとんどのメンバーはアイオワ州出身です。
そんな彼らのセカンドアルバムのタイトルは「Iowa」だったりします。
アイオワ州 エイムズ市 - Wikipedia
住宅の価格を予測する問題なので、「敷地面積とか築年数は住宅価格と関係ありそうだな!」といった予想がたてられるくらいにはイメージしやすい問題だと思います。
予測に使う変数の簡単な説明は実際のデータを読み込んだ後にやるものとして、まずは準備を進めていきます。
ちなみに僕のNotebookはこちらで公開しています。(この記事とほとんど同じ内容ですが、よかったら投票をお願いいたします。)
#目次
1 ライブラリインポート
2 学習データの読み込み
3 データ前処理
4 学習モデルの作成
5 住宅価格を予測する
6 提出用csvを作成し、提出
#1 ライブラリインポート
# numpy , pandas
import numpy as np
import pandas as pd
# scikit-learn
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Lasso
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
# 可視化用ライブラリ
import matplotlib.pyplot as plt
import seaborn as sns
#pandasのカラムが100列まで見れるようにする
pd.set_option('display.max_columns', 100)
#2 学習データの読み込み
# 学習データの読み込み
train_df = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/train.csv',index_col=0)
# 先頭5行をみてみる。
train_df.head()
変数について
データの詳細な定義についてはHouse Prices / Dataに記載されています。
目的変数(本課題で予測する値) - 1個
- SalePrice : 住宅価格
説明変数(予測に使う変数) - 79個 一部抜粋
住宅価格に大きく影響しそうな変数
- LotArea : 敷地面積
- YearBuilt : 建築年
- YearRemodAdd : リフォーム年
住宅価格にある程度影響しそうな変数
- KitchenQual : キッチンの状態
- RoofMatl : 屋根の材質
めっちゃアメリカっぽいやつも。。。
- HeatingQC : 暖炉の状態
- PoolArea : プールの広さ
2.1 住宅価格について見てみる
目的変数である住宅価格(SalePrice)について見ていきましょう。
# 売却価格のヒストグラム
sns.distplot(train_df['SalePrice'])
# 売却価格の概要をみてみる
print(train_df["SalePrice"].describe())
print(f"歪度: {round(train_df['SalePrice'].skew(),4)}" )
print(f"尖度: {round(train_df['SalePrice'].kurt(),4)}" )
価格の安い物件が大多数を占めていそう。。。
大多数の一般的市民とごく少数の富裕層の物件があるイメージですね。
このような正規分布に従わないデータの処理については次章で説明します。
ちょっとだけ補足
分布の正規分布との形状の近さをみるときは、歪度と尖度という指標を使います。
歪度(わいど)
分布を正規分布と比較したときの偏り度合い(非対称性)を示す尺度。0に近いほど左右対称。
尖度(せんど)
分布を正規分布と比較したときの尖り具合を示す尺度。0に近いほど正規分布に近い。
また、正の大きい値をとるほど、ピーク付近に標本が集中しており、負の大きい値をとるほど、分布が散らばっている。
#3 データ前処理
データ前処理は機械学習のモジュールに合わせたデータ形式にしたり、予測の精度を上げるためにデータを整形する行程です。
今回、学習用データセットと予測用データセットを合わせて同時に前処理を行います。
(ただし、予測データは全ての行を残す必要があるため、外れ値の除去はしない。)
# 予測用データセットの読み込み
test_df = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/test.csv',index_col=0)
# 学習データの説明変数と、予測用データを結合
all_df = pd.concat([train_df.drop(columns='SalePrice'),test_df])
3.1 一部の数字が入っている変数を文字列に変換
数字が入っているが、数字の大小関係が予測に影響を与えない方が良いものはカテゴリ変数にすべきです。
数字をそのまま文字列に変換してしまいましょう。
例:血液型の変数が
{'A' : 0 ,'B' : 1 , 'O' : 2 , 'AB' : 3 }と割り当てられている場合、数の大小に意味はないので
{'A' : '0' ,'B' : '1' , 'O' : '2' , 'AB' : '3' }に変換する。
num2str_list = ['MSSubClass','YrSold','MoSold']
for column in num2str_list:
all_df[column] = all_df[column].astype(str)
3.2 欠損値の処理
欠損値の処理には、主に以下の2つがあります。
- 欠損値を含む行または列の削除
- 欠損値を代表値などで埋める
今回のデータにおいて欠損値となっている箇所は、その設備がないという場合が多いようです。
(例:PoolAreaが欠損している物件はそもそもプールがない。)
そこで今回は、文字列の変数の欠損は「'None'」、数字の変数の欠損は「0」で埋めるものとします。
# 変数の型ごとに欠損値の扱いが異なるため、変数ごとに処理
for column in all_df.columns:
# dtypeがobjectの場合、文字列の変数
if all_df[column].dtype=='O':
all_df[column] = all_df[column].fillna('None')
# dtypeがint , floatの場合、数字の変数
else:
all_df[column] = all_df[column].fillna(0)
3.3 特徴量エンジニアリング
特徴量エンジニアリングでは予測の精度向上に貢献しそうな新しい変数を追加していきます。
例:建物内の総面積 = 1階の面積 + 2階の面積 + 地下の面積
ここは自分のアイデアが入り込む余地がある大いにある処理ですね。
まあ、一旦パクりましたけど。。。
# 特徴量エンジニアリングによりカラムを追加する関数
def add_new_columns(df):
# 建物内の総面積 = 1階の面積 + 2階の面積 + 地下の面積
df["TotalSF"] = df["1stFlrSF"] + df["2ndFlrSF"] + df["TotalBsmtSF"]
# 一部屋あたりの平均面積 = 建物の総面積 / 部屋数
df['AreaPerRoom'] = df['TotalSF']/df['TotRmsAbvGrd']
# 築年数 + 最新リフォーム年 : この値が大きいほど値段が高くなりそう
df['YearBuiltPlusRemod']=df['YearBuilt']+df['YearRemodAdd']
# お風呂の総面積
# Full bath : 浴槽、シャワー、洗面台、便器全てが備わったバスルーム
# Half bath : 洗面台、便器が備わった部屋)(シャワールームがある場合もある)
# シャワーがない場合を想定してHalf Bathには0.5の係数をつける
df['TotalBathrooms'] = (df['FullBath'] + (0.5 * df['HalfBath']) + df['BsmtFullBath'] + (0.5 * df['BsmtHalfBath']))
# 合計の屋根付きの玄関の総面積
# Porch : 屋根付きの玄関 日本風にいうと縁側
df['TotalPorchSF'] = (df['OpenPorchSF'] + df['3SsnPorch'] + df['EnclosedPorch'] + df['ScreenPorch'] + df['WoodDeckSF'])
# プールの有無
df['HasPool'] = df['PoolArea'].apply(lambda x: 1 if x > 0 else 0)
# 2階の有無
df['Has2ndFloor'] = df['2ndFlrSF'].apply(lambda x: 1 if x > 0 else 0)
# ガレージの有無
df['HasGarage'] = df['GarageArea'].apply(lambda x: 1 if x > 0 else 0)
# 地下室の有無
df['HasBsmt'] = df['TotalBsmtSF'].apply(lambda x: 1 if x > 0 else 0)
# 暖炉の有無
df['HasFireplace'] = df['Fireplaces'].apply(lambda x: 1 if x > 0 else 0)
# カラムを追加
add_new_columns(all_df)
3.4 文字列のカテゴリ変数化(One-Hot-Encoding)
文字列が入っている変数はそのままでは、正規化したり機械学習にかけられないので、One-Hot-Encodingによりカテゴリ変数化します。
出典:Building a One Hot Encoding Layer with TensorFlow
# pd.get_dummiesを使うとカテゴリ変数化できる。
all_df = pd.get_dummies(all_df)
all_df.head()
元の変数名_値 の形式で列が生成されました。
これにより、説明変数の数が89個から350個に増え、文字列が入っていた列がカテゴリ変数化されました。
ちょっとだけ補足
カテゴリ変数に大小の意味を持たせたい場合はOrdinal Encoding(Label Encoding)を用います。
例: テストの結果に基づいた評定 { '90点以上' :'S', '80 ~ 89点' :'A', '70 ~ 79点' :'B', '70点未満' : 'C' } の時は大小関係に意味がある。
参考:カテゴリカル変数はなんでもダミー変換すればよいのか?-アルゴリズムに応じたOne Hot EncodingとLabel Encodingの使い分け
One-Hot-Encodingによって、変数が過剰に増えて、計算リソースが増えてしまうことを防ぐために、Ordinal Encodingを用いることもあるようです。
3.5 外れ値の除去
学習データから特異な物件のデータは削除します。
例:一般の住宅のデータにハリウッドスターの物件のデータが入っていると、そのデータに予測が影響を受けてしまう。
(予測データの物件は削除できないので、学習データに対してのみ実施)
# 学習データと予測データに分割して元のデータフレームに戻す。
train_df = pd.merge(all_df.iloc[train_df.index[0]:train_df.index[-1]],train_df['SalePrice'],left_index=True,right_index=True)
test_df = all_df.iloc[train_df.index[-1]:]
以下の条件に当てはまるものを外れ値として定義し、該当する物件のデータを削除します。
- 物件の価格が400,000ドル以上
- 敷地面積が20,000平方メートル以上
- 建築年が1920年より前
train_df = train_df[(train_df['LotArea']<20000) & (train_df['SalePrice']<400000)& (train_df['YearBuilt']>1920)]
##3.6 住宅価格を対数変換
多くの機械学習アルゴリズムは正規分布のデータを想定しているため、正規分布ではないデータに対して精度が出ない場合が多いです。
そのとき、正規分布に従わない変数の対数変換をすると正規分布に近付くことがあります。
(対数変換を行うことで正規分布になるような分布を、対数正規分布と呼ぶ。)
# 対数変換前のヒストグラム、歪度、尖度
sns.distplot(train_df['SalePrice'])
print(f"歪度: {round(train_df['SalePrice'].skew(),4)}" )
print(f"尖度: {round(train_df['SalePrice'].kurt(),4)}" )
# SalePriceLogに対数変換した値を入れる。説明の都合上新たなカラムを作るが、基本的にそのまま代入して良い。
# np.log()は底がeの対数変換を行う。
train_df['SalePriceLog'] = np.log(train_df['SalePrice'])
# 対数変換後のヒストグラム、歪度、尖度
sns.distplot(train_df['SalePriceLog'])
print(f"歪度: {round(train_df['SalePriceLog'].skew(),4)}" )
print(f"尖度: {round(train_df['SalePriceLog'].kurt(),4)}" )
ぱっと見わかりづらいですが、歪度、尖度をみるとかなり正規分布に近づいたことがわかりますね。
ちょっとだけ補足
元の目的変数の分布によっては対数変換が意味をなさない、もしくは悪影響を与える場合があるようです。
参考:対数変換が適さない場合がある!?対数変換すると結果が悪くなる例の紹介
3.7 学習データの説明変数と目的変数、予測データの説明変数にデータフレームを分割する。
# 学習データ、説明変数
train_X = train_df.drop(columns = ['SalePrice','SalePriceLog'])
# 学習データ、目的変数
train_y = train_df['SalePriceLog']
# 予測データ、目的変数
test_X = test_df
ここまででデータの前処理は終了です。
#4 学習モデルの作成
学習データの傾向を学習し、モデルを作成します。
今回は機械学習アルゴリズムとしてラッソ回帰(Lasso Regression)を用いる。
ラッソ回帰 - Wikipedia
ラッソ回帰を選んだ理由
- 変数の数が多く、少数の変数が予測に大きな影響を与えそうな場合に有効 (敷地面積とか築年数が影響を与えそう。。。)
- ハイパーパラメータがalphaだけ調整すれば行けそうな雰囲気だった。(チューニングがチョロそう)
4.1 ハイパーパラメータチューニング
交差検証(Grid Search)を用いて、ハイパーパラメータのalphaをチューニングします。
lasso_tuning -> 学習データ内で交差検証をして、最も精度の高かったalphaのパラメータを返す。
def lasso_tuning(train_x,train_y):
# alphaパラメータのリスト
param_list = [0.001, 0.01, 0.1, 1.0, 10.0,100.0,1000.0]
for cnt,alpha in enumerate(param_list):
# パラメータを設定したラッソ回帰モデル
lasso = Lasso(alpha=alpha)
# パイプライン生成
pipeline = make_pipeline(StandardScaler(), lasso)
# 学習データ内でホールドアウト検証のために分割 テストデータの割合は0.3 seed値を0に固定
X_train, X_test, y_train, y_test = train_test_split(train_x, train_y, test_size=0.3, random_state=0)
# 学習
pipeline.fit(X_train,y_train)
# RMSE(平均誤差)を計算
train_rmse = np.sqrt(mean_squared_error(y_train, pipeline.predict(X_train)))
test_rmse = np.sqrt(mean_squared_error(y_test, pipeline.predict(X_test)))
# ベストパラメータを更新
if cnt == 0:
best_score = test_rmse
best_param = alpha
elif best_score > test_rmse:
best_score = test_rmse
best_param = alpha
# ベストパラメータのalphaと、そのときのMSEを出力
print('alpha : ' + str(best_param))
print('test score is : ' +str(round(best_score,4)))
# ベストパラメータを返却
return best_param
# best_alphaにベストパラメータのalphaが渡される。
best_alpha = lasso_tuning(train_X,train_y)
最も精度が出たalphaは「0.01」、RMSE(平均二乗平方根誤差)は「0.0985」となりました。
4.2 ベストパラメータを用いてモデル作成
交差検証の結果をもとに、最優秀alphaをパラメータに設定して本番の学習をしてきます。
# ラッソ回帰モデルにベストパラメータを設定
lasso = Lasso(alpha = best_alpha)
# パイプラインの作成
pipeline = make_pipeline(StandardScaler(), lasso)
# 学習
pipeline.fit(train_X,train_y)
これで学習は完了です。scikit-learnはシンプルに記述できていいですね。
#5 住宅価格を予測する
それではいよいよ住宅価格を予測していきます。
5.1 作成したモデルで住宅価格を予測する。
# 結果を予測
pred = pipeline.predict(test_X)
拍子抜けかもしれませんがたった1行で終了です。scikit-learnはシンプ(ry
5.2 予測結果を指数変換
# 予測結果のプロット
sns.distplot(pred)
# 歪度と尖度
print(f"歪度: {round(pd.Series(pred).skew(),4)}" )
print(f"尖度: {round(pd.Series(pred).kurt(),4)}" )
おお、正規分布っぽい!
ですが、「3.6 住宅価格を対数変換」で対数変換した値を目的変数としてモデルを構築したので、予測された値も対数変換された値になっています。
そこで、予測した値を指数変換してあげる必要があります。
先ほどは底がeの対数で変換したので、np.exp()をとって戻してあげましょう。
# 指数変換
pred_exp = np.exp(pred)
# 指数変換した予測結果をプロット
sns.distplot(pred_exp)
# 歪度と尖度
print(f"歪度: {round(pd.Series(pred_exp).skew(),4)}" )
print(f"尖度: {round(pd.Series(pred_exp).kurt(),4)}" )
まさにハリウッドスターの物件ような高価な物件の予測結果があって見づらいですね。
予測結果としては必要ですが、一旦学習データ同様に400,000ドル以上の物件は外れ値とみなして描画してみましょう。
# 400,000より高い物件は除去
pred_exp_ex_outliars = pred_exp[pred_exp<400000]
# 指数変換した予測結果をプロット
sns.distplot(pred_exp_ex_outliars)
# 歪度と尖度
print(f"歪度: {round(pd.Series(pred_exp_ex_outliars).skew(),4)}" )
print(f"尖度: {round(pd.Series(pred_exp_ex_outliars).kurt(),4)}" )
# 学習データの住宅価格をプロット(外れ値除去済み)
sns.distplot(train_df['SalePrice'])
# 歪度と尖度
print(f"歪度: {round(pd.Series(train_df['SalePrice']).skew(),4)}" )
print(f"尖度: {round(pd.Series(train_df['SalePrice']).kurt(),4)}" )
#6 提出用csvを作成
sample_submission.csvのフォーマットに従ってデータフレームを作成し、csv形式で出力して提出します。
# sample_submission.csvの読み込み
submission_df = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/sample_submission.csv')
# sample_submission.csvの形式を確認するために先頭五行を見てみる。
submission_df.head()
1列目にID、2列目に予測した住宅価格のサンプルが入ってますね。
基本はSalePrice列に予測結果を代入するだけでOKです。
予測結果の個数が間違っている(例えば外れ値の削除などで行が誤って削除されている)とエラーになるので、フォーマットチェックにもなります。
# 指数変換した値を代入
submission_df['SalePrice'] = pred_exp
エラー出ませんでした!大丈夫そうです!
最後にcsvとして出力しましょう。
ここで、必ず引数index=Falseとしておいてください。(または読み込み時にIDをindexに設定しておいてもOK)
行番号も出力されてしまい、提出フォーマットとずれてしまいます。
# submission.csvを出力
submission_df.to_csv('submission.csv',index=False)
あとは出力されたcsvファイルを提出するだけです。
(提出の方法はいくつかありますが、ここでは割愛します。)
フォーマットにも問題なく、提出完了です!
スコアは0.14125、パラメータチューニング時点でのベストスコアは0.0985だったので過学習気味かもしれませんね。
順位は5112人中2383位(上位47%以内)でした。
予測結果としてはまだまだですが、回帰分析の基本的な流れを学ぶことができました。
この課題でのこれ以上の精度向上は狙わず、今回学んだ手法が他の課題でどれだけ活かせるのか、別の課題に挑戦してみたいと思います。
それでは!