1. 分析
ここからは実際にモデルを構築し、分析と性能評価を行う。
前処理、EDA編はこちら。
今回は中古マンションの取引価格を予測する回帰のタスク。
線形回帰、非ディープラーニング、ディープラーニングの3つの手法を用いて精度を比較する。
それぞれ、重回帰分析、LightGBM、Tabnetを使ってモデルを構築する。
-
重回帰分析
2以上の説明変数を用いて目的変数を予測。 -
LightGBM
テーブルデータに対し決定木をベースに、勾配ブースティングによりアンサンブル学習した手法。テーブルデータの分析コンペ等で使用され高速で精度が高い。
- Tabnet
テーブルデータ向けのDNN(Deep Neural Network)モデル。LightGBM等と同様にテーブルデータの分析コンペ等で使用されている。
先ずは必要なライブラリのインポートから。(今後分析で使用するライブラリを一括して表示)
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.linear_model import SGDRegressor, LinearRegression
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.model_selection import KFold, StratifiedKFold, GridSearchCV
from sklearn import preprocessing
import warnings
import lightgbm as lgb
pip install japanize-matplotlib
import japanize_matplotlib
前処理、EDA編で整形し保存したデータを読み込む。
df = pd.read_csv("mid.csv",encoding ='UTF-8')
1. 線形回帰
特徴量を追加してモデルを構築する。
1. 説明変数の絞り込み
中身を確認し影響がありそうな説明変数に絞る。
質的データについてはヒートマップで相関を可視化しており、内容を見てもそれなりに相関がありそうなので全て使用する。
量的データについて見ていく。
種類と地道府県は全て同じデータの為、削除。
df = df.drop(["種類", "都道府県名"], axis=1)
今後の利用目的と取引の事情等が気になったので内容を見てみる。
先ずは、今後の利用目的から。
df["今後の利用目的"].value_counts()
住宅 41803
その他 2106
事務所 673
店舗 220
Name: 今後の利用目的, dtype: int64
カラム名のとおり、将来的な利用予定目的が表されていた。
未定(欠損)項目もあるが、現在の利用目的を表している用途カラムと強い連関がありそうなので説明変数から除外する。
次に取引の事情等。
df["取引の事情等"].value_counts()
調停・競売等 1608
関係者間取引 155
瑕疵有りの可能性 12
その他事情有り 8
他の権利・負担付き 5
関係者間取引、調停・競売等 1
関係者間取引、瑕疵有りの可能性 1
Name: 取引の事情等, dtype: int64
価格に影響がある追加情報が表示されていた。
内容については、調停・競売等や瑕疵あり物件等であり、極めて特殊性が強く一般的な売買とは言いがたいため説明変数から除外する。
df = df.drop(["今後の利用目的", "取引の事情等"], axis=1)
種類、都道府県名、今後の利用目的、取引の事情等を削除。
2. 特徴量の追加
特徴量の追加ライブラリであるfeaturetoolsを用いて総当り特徴量エンジニアリングを行う。
今回は四則演算系の列同士の各行を計算させるtransformタイプを使用し列同士の足し算、掛け算を行い特徴量を追加する。
pip install featuretools
import featuretools as ft
es = ft.EntitySet(id='example')
es.add_dataframe(dataframe_name='locations',
dataframe=df,
index='name')
trans_primitives=['add_numeric', 'multiply_numeric']
agg_primitives=[]
feature_matrix, feature_defs = ft.dfs(entityset=es,
target_dataframe_name='locations',
trans_primitives=trans_primitives,
agg_primitives=agg_primitives,
max_depth=1)
特徴量の確認
<class 'pandas.core.frame.DataFrame'>
Int64Index: 86538 entries, 0 to 86537
Data columns (total 57 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 市区町村名 86538 non-null category
1 地区名 86535 non-null category
2 最寄駅：名称 86393 non-null category
3 最寄駅：距離（分） 84807 non-null Int64
4 取引価格（総額） 86538 non-null int64
5 間取り 81421 non-null category
6 面積（㎡） 86538 non-null float64
7 建物の構造 85427 non-null category
8 用途 78688 non-null category
9 都市計画 85829 non-null category
10 建ぺい率（％） 85367 non-null Int64
11 容積率（％） 85367 non-null Int64
12 取引時点 86538 non-null float64
13 改装 75347 non-null category
14 築年数 85649 non-null Int64
15 取引価格（総額） + 取引時点 86538 non-null float64
16 取引価格（総額） + 容積率（％） 85367 non-null float64
17 取引価格（総額） + 建ぺい率（％） 85367 non-null float64
18 取引価格（総額） + 最寄駅：距離（分） 84807 non-null float64
19 取引価格（総額） + 築年数 85649 non-null float64
20 取引価格（総額） + 面積（㎡） 86538 non-null float64
21 取引時点 + 容積率（％） 85367 non-null float64
22 取引時点 + 建ぺい率（％） 85367 non-null float64
23 取引時点 + 最寄駅：距離（分） 84807 non-null float64
24 取引時点 + 築年数 85649 non-null float64
25 取引時点 + 面積（㎡） 86538 non-null float64
26 容積率（％） + 建ぺい率（％） 85367 non-null float64
27 容積率（％） + 最寄駅：距離（分） 84058 non-null float64
28 容積率（％） + 築年数 84628 non-null float64
29 容積率（％） + 面積（㎡） 85367 non-null float64
30 建ぺい率（％） + 最寄駅：距離（分） 84058 non-null float64
31 建ぺい率（％） + 築年数 84628 non-null float64
32 建ぺい率（％） + 面積（㎡） 85367 non-null float64
33 最寄駅：距離（分） + 築年数 84015 non-null float64
34 最寄駅：距離（分） + 面積（㎡） 84807 non-null float64
35 築年数 + 面積（㎡） 85649 non-null float64
36 取引価格（総額） * 取引時点 86538 non-null float64
37 取引価格（総額） * 容積率（％） 85367 non-null float64
38 取引価格（総額） * 建ぺい率（％） 85367 non-null float64
39 取引価格（総額） * 最寄駅：距離（分） 84807 non-null float64
40 取引価格（総額） * 築年数 85649 non-null float64
41 取引価格（総額） * 面積（㎡） 86538 non-null float64
42 取引時点 * 容積率（％） 85367 non-null float64
43 取引時点 * 建ぺい率（％） 85367 non-null float64
44 取引時点 * 最寄駅：距離（分） 84807 non-null float64
45 取引時点 * 築年数 85649 non-null float64
46 取引時点 * 面積（㎡） 86538 non-null float64
47 容積率（％） * 建ぺい率（％） 85367 non-null float64
48 容積率（％） * 最寄駅：距離（分） 84058 non-null float64
49 容積率（％） * 築年数 84628 non-null float64
50 容積率（％） * 面積（㎡） 85367 non-null float64
51 建ぺい率（％） * 最寄駅：距離（分） 84058 non-null float64
52 建ぺい率（％） * 築年数 84628 non-null float64
53 建ぺい率（％） * 面積（㎡） 85367 non-null float64
54 最寄駅：距離（分） * 築年数 84015 non-null float64
55 最寄駅：距離（分） * 面積（㎡） 84807 non-null float64
56 築年数 * 面積（㎡） 85649 non-null float64
dtypes: Int64(4), category(8), float64(44), int64(1)
memory usage: 34.2 MB
feature_matrix.shape
(86538, 57)
特徴量が15個から57個まで増加。
但し、目的変数を基に作成された特徴量は相関が強いと判断して削除する。
feature_matrix = feature_matrix.drop(["取引価格（総額） + 取引時点",
"取引価格（総額） + 容積率（％）",
"取引価格（総額） + 建ぺい率（％）",
"取引価格（総額） + 最寄駅：距離（分）",
"取引価格（総額） + 築年数",
"取引価格（総額） + 面積（㎡）",
"取引価格（総額） * 取引時点",
"取引価格（総額） * 容積率（％）",
"取引価格（総額） * 建ぺい率（％）",
"取引価格（総額） * 最寄駅：距離（分）",
"取引価格（総額） * 築年数",
"取引価格（総額） * 面積（㎡）"], axis=1)
3. 欠損値処理
重回帰分析を行うにあたり、欠損値は何らかの値で補完または削除する必要があるが、補完する場合、平均値や最頻値等では値の根拠に妥当性を欠くので、今回は欠損値を削除する。
欠損値確認、削除
feature_matrix.isnull().sum()
市区町村名 0
地区名 3
最寄駅：名称 145
最寄駅：距離（分） 1731
取引価格（総額） 0
間取り 5117
面積（㎡） 0
建物の構造 1111
用途 7850
都市計画 709
建ぺい率（％） 1171
容積率（％） 1171
取引時点 0
改装 11191
築年数 889
取引時点 + 容積率（％） 1171
取引時点 + 建ぺい率（％） 1171
取引時点 + 最寄駅：距離（分） 1731
取引時点 + 築年数 889
取引時点 + 面積（㎡） 0
容積率（％） + 建ぺい率（％） 1171
容積率（％） + 最寄駅：距離（分） 2480
容積率（％） + 築年数 1910
容積率（％） + 面積（㎡） 1171
建ぺい率（％） + 最寄駅：距離（分） 2480
建ぺい率（％） + 築年数 1910
建ぺい率（％） + 面積（㎡） 1171
最寄駅：距離（分） + 築年数 2523
最寄駅：距離（分） + 面積（㎡） 1731
築年数 + 面積（㎡） 889
取引時点 * 容積率（％） 1171
取引時点 * 建ぺい率（％） 1171
取引時点 * 最寄駅：距離（分） 1731
取引時点 * 築年数 889
取引時点 * 面積（㎡） 0
容積率（％） * 建ぺい率（％） 1171
容積率（％） * 最寄駅：距離（分） 2480
容積率（％） * 築年数 1910
容積率（％） * 面積（㎡） 1171
建ぺい率（％） * 最寄駅：距離（分） 2480
建ぺい率（％） * 築年数 1910
建ぺい率（％） * 面積（㎡） 1171
最寄駅：距離（分） * 築年数 2523
最寄駅：距離（分） * 面積（㎡） 1731
築年数 * 面積（㎡） 889
dtype: int64
feature_matrix = feature_matrix.dropna()
feature_matrix.shape
(65370, 45)
特徴量が57個から45個まで減少。
4. モデル構築
前処理としてカテゴリー変数を数値化する。
One hot Encoding等の手法があるが、今回は各カテゴリに対して数値ラベルを割り振るLabel Encodingにより処理を行う。(One hot Encodingも試したがダミー変数化するとカラム数が2000を超え増えすぎてしまいモデル構築がうまくいかなかった為、今回はLabel Encodingを採用。)
def data_preprocess(df):
cat_fetures = ["市区町村名", "地区名", "最寄駅：名称", "間取り", "建物の構造", "用途", "都市計画", "改装"]
for col in cat_fetures:
lbl = preprocessing.LabelEncoder()
lbl.fit(df[col])
lbl.transform(df[col])
df[col] = lbl.transform(df[col])
return df
feature_matrix =data_preprocess(feature_matrix)
説明変数と目的変数をそれぞれX、yと定義、hold-out法により学習データとテストデータにXとyを分割。結果が同じになるように乱数シードは固定。
X = feature_matrix.drop("取引価格（総額）", axis=1)
y = feature_matrix["取引価格（総額）"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
正規方程式を用いたLinearRegression、勾配降下法を用いたSGDRegressorの2種類の手法があるが、今回はパラメータ調整もいらず速度も早いLinearRegressionを用いてモデルを構築、学習、予測を行う。
model = LinearRegression()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
5. 評価・可視化
予測値の主要な各評価指標と、汎化性を算出。
print("平均二乗誤差 (MSE): ", mean_squared_error(y_test, y_pred))
print("二乗平均平方根誤差(RMSE): ", mean_squared_error(y_test, y_pred, squared=False))
print("平均絶対誤差 (MAE): ", mean_absolute_error(y_test, y_pred))
print("決定係数(R2): ", r2_score(y_test, y_pred))
print("平均絶対誤差率(MAPE): ", np.mean(np.abs((y_pred - y_test) / y_test)) * 100)
平均二乗誤差 (MSE): 30523098029306.727
二乗平均平方根誤差(RMSE): 5524771.310136441
平均絶対誤差 (MAE): 4217789.222480547
決定係数(R2): 0.6646073829056789
平均絶対誤差率(MAPE): 172.27180488953232
print("学習データに対する決定係数(R2): ", model.score(X_train ,y_train))
print("学習データに対する決定係数(R2): ", model.score(X_test, y_test))
学習データに対する決定係数(R2): 0.6498311050316787
学習データに対する決定係数(R2): 0.6646073829056789
MSE等は値が大きすぎてどう誤差があるのかが分かりにくいため、直感的に分かりやすい回帰モデルの当てはまりの良さを表す決定係数（R2）を基準に比較していく。
可視化して予測と実測を比較。
plot_data = pd.DataFrame()
plot_data['observed'] = y_test[:]
plot_data['predict'] = y_pred[:]
plt.figure(figsize=(10,10))
sns.jointplot(x='observed', y='predict', data=plot_data, kind='reg')
plt.show()
plt.figure(figsize=(5, 5))
#plt.rcParams['agg.path.chunksize'] = 1000000
plt.plot([0,y_test.max()], [0,y_test.max()], label="y=x", color="black")
plt.scatter(y_pred, y_test)
plt.legend()
plt.xlim(0, y_test.max())
plt.ylim(0, y_test.max())
plt.xlabel("predict")
plt.ylabel("true")
plt.show()
plt.hist(y_test, bins=100, alpha=0.5, label="test")
plt.hist(y_pred, bins=100, alpha=0.5, label="predict")
plt.legend()
plt.show()
似たような分布であるが、ヒストグラムを見てもなんとも言えないというのが正直な所である。
| 線形回帰(featuretools) | R2 |
|---|---|
| R2 | 0.6646 |
参考までにfeaturetoolsを使用せず線形回帰を行うと、R2は0.639であり特徴量エンジニアリングを行った方が若干精度は向上した。
2. lightGBM
1. 特徴量追加・チューニングの無し場合
1. 前処理
データの読み込み、カラムの削除、カテゴリーデータのラベルエンコーディング。(コードは同じなので省略)
2. 欠損値処理
lightGBMでは欠損値があっても学習が進むので、削除等の処理は行わない。
3. モデル構築
線形回帰同様に説明変数と目的変数をそれぞれX、yと定義、hold-out法により学習データとテストデータにXとyを分割、モデル構築。(コードは同じなので省略)
model = lgb.LGBMRegressor(
boosting_type='gbdt',
objective='regression',
n_estimators=1000,
)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
4. 評価・可視化
予測値の主要な各評価指標と、汎化性を算出。
print("平均二乗誤差 (MSE): ", mean_squared_error(y_test, y_pred))
print("二乗平均平方根誤差(RMSE): ", mean_squared_error(y_test, y_pred, squared=False))
print("平均絶対誤差 (MAE): ", mean_absolute_error(y_test, y_pred))
print("決定係数(R2): ", r2_score(y_test, y_pred))
print("平均絶対誤差率(MAPE): ", np.mean(np.abs((y_pred - y_test) / y_test)) * 100)
平均二乗誤差 (MSE): 11285666120287.799
二乗平均平方根誤差(RMSE): 3359414.5502286255
平均絶対誤差 (MAE): 2319922.6782098874
決定係数(R2): 0.8775596649344537
平均絶対誤差率(MAPE): 129.0757449830768
print("学習データに対する決定係数(R2): ", model.score(X_train ,y_train))
print("学習データに対する決定係数(R2): ", model.score(X_test, y_test))
学習データに対する決定係数(R2): 0.9227912206040824
テストデータに対する決定係数(R2): 0.8775596649344537
やや過学習気味であるが、線形回帰と比べて高い精度となった。
| 線形回帰(featuretools) | lightGBM | |
|---|---|---|
| R2 | 0.6646 | 0.8775 |
各特徴量の重要度を可視化。
lgb.plot_importance(model, importance_type='split', figsize=(10,10))
plt.show()
| 重要特徴量 | Hold-out法 |
|---|---|
| 1位 | 築年数 |
| 2位 | 取引時点 |
| 3位 | 最寄駅：名称 |
| 4位 | 地区名 |
| 5位 | 面積（m2） |
面積がもう少し上位に位置しているかと思ったが、概ね思っていた結果であった。
線形回帰同様にグラフで可視化。
plt.figure(figsize=(5, 5))
plt.plot([0,y_test.max()], [0,y_test.max()], label="y=x", color="black")
plt.scatter(y_pred, y_test)
plt.legend()
plt.xlim(0, y_test.max())
plt.ylim(0, y_test.max())
plt.xlabel("predict")
plt.ylabel("true")
plt.show()
plt.hist(y_test, bins=100, alpha=0.5, label="test")
plt.hist(y_pred, bins=100, alpha=0.5, label="predict")
plt.legend()
plt.show()
2. 特徴量を追加をした場合
特徴量を追加してモデルを構築する。今回は特徴量を自動で作成するライブラリを2つ用いて内容を比較してみる。
1.AutoFeat
先ずはAutoFeatを用いて特徴量を追加する。
前処理、ラベルエンコーディングはこれまでと同様なので割愛、ライブラリをインポート。
pip install autofeat
from autofeat import AutoFeatRegressor
autofeatを使用するため欠損値は削除。
df = df.dropna()
説明変数と目的変数の定義、hold-out法により学習データとテストデータにXとyを分割。（コードはこれまでと同様なので省略）
autofeatにより訓練データに対し特徴量を追加。
af = AutoFeatRegressor(verbose=1)
X_train_feature_creation = af.fit_transform(X_train, y_train)
print("number of features in X_train:",
X_train.shape[1])
print("number of features in X_train_feature_creation:",
X_train_feature_creation.shape[1])
number of features in X_train: 14
number of features in X_train_feature_creation: 94
80個の特徴量が生成、テストデータに対しても同様の処理を行う。
X_test_feature_creation = af.transform(X_test)
[AutoFeat] Computing 80 new features.
[AutoFeat] 80/ 80 new features ...done.
特徴量追加後のデータを用いてモデル構築、評価を行う。
model = lgb.LGBMRegressor(
boosting_type='gbdt',
objective='regression',
n_estimators=1000,
)
model.fit(X_train_feature_creation, y_train)
y_pred = model.predict(X_test_feature_creation)
予測値の主要な各評価指標と、汎化性を算出。
print("平均二乗誤差 (MSE): ", mean_squared_error(y_test, y_pred))
print("二乗平均平方根誤差(RMSE): ", mean_squared_error(y_test, y_pred, squared=False))
print("平均絶対誤差 (MAE): ", mean_absolute_error(y_test, y_pred))
print("決定係数(R2): ", r2_score(y_test, y_pred))
print("平均絶対誤差率(MAPE): ", np.mean(np.abs((y_pred - y_test) / y_test)) * 100)
平均二乗誤差 (MSE): 11273346609519.227
二乗平均平方根誤差(RMSE): 3357580.4695523274
平均絶対誤差 (MAE): 2314139.5228874027
決定係数(R2): 0.8764455053951168
平均絶対誤差率(MAPE): 60.42539249346315
print('学習データに対する決定係数(R2): {}'.format(model.score(X_train_feature_creation, y_train)))
print('テストデータに対する決定係数(R2): {}'.format(model.score(X_test_feature_creation, y_test)))
学習データに対する決定係数(R2): 0.9422574227255046
テストデータに対する決定係数(R2): 0.8764455053951168
| lightGBM | Hold-out法 | autofeat |
|---|---|---|
| R2 | 0.8775 | 0.8764 |
過学習気味であり精度が若干低下した。あまり効果的な特徴量が追加されなかったのであろうか。
特徴量の重要度を可視化。
import japanize_matplotlib
lgb.plot_importance(model, importance_type='split', figsize=(10,30))
plt.show()
| 重要特徴量 | Hold-out法 | autofeat |
|---|---|---|
| 1位 | 築年数 | 取引時点 |
| 2位 | 取引時点 | 面積/築年数 |
| 3位 | 最寄駅：名称 | 地区名**3 * 最寄駅距離分 **2 |
| 4位 | 地区名 | sqrt(最寄駅距離分)/面積 |
| 5位 | 面積（m2） | log(築年数)/面積 |
追加した特徴量が効いているようであるが、ラベルエンコーディング後に特徴量が作成される為、カテゴリーデータを基に作成された特徴量も存在し、一概に似たような特徴量重要度であるとは言えないようである。（例えば3位）
単独特徴量で見みると市町村名と最寄駅：名称が比較的上位に存在した。
2. featuretools
次にfeaturetoolsを用いて特徴量を追加する、前処理等は同様なので省略。
線形回帰の場合と同様に、add_numericで特徴量間の和を、multiply_numericで特徴量間の積を追加。
pip install featuretools
import featuretools as ft
es = ft.EntitySet(id='example')
es.add_dataframe(dataframe_name='locations',
dataframe=df,
index='name')
trans_primitives=['add_numeric', 'multiply_numeric']
agg_primitives=[]
feature_matrix, feature_defs = ft.dfs(entityset=es,
target_dataframe_name='locations',
trans_primitives=trans_primitives,
agg_primitives=agg_primitives,
max_depth=1)
特徴量が57個まで増加。
特徴量の確認
<class 'pandas.core.frame.DataFrame'>
Int64Index: 86538 entries, 0 to 86537
Data columns (total 57 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 市区町村名 86538 non-null category
1 地区名 86535 non-null category
2 最寄駅：名称 86393 non-null category
3 最寄駅：距離（分） 84807 non-null Int64
4 取引価格（総額） 86538 non-null int64
5 間取り 81421 non-null category
6 面積（㎡） 86538 non-null float64
7 建物の構造 85427 non-null category
8 用途 78688 non-null category
9 都市計画 85829 non-null category
10 建ぺい率（％） 85367 non-null Int64
11 容積率（％） 85367 non-null Int64
12 取引時点 86538 non-null float64
13 改装 75347 non-null category
14 築年数 85649 non-null Int64
15 取引価格（総額） + 取引時点 86538 non-null float64
16 取引価格（総額） + 容積率（％） 85367 non-null float64
17 取引価格（総額） + 建ぺい率（％） 85367 non-null float64
18 取引価格（総額） + 最寄駅：距離（分） 84807 non-null float64
19 取引価格（総額） + 築年数 85649 non-null float64
20 取引価格（総額） + 面積（㎡） 86538 non-null float64
21 取引時点 + 容積率（％） 85367 non-null float64
22 取引時点 + 建ぺい率（％） 85367 non-null float64
23 取引時点 + 最寄駅：距離（分） 84807 non-null float64
24 取引時点 + 築年数 85649 non-null float64
25 取引時点 + 面積（㎡） 86538 non-null float64
26 容積率（％） + 建ぺい率（％） 85367 non-null float64
27 容積率（％） + 最寄駅：距離（分） 84058 non-null float64
28 容積率（％） + 築年数 84628 non-null float64
29 容積率（％） + 面積（㎡） 85367 non-null float64
30 建ぺい率（％） + 最寄駅：距離（分） 84058 non-null float64
31 建ぺい率（％） + 築年数 84628 non-null float64
32 建ぺい率（％） + 面積（㎡） 85367 non-null float64
33 最寄駅：距離（分） + 築年数 84015 non-null float64
34 最寄駅：距離（分） + 面積（㎡） 84807 non-null float64
35 築年数 + 面積（㎡） 85649 non-null float64
36 取引価格（総額） * 取引時点 86538 non-null float64
37 取引価格（総額） * 容積率（％） 85367 non-null float64
38 取引価格（総額） * 建ぺい率（％） 85367 non-null float64
39 取引価格（総額） * 最寄駅：距離（分） 84807 non-null float64
40 取引価格（総額） * 築年数 85649 non-null float64
41 取引価格（総額） * 面積（㎡） 86538 non-null float64
42 取引時点 * 容積率（％） 85367 non-null float64
43 取引時点 * 建ぺい率（％） 85367 non-null float64
44 取引時点 * 最寄駅：距離（分） 84807 non-null float64
45 取引時点 * 築年数 85649 non-null float64
46 取引時点 * 面積（㎡） 86538 non-null float64
47 容積率（％） * 建ぺい率（％） 85367 non-null float64
48 容積率（％） * 最寄駅：距離（分） 84058 non-null float64
49 容積率（％） * 築年数 84628 non-null float64
50 容積率（％） * 面積（㎡） 85367 non-null float64
51 建ぺい率（％） * 最寄駅：距離（分） 84058 non-null float64
52 建ぺい率（％） * 築年数 84628 non-null float64
53 建ぺい率（％） * 面積（㎡） 85367 non-null float64
54 最寄駅：距離（分） * 築年数 84015 non-null float64
55 最寄駅：距離（分） * 面積（㎡） 84807 non-null float64
56 築年数 * 面積（㎡） 85649 non-null float64
dtypes: Int64(4), category(8), float64(44), int64(1)
memory usage: 34.2 MB
feature_matrix.shape
(86538, 57)
線形回帰同様、目的変数を元に作成された特徴量は強い相関を持つと判断し、それらの特徴量は削除する。
特徴量を削除
```py
feature_matrix = feature_matrix.drop(["取引価格（総額） + 取引時点",
"取引価格（総額） + 容積率（％）",
"取引価格（総額） + 建ぺい率（％）",
"取引価格（総額） + 最寄駅：距離（分）",
"取引価格（総額） + 築年数",
"取引価格（総額） + 面積（㎡）",
"取引価格（総額） * 取引時点",
"取引価格（総額） * 容積率（％）",
"取引価格（総額） * 建ぺい率（％）",
"取引価格（総額） * 最寄駅：距離（分）",
"取引価格（総額） * 築年数",
"取引価格（総額） * 面積（㎡）"
], axis=1
```
特徴量追加後のデータを用いてモデル構築、評価。
ラベルエンコーディング、説明変数と目的変数の定義、hold-out法により学習データとテストデータにXとyを分割しモデル構築。（コードはこれまでと同様なので省略）
予測値の主要な各評価指標と、汎化性を算出。
print("平均二乗誤差 (MSE): ", mean_squared_error(y_test, y_pred))
print("二乗平均平方根誤差(RMSE): ", mean_squared_error(y_test, y_pred, squared=False))
print("平均絶対誤差 (MAE): ", mean_absolute_error(y_test, y_pred))
print("決定係数(R2): ", r2_score(y_test, y_pred))
print("平均絶対誤差率(MAPE): ", np.mean(np.abs((y_pred - y_test) / y_test)) * 100)
平均二乗誤差 (MSE): 11582633831816.389
二乗平均平方根誤差(RMSE): 3403326.877015546
平均絶対誤差 (MAE): 2354495.7853116845
決定係数(R2): 0.8743378058332133
平均絶対誤差率(MAPE): 128.7027955459018
print('学習データに対する決定係数(R2): {}'.format(model.score(X_train, y_train)))
print('テストデータに対する決定係数(R2): {}'.format(model.score(X_test, y_test)))
学習データに対する決定係数(R2): 0.9331375546303661
テストデータに対する決定係数(R2): 0.8743378058332133
| lightGBM | hold-out法 | autofeat | featuretools |
|---|---|---|---|
| R2 | 0.8775 | 0.8764 | 0.8743 |
autofeatの際と同様に、過学習気味であり特徴量を追加せずにhold-out法で実装した方が精度が良かった。
特徴量の重要度を可視化。
lgb.plot_importance(model, importance_type='split', figsize=(20,10))
plt.show()
| 重要特徴量 | Hold-out法 | autofeat | featuretools |
|---|---|---|---|
| 1位 | 築年数 | 取引時点 | 市町村名 |
| 2位 | 取引時点 | 面積/築年数 | 最寄駅:名称 |
| 3位 | 最寄駅：名称 | 地区名**3 * 最寄駅距離分 **2 | 地区名 |
| 4位 | 地区名 | sqrt(最寄駅距離分)/面積 | 取引時点_+_面積（m2） |
| 5位 | 面積（m2） | log(築年数)/面積 | 取引時点_+_築年数 |
市町村名が1位であったのは意外であったが、概ね似たような特徴量が上位に存在した。築年数、面積が比較的下位に存在しているのが気になった。
3. 交差検証とハイパーパラメーターのチューニング
lightGBMではどのモデルもR2が0.88付近となり高精度な結果となった。
但し、hold-out法により訓練データとテストデータの分割を行った際に、たまたま高精度となる分割となった可能性を否定できない為、特徴量追加を行わなかったモデルに対して交差検証を行い精度を比較する。
今回は、交差検証を内側のクロスバリデーション(innerCV)と外側のクロスバリデーション(outerCV)に分け、innerCVではハイパーパラメータの調整に注力し、選択されたハイパーパラメータを用いてouterCVでモデルを5分割し精度を評価するダブルクロスバリデーションを実施する。
学習データとテストデータにXとyを分割するまでのコードは同じである為省略。
model = lgb.LGBMRegressor(
objective='regression',
n_estimators=1000,
)
# 内側のクロスバリデーション（GridSearchCV）のハイパーパラメーターの条件
params = {
'num_leaves' : [10, 100],
'learning_rate' : [0.001, 0.01, 0.1],
'max_depth' : [7, 127, 511],
}
skf_val = StratifiedKFold(n_splits=3) #内側の3分割設定
# 内側のクロスバリデーション（GridSearchCV）
print('innerCV START')
gscv = GridSearchCV(model, params, cv=skf_val, verbose=2, scoring='r2')
gscv.fit(X_train, y_train, verbose= 1)
print('innerCV best params',gscv.best_params_)
print('outerCV START')
skf = StratifiedKFold(n_splits=5)#外側の交差検証数
scores = cross_val_score(gscv, X, y, cv=skf)# 外側の交差検証5分割
print(scores)#外側5回分の交差検証の結果
scores_nested_cv = scores.mean()
print(scores_nested_cv)#外側5回分の交差検証の結果の平均
[0.78437895 0.82667577 0.84881886 0.84450076 0.77004138]
0.8148831435457533
y_pred = gscv.best_estimator_.predict(X_test)
print('lightGBM bestparameter', gscv.best_params_)
lightGBM bestparameter {'learning_rate': 0.1, 'max_depth': 7, 'num_leaves': 100}
| lightGBM | hold-out法 | autofeat | featuretools | ダブルクロスバリデーション |
|---|---|---|---|---|
| R2 | 0.8775 | 0.8764 | 0.8743 | 0.8148 |
ダブルクロスバリデーションの結果、R2の平均は0.8148となった。
hold-out法の際には、たまたま精度が高い分割となっていたと考えられる。
autofeat、featuretools共に、R2がダブルクロスバリデーションの結果よりも大きく、特徴量の追加は精度向上に寄与したと言ってもよさそうである。
3. tabnet
テーブルデータ向けのニューラルネットワークモデルを用いる。
tabnetは決定木ベースのモデルの解釈可能性を持ちつつ、 大規模なテーブルデータに対して高精度な学習が可能であり、特徴量の作成を必要としないエンドツーエンドの学習モデルとなっている。
必要なライブラリの読み込み。
!pip install pytorch-tabnet
import torch
from pytorch_tabnet.tab_model import TabNetRegressor
from pytorch_tabnet.pretraining import TabNetPretrainer
from sklearn import preprocessing
import numpy as np
import pandas as pd
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split
これまでと同様、一連の前処理を行う。
df = pd.read_csv("mid.csv",encoding ='UTF-8')
df = df.drop(["種類", "都道府県名"], axis=1)
df = df.drop(["今後の利用目的", "取引の事情等"], axis=1)
df = df.dropna()
def data_preprocess(df):
cat_fetures = ["市区町村名", "地区名", "最寄駅：名称", "間取り", "建物の構造", "用途", "都市計画", "改装"]
for col in cat_fetures:
lbl = preprocessing.LabelEncoder()
lbl.fit(df[col])
lbl.transform(df[col])
df[col] = lbl.transform(df[col])
return df
df =data_preprocess(df)
X = df.drop(["取引価格（総額）"], axis=1).values
y = df["取引価格（総額）"].values.reshape(-1,1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
X_train, X_valid, y_train, y_valid = train_test_split(X_train, y_train, test_size=.3, random_state=0)
事前学習なしでモデルを作成し、性能を見てみる。
tabnet_params = dict(
seed = 0,
optimizer_params = dict(lr = 1e-2),
n_steps=3,
verbose = 10,
)
model = TabNetRegressor(**tabnet_params)
model.fit(
X_train, y_train,
eval_set = [(X_valid, y_valid)],
batch_size = 32,
max_epochs = 100,
patience = 100,
)
性能評価
preds = model.predict(X_test)
print("RMSE: ",np.sqrt(mean_squared_error(y_test, preds)))
print("R2: ",r2_score(y_test, preds))
RMSE: 4549691.8667240115
R2: 0.772548473442036
続いて自己事前学習ありの場合でモデルを作成し、性能がどう変化するか見てみる。
先ずは事前学習を実施。
pretrainer = TabNetPretrainer(**tabnet_params)
pretrainer.fit(
X_train,
eval_set=[X_valid],
max_epochs=5000,
patience=100,
)
事前学習した結果を from_unsupervised=pretrainer として利用して学習する。
model = TabNetRegressor(**tabnet_params)
model.fit(
X_train, y_train,
eval_set=[(X_valid, y_valid)],
max_epochs=5000,
patience=100,
eval_metric={'rmse'},
from_unsupervised=pretrainer
)
性能評価
preds = model.predict(X_test)
print("RMSE: ",np.sqrt(mean_squared_error(y_test, preds)))
print("R2: ",r2_score(y_test, preds))
RMSE: 4111006.248914816
R2: 0.8142960417553144
| tabnet | 事前学習なし | 事前学習あり |
|---|---|---|
| R2 | 0.7725 | 0.8142 |
事前学習を行った事で精度は向上した。
今回はたまたま精度が向上したのかもしれないが、パラメーターチューニング等により高精度を算出できそうなモデルであった。
特徴量の重要度を可視化。
import matplotlib.pyplot as plt
plt.barh(list(df.columns[:-1])[::-1], model.feature_importances_[::-1])
plt.grid()
plt.show()
| 重要特徴量 | Hold-out法 | autofeat | featuretools | tabnet |
|---|---|---|---|---|
| 1位 | 築年数 | 取引時点 | 市町村名 | 都市計画 |
| 2位 | 取引時点 | 面積/築年数 | 最寄駅:名称 | 面積（m2） |
| 3位 | 最寄駅：名称 | 地区名**3 * 最寄駅距離分 **2 | 地区名 | 用途 |
| 4位 | 地区名 | sqrt(最寄駅距離)/面積 | 取引時点_+_面積（m2） | 最寄駅：名称 |
| 5位 | 面積（m2） | log(築年数)/面積 | 取引時点_+_築年数 | 改装 |
lightGMBとは違った毛色の重要特徴量順位であった。ディープラーニングモデルであるが故であろうか不明であるが、今回は比較だけに留めておく。
2. 終わりに
1. 各モデルの性能評価比較
構築した3つのモデルで最も良かった決定係数（R2）を比較。
| 線形回帰 | lightGBM（autofeat） | tabnet | |
|---|---|---|---|
| R2 | 0.6646 | 0.8764 | 0.8142 |
autofeatを使用し特徴量を追加したlightGBMが一番高い精度となった。
予想していた通り線形回帰が一番低い精度であったが、一概にディープラーニングを用いた方が良い結果に繋がるとは言えないようである。
lightGBMの結果を基に、実測値と予測値の差を見てみる。
pd.options.display.float_format = '{:.2f}'.format
ctual_pred_df_test = pd.DataFrame({"actual": (y_test), "pred": (y_pred)})
actual_pred_df_test["差額"] = actual_pred_df_test["pred"] - actual_pred_df_test["actual"]
actual_pred_df_test
数字だけでは分かりにくいので差額をヒストグラムで表示。
plt.hist(actual_pred_df_test["差額"], bins=100, label="差額")
plt.legend()
plt.show()
外れ値は存在するものの差額の分布は0付近が一番多く、予想通り正規分布となりそこまでおかしな結果ではないようであった。
今回は、フレームワークを用いて特徴量エンジニアリングを実施し精度は向上したが、大きな変化は見られなかった。精度向上にあたっては、データに対するドメイン知識を元に特徴量エンジニアリングを行う重要性を実感した。
使用データは、国土交通省から公開されているものを用いたが、不動産会社から公開されている賃貸物件の情報をwebスクレイピングにより取得し、賃貸価格の予測を行ってもイメージが湧きやすく面白い分析ができそうであると感じた。
今後も知見を広げ、技術向上に努めていきたい。
2. 所感
今回の分析については、データ選定、前処理、EDA、モデル作成、性能評価と機械学習の一連の流れに注力して取り組んだ。
全工程で何をどこまで実施するか全て自分で定義して取り組んだ為、自由研究的な要素が強く、振り返ると必要以上に時間を費やし思うように前に進まない場面があった事も否定できなかった。
実務に取り組む際は、限られた労力や時間でどれだけの成果を生み出せるかという生産性が問われる為、ある意味良い経験となった。
また、pythonを使用したデータ分析の知識だけではなく、並行してHTML、CSS基礎・SQL・GITを使ったバージョン管理・クリティカルシンキング・マーケティングの知識についても必要性を感じ、学習に取り組んだ。
スクールでの学びを活かすと共に、今後も幅広く自分の知見を広げ、社会へ貢献できる人材となれるよう努めていきたい。
3. コード詳細・ディレクトリ構成・HP制作
コードはGitHubを参照。
ディレクトリ構成は以下のとおり。
※ダウンロードデータについては容量が大きい為、整形データのみを掲載。
├── README.md
├── data
│ └── interim
│ └── mid.csv
│
└── notebooks
├── EDA.ipynb
├── linear.ipynb
├── lightGBM_simple.ipynb
├── lightGBM_autofeat.ipynb
├── lightGBM_featuretools.ipynb
├── lightGBM_gridsearch.ipynb
└── tabnet.ipynb
本分析と並行して取り組んだ簡易なHPの作成。