R|階層線形モデルで渋谷区の賃貸価格を予想する
以前の重回帰による賃貸価格の予測では、『最寄り駅ごとに賃貸価格の母集団分布(例えば各説明変数の母回帰係数)が異なりそう』という仮説を立て、1駅ずつ重回帰モデルを推定していきました。
ただ、駅が変わるとモデルの作り直しをしたりですんごいめんどくさい。
今回はそうした「最寄り駅」ごとに階層構造になっているデータに対する分析手法として知られている「階層線形モデル」を試してそうしたお悩みの解消を試みたいと思います( ˆoˆ )
■ 今回やりたいこと
■ 階層線形モデルとは
切片や傾きが各グループ(今回であれば各駅)ごとに正規分布に従って得られると仮定してモデルを推定します。
階層線形モデルの嬉しいところ
■ lme4::lmer使い方
階層線形モデル作成のために lme4 パッケージ を使用します。
1. 基本
> library(lme4)
> model1 <- lmer(y ~ x1 + x2 + x3 + (1 + x1 | factor), data=data)
> model2 <- lmer(y ~ x1 + x2 + x3 + (0 + x1 | factor), data=data)
2. ランダム効果変数が複数の場合
> model3 <- lmer(y ~ x1 + x2 + x3 + (1 + x2 + x3 | factor), data=data)
3. 回帰係数間の相関関係がないと仮定する場合
階層線形モデルでは、説明変数の回帰係数の間に相関関係を考えることが出来ます。例えば、敷地面積の回帰係数が大きければ駅徒歩の回帰係数の大きい傾向にあるなど。このような関係がないと仮定したいときは、以下のようにformulaを書きます。
> model4 <- lmer(y ~ x1 + x2 + x3 + (1 + x1 | factor), data=data)
> model5 <- lmer(y ~ x1 + x2 + x3 + (1 + x1 || factor), data=data)
■ データ観察
1. 駅の数と駅ごとのサンプルサイズを調べる
> rooms %>% group_by(station) %>% summarise(n())
# A tibble: 18 x 2
station `n()`
<chr> <int>
1 代々木上原駅 405
2 代々木八幡駅 369
3 代々木公園駅 353
4 代々木駅 304
5 代官山駅 219
6 初台駅 608
7 北参道駅 168
8 千駄ヶ谷駅 50
9 南新宿駅 56
10 原宿駅 96
11 参宮橋駅 248
12 国立競技場駅 4
13 幡ヶ谷駅 730
14 恵比寿駅 1118
15 明治神宮前駅 37
16 渋谷駅 651
17 神泉駅 156
18 笹塚駅 787
国立競技場駅はサンプルサイズが4件とあまりにも小さく、適切なモデルの推定が出来ないと考えたため、今回は削除することにしました。
2. 駅ごとに面積と家賃の散布図を描いて、外れ値そうなのがないかチェック
> ggplot(rooms, aes(x=space, y=price)) + facet_wrap( ~ station) + geom_point() + stat_smooth(method = "lm", se = FALSE) + theme_gray(base_family="HiraKakuProN-W3")
3. 極端な外れ値を除去
4. 対数を取らない場合と対数を取った場合の散布図を確認
■ 階層線形モデルでモデル作成
実際の賃貸データを使ってモデルを作っていきます。
検証用にまずはデータを訓練データとテストデータに分割します。
1. データの分割
> library(tidymodels)
> splited_data <- initial_split(rooms, p = 0.75, strata = c('station'))
> training_data <- training(splited_data)
> test_data <- testing(splited_data)
各駅ごとにテストデータが1/4の割合になるよう、層別サンプリング(stratified shuffle split)を行う。
2. lmer関数で階層線形モデルを推定
> model <- lmer(log_price ~ log_space + (1 + log_space || station), data=training_data)
3. summary結果の確認
> summary(model)
Linear mixed model fit by REML ['lmerMod']
Formula: log_price ~ log_space + ((1 | station) + (0 + log_space | station))
Data: training_data
REML criterion at convergence: -3040.5
Scaled residuals:
Min 1Q Median 3Q Max
-3.8413 -0.6406 0.0344 0.6507 4.6524
Random effects:
Groups Name Variance Std.Dev.
station (Intercept) 0.059952 0.24485
station.1 log_space 0.008947 0.09459
Residual 0.029779 0.17257
Number of obs: 4725, groups: station, 17
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.71476 0.06575 -10.87
log_space 0.95903 0.02425 39.55
Correlation of Fixed Effects:
(Intr)
log_space -0.137
専有面積の固定効果(Fixed effects)の回帰係数が正の値になっているため、面積が広いほど価格が高そうというドメイン知識と合っていそうです。
↑で得た恵比寿の切片とlog_spaceの係数で回帰直線を確認します。
> plot(data_without_kokuritsu[data_without_kokuritsu$station == '恵比寿駅', c('log_space', 'log_price')])
> abline(-0.711884, 0.9951441)
■ モデル評価
1. 決定係数の確認
各駅ごとに決定係数を確認したいと思います。
# A tibble: 18 x 8
station training test log_training log_test training_size test_size all_size
<chr> <dbl> <dbl> <dbl> <dbl> <int> <int> <int>
1 渋谷駅 0.825 0.842 0.883 0.902 484 161 645
2 千駄ヶ谷駅 0.823 0.775 0.914 0.855 34 15 49
3 代々木上原駅 0.925 0.915 0.932 0.910 299 103 402
4 原宿駅 0.784 0.956 0.845 0.902 67 29 96
5 明治神宮前駅 0.868 0.826 0.924 0.778 24 13 37
6 代官山駅 0.799 0.839 0.835 0.880 159 51 210
7 代々木駅 0.815 0.817 0.774 0.785 238 65 303
8 初台駅 0.834 0.802 0.842 0.845 462 139 601
9 笹塚駅 0.769 0.753 0.775 0.723 581 200 781
10 代々木公園駅 0.855 0.843 0.874 0.834 252 95 347
11 神泉駅 0.909 0.883 0.900 0.871 121 35 156
12 南新宿駅 0.918 0.790 0.849 0.771 44 12 56
13 代々木八幡駅 0.819 0.840 0.871 0.865 278 85 363
14 北参道駅 0.774 0.774 0.836 0.829 127 40 167
15 幡ヶ谷駅 0.807 0.838 0.796 0.833 538 186 724
16 恵比寿駅 0.844 0.813 0.870 0.860 836 278 1114
17 参宮橋駅 0.914 0.876 0.904 0.883 181 67 248
18 渋谷区全体 0.871 0.868 0.894 0.892 4725 1574 6303
2. 残差プロットの確認
上下のばらつきも偏りがなく、0 を中心としてまとまっているのでよさそう。
3. Q-Qプロットの確認
HLMdiagパッケージで物件レベルのQ-Qプロットを確認します。
> library(HLMdiag)
> resid1 <- HLMresid(model, level = 1, type = 'LS', standardize = TRUE)
> ggplot_qqnorm(x = na.omit(resid1$LS.resid), line = 'rlm')
若干右上が外れていますが概ねよさそうです。
以上から今回のモデルは妥当だと言えそうです。
■ モデルの考察
渋谷区全体
渋谷区全体では、0.96乗を1乗として見たときに、面積が1㎡広くなると賃料は4,900円増加するようです。
※ spaceの指数: 面積の価格に対する弾力性という
各駅ごと
(Intercept) log_space space_elasticity
代々木上原駅 -5.246001e-02 -0.003368677 0.9556636
代々木八幡駅 -9.618074e-02 0.028776319 0.9878086
代々木公園駅 -3.121204e-02 0.035172277 0.9942046
代々木駅 1.579823e-01 -0.042189230 0.9168431
代官山駅 -5.421307e-02 0.051477602 1.0105099
初台駅 1.395273e-01 -0.093421536 0.8656108
北参道駅 -6.300289e-01 0.215660779 1.1746931
千駄ヶ谷駅 -4.795593e-02 0.002003871 0.9610362
南新宿駅 9.810516e-02 -0.046588307 0.9124440
原宿駅 -1.255493e-01 0.065729577 1.0247619
参宮橋駅 2.393348e-01 -0.089726112 0.8693062
幡ヶ谷駅 2.133135e-01 -0.127262193 0.8317701
恵比寿駅 2.875997e-03 0.036114078 0.9951464
明治神宮前駅 1.183978e-05 0.035050543 0.9940828
渋谷駅 -1.756146e-01 0.080525533 1.0395578
神泉駅 -5.514785e-02 0.021029193 0.9800615
笹塚駅 4.172115e-01 -0.168983715 0.7900486
各駅の弾力性を地図にプロットしてみましょう。
笹塚の弾力性が低く、北参道の弾力性が高いですね。
北参道駅と笹塚駅が対象的なのでプロットして比較してみたいと思います。
北参道は弾力性が高いので専有面積が広くなるにつれ笹塚の賃料と差が開いていく傾向にある。安く広い賃貸に住みたいなら笹塚がよさそうです。
※ 北参道駅が最寄りの賃貸は最低18㎡〜だったのでx軸の値は18から描画しています。
■ まとめ・今後やりたいこと
1つのモデルで渋谷区全体の賃貸を予測・説明することができるようになったので、上で触れたように駅ごとにどのような賃料の傾向があるのかもわかりやすくなってどの地域に住みたいか判断材料にもなりそうですね!
今回専有面積のみ説明変数に使用したので、次は他の変数を足すなどして決定係数が 0.9 以上になるような精度の高いモデルが作れそうかモデル比較をしてみたい...!
この記事が少しでも階層線形モデルをRでやってみたいな、という方の役に立ったらうれしいです 😌