[Python] その2 ボストンの住宅価格を重回帰分析で予測してみた
下記記事の続編になります。
前回の記事では、ボストンの住宅価格 (Boston housing prices)のデータセットの説明を書いたので、データの理解をしたい方はこちらへどうぞ
前回は単回帰での住宅価格の予測ですが、今回は重回帰で予測分析をしました。
しばらくはボストンの住宅価格の予測モデル向上を目指して分析していくことにします。
Checking Boston housing prices corrected dataset 2
パッケージのインポート
import seaborn as sns import statsmodels.api as sm import pandas as pd
データの読み込み
df = pd.read_csv("http://lib.stat.cmu.edu/datasets/boston_corrected.txt",skiprows=9,sep="\t")
※ 22年2月1日更新: UnicodeDecodeErrorが発生する方はencoding='Windows-1252'を追加してください。 詳細は下記新ブログをご参照ください。
pd.DataFrame(df.columns)
Out[0]
| 0 | |
|---|---|
| 0 | OBS. |
| 1 | TOWN |
| 2 | TOWN# |
| 3 | TRACT |
| 4 | LON |
| 5 | LAT |
| 6 | MEDV |
| 7 | CMEDV |
| 8 | CRIM |
| 9 | ZN |
| 10 | INDUS |
| 11 | CHAS |
| 12 | NOX |
| 13 | RM |
| 14 | AGE |
| 15 | DIS |
| 16 | RAD |
| 17 | TAX |
| 18 | PTRATIO |
| 19 | B |
| 20 | LSTAT |
df.head(3)
Out[0]
| OBS. | TOWN | TOWN# | TRACT | LON | LAT | MEDV | CMEDV | CRIM | ZN | INDUS | CHAS | NOX | RM | AGE | DIS | RAD | TAX | PTRATIO | B | LSTAT | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | Nahant | 0 | 2011 | -70.955 | 42.2550 | 24.0 | 24.0 | 0.00632 | 18.0 | 2.31 | 0 | 0.538 | 6.575 | 65.2 | 4.0900 | 1 | 296 | 15.3 | 396.90 | 4.98 |
| 1 | 2 | Swampscott | 1 | 2021 | -70.950 | 42.2875 | 21.6 | 21.6 | 0.02731 | 0.0 | 7.07 | 0 | 0.469 | 6.421 | 78.9 | 4.9671 | 2 | 242 | 17.8 | 396.90 | 9.14 |
| 2 | 3 | Swampscott | 1 | 2022 | -70.936 | 42.2830 | 34.7 | 34.7 | 0.02729 | 0.0 | 7.07 | 0 | 0.469 | 7.185 | 61.1 | 4.9671 | 2 | 242 | 17.8 | 392.83 | 4.03 |
モデリング
前回は単回帰を試したので、今回は重回帰分析を試す
サンプリング (訓練データ、テストデータ分割)
# 訓練データとテストデータに分割する。 # 本当は町ごとにサンプリングした方がいいと思うが、TODOにしておく。 from sklearn.model_selection import train_test_split # TODO:層別サンプリング train, test = train_test_split(df, test_size=0.20, stratify=df["町区分"], random_state=100) train, test = train_test_split(df, test_size=0.20,random_state=100)
# 訓練データの件数確認 train.count()["OBS."]
Out[0]
404
# テストデータの件数確認 test.count()["OBS."]
Out[0]
102
訓練データ、テストデータ作成
X_train = train[["RM","LSTAT"]] #部屋数と身分が低い人口割合での重回帰 Y_train = train["CMEDV"] #修正済み住宅価格中央値 X_test = test[["RM","LSTAT"]] Y_test = test["CMEDV"]
モデル作成
import numpy as np from sklearn.linear_model import LinearRegression model = LinearRegression() model.fit(X_train,Y_train)
Out[0]
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)
# 係数逆転現象の確認 pd.DataFrame({"name":X_train.columns,"coefficients":model.coef_})
Out[0]
| name | coefficients | |
|---|---|---|
| 0 | RM | 4.842410 |
| 1 | LSTAT | -0.618469 |
住宅価格 = 4.8*部屋数 + -0.61貧困人口割合 + 切片 という結果になった。部屋数はプラスに働き、貧困人口割合はマイナスに働いていることがわかる。人の感覚とも前回のデータ確認結果と比べてもおかしくないので問題なさそう。
精度確認
関数定義
# 自由度調整済みr2を算出 def adjusted_r2(X,Y,model): from sklearn.metrics import r2_score import numpy as np r_squared = r2_score(Y, model.predict(X)) adjusted_r2 = 1 - (1-r_squared)*(len(Y)-1)/(len(Y)-X.shape[1]-1) #yhat = model.predict(X) \ #SS_Residual = sum((Y-yhat)**2) \ #SS_Total = sum((Y-np.mean(Y))**2) #r_squared = 1 - (float(SS_Residual))/ SS_Total return adjusted_r2 # 予測モデルの精度確認の各種指標を算出 def get_model_evaluations(X_train,Y_train,X_test,Y_test,model): from sklearn.metrics import explained_variance_score from sklearn.metrics import mean_absolute_error from sklearn.metrics import mean_squared_error from sklearn.metrics import mean_squared_log_error from sklearn.metrics import median_absolute_error # 評価指標確認 # 参考: https://funatsu-lab.github.io/open-course-ware/basic-theory/accuracy-index/ yhat_test = model.predict(X_test) return "adjusted_r2(train) :" + str(adjusted_r2(X_train,Y_train,model)) \ , "adjusted_r2(test) :" + str(adjusted_r2(X_test,Y_test,model)) \ , "平均誤差率(test) :" + str(np.mean(abs(Y_test / yhat_test - 1))) \ , "MAE(test) :" + str(mean_absolute_error(Y_test, yhat_test)) \ , "MedianAE(test) :" + str(median_absolute_error(Y_test, yhat_test)) \ , "RMSE(test) :" + str(np.sqrt(mean_squared_error(Y_test, yhat_test))) \ , "RMSE(test) / MAE(test) :" + str(np.sqrt(mean_squared_error(Y_test, yhat_test)) / mean_absolute_error(Y_test, yhat_test)) #better if result = 1.253
get_model_evaluations(X_train,Y_train,X_test,Y_test,model)
Out[0]
('adjusted_r2(train) :0.6306822565643309',
'adjusted_r2(test) :0.6621931601574613',
'平均誤差率(test) :0.20928344610626637',
'MAE(test) :4.057782817379187',
'MedianAE(test) :3.180979173809634',
'RMSE(test) :5.652366792313536',
'RMSE(test) / MAE(test) :1.3929692757593783')
MAEで確認すると4,058ドル誤差がある。 貧困人口割合のみの単回帰分析よりも精度は改善した。どこまで改善するか挑戦
# 描画設定 from matplotlib import rcParams rcParams['xtick.labelsize'] = 12 # x軸のラベルのフォントサイズ rcParams['ytick.labelsize'] = 12 # y軸のラベルのフォントサイズ rcParams['figure.figsize'] = 18,8 # 画像サイズの変更(inch) import matplotlib.pyplot as plt from matplotlib import ticker sns.set_style("whitegrid") # seabornのスタイルセットの一つ sns.set_color_codes() # デフォルトカラー設定 (deepになってる) plt.figure() ax = sns.regplot(x=Y_test, y=model.predict(X_test), fit_reg=False,color='#4F81BD') ax.set_xlabel(u"CMEDV") ax.set_ylabel(u"(Predicted) CMEDV") ax.get_xaxis().set_major_formatter(ticker.FuncFormatter(lambda x, p: format(int(x), ','))) ax.get_yaxis().set_major_formatter(ticker.FuncFormatter(lambda y, p: format(int(y), ','))) ax.plot([0,10,20,30,40,50],[0,10,20,30,40,50], linewidth=2, color="#C0504D",ls="--")
Out[0]
ほぼぴったり予測できている物件もある。値段が高い住宅は予測値から大きく外れているので他の要因がありそう
属性を増やして再モデリングを実施
X_train = train[["RM","LSTAT","PTRATIO","B"]] #部屋数、身分が低い人口割合、生徒と先生比率、人種スコアでの重回帰 Y_train = train["CMEDV"] #修正済み住宅価格中央値 X_test = test[["RM","LSTAT","PTRATIO","B"]] Y_test = test["CMEDV"]
import numpy as np from sklearn.linear_model import LinearRegression model = LinearRegression() model.fit(X_train,Y_train)
Out[0]
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)
# 係数逆転現象の確認 pd.DataFrame({"name":X_train.columns,"coefficients":model.coef_})
Out[0]
| name | coefficients | |
|---|---|---|
| 0 | RM | 4.565370 |
| 1 | LSTAT | -0.494715 |
| 2 | PTRATIO | -0.914036 |
| 3 | B | 0.010658 |
精度確認
get_model_evaluations(X_train,Y_train,X_test,Y_test,model)
Out[0]
('adjusted_r2(train) :0.6790168470285596',
'adjusted_r2(test) :0.7018704528340367',
'平均誤差率(test) :0.18138429186748675',
'MAE(test) :3.541711639134873',
'MedianAE(test) :2.4934473648984934',
'RMSE(test) :5.256139667400787',
'RMSE(test) / MAE(test) :1.4840676494726415')
MAEで確認すると3,541ドル誤差がある。 平均誤差率も低くなったので精度は向上している。PTRATIOが大きくなるとマイナスに働いているのは先生よりも生徒の比重がより多い理由は例えば高級住宅街ほど子供も人数が多くなく結果として先生一人あたりの生徒の数が少ないと考えられる。Bのスコアがプラスに出ているのは、白人が多いエリアほど大きくプラスになるようなスコアになっており、(当時は)白人が多いエリアは住宅価格が大きい物件だと思われるため。
モデリング結果描画
# 描画設定 from matplotlib import rcParams rcParams['xtick.labelsize'] = 12 # x軸のラベルのフォントサイズ rcParams['ytick.labelsize'] = 12 # y軸のラベルのフォントサイズ rcParams['figure.figsize'] = 18,8 # 画像サイズの変更(inch) import matplotlib.pyplot as plt from matplotlib import ticker sns.set_style("whitegrid") # seabornのスタイルセットの一つ sns.set_color_codes() # デフォルトカラー設定 (deepになってる) plt.figure() ax = sns.regplot(x=Y_test, y=model.predict(X_test), fit_reg=False,color='#4F81BD') ax.set_xlabel(u"CMEDV") ax.set_ylabel(u"(Predicted) CMEDV") ax.get_xaxis().set_major_formatter(ticker.FuncFormatter(lambda x, p: format(int(x), ','))) ax.get_yaxis().set_major_formatter(ticker.FuncFormatter(lambda y, p: format(int(y), ','))) ax.plot([0,10,20,30,40,50],[0,10,20,30,40,50], linewidth=2, color="#C0504D",ls="--")
Out[0]
描画するとより赤線に集約したように見える。しかしまだまだ特徴を捉えきれていない物件が多い
