t分布を用いたロバストな家賃相場の階層ベイズモデリング
はじめに
このブログでは、SUUMOのデータを用いて東京23区の賃貸マンションの家賃相場を推定する階層ベイズモデルの開発に取り組んできました。
家賃データには時々外れ値のようなレコードが観測されます。推定されるモデルのパラメータが外れ値に影響を受けないように、モデルにt分布を導入してみました。実装は引き続きRとStanで行いました。
データ
これまでの記事と同様に、SUUMOから2024年12月にPythonでスクレイピングした東京23区の賃貸マンションの家賃相場データを用います。
こういうデータフレームです。
用いるデータは上の画像に載せた東京23区の物件のデータのうち、以下を満たす約62万件の物件の家賃データです。
- マンションのみ（アパートや一戸建ては対象外）
- 間取りは1K, 1DK, 1LDK, 2K, 2DK, 2LDK
- 家賃+管理費が100万円以下
- 面積は20m2～100m2
- 築40年以内
- 駅から徒歩20分以内
- マンションの階数（各物件の階数ではなく、建物自体の階数）は地上15階まで、かつ地下階はないか地下1階まで
ワンルームは部屋の設備が簡略化されているなど1K以上とは異なる特性があると考えられるため、除外しました。
また3LDK以上の物件は賃貸では物件数が少なく希少であり、2LDK以下の物件とは同一のモデルでとらえられないと思ったのでこちらも除外しました。この2年くらい家賃相場が上がっていますが、特に3LDKは家賃の上がり方が激しいです。もともと3LDKは賃貸の供給が少ないですが、昨今の分譲マンション価格の高騰を受けて分譲マンションの購入と賃貸を天秤にかけた結果賃貸を選択したファミリー層の流入によって需要が増えたことが一つの理由だと聞きます。
外れ値の存在
x軸に面積の対数、y軸に家賃+管理費（以下家賃と呼びます）の対数をとり、築年数で色分けしてプロットしてみました。
面積の対数と家賃の対数はほぼ線形に並んでいます。このことから、面積の対数と家賃の対数を線形回帰することは妥当だと思われますが、一部外れ値のような点があります。
たとえばx=3.25, y=4.5くらいにあるいくつかの点です。面積はexp(3.25)=25m2、家賃はexp(4.5)=90万円くらいということです。
SUUMOの掲載ページを確認したところ、どうやら家賃の桁が間違えて一桁多く入力されているようでした。外れ値の検出手法で除外することもできますが少し面倒です。
また、線形のライン状に密集しているエリアの少しだけ上と下にも点があります。これらは誤記入ではないかもしれませんが、他の物件とは並外れた特徴がある物件かもしれません。
前者のような誤記入データが入っていると回帰直線が外れ値に引きずられます。また、仮に後者のような少し外れた特徴がある物件がそれなりにあるとするならば、現実のデータ生成過程を抽象して表現するのが統計モデリングですから、回帰直線の誤差項は通常の正規分布ではなく、それより裾が厚い確率変数であるというデータ生成過程を仮定することが望ましいです1。
線形回帰の誤差項を正規分布より裾が厚いt分布やコーシー分布とすることで、どちらの可能性にも対応でき、ロバストな回帰を行うことができます。こちらの記事が分かりやすいです。尤度関数におけるガウス分布とスチューデントのt分布の比較 - suzuzusu日記
散布図を見る限りはそこまで裾が重そうな分布ではないので、コーシー分布とするとやりすぎかなというのと、自由度が1のt分布はコーシー分布でありt分布はコーシー分布を包含するため、この記事ではt分布を適用してみます。
モデル
物件$i(1, \dots, N)$の最寄り駅を$sta[i] (1, \dots, S)$とします。
$$ \begin{align*} \log{y_{i}} & \sim student\_t(\nu, \mu_{i}, \sigma) \\\ \mu_{i} &= a_{sta[i]} + b_{sta[i]} \log{\mathrm{area}_{i}} \\\ &+ \beta_{\mathrm{age}} \mathrm{age}_{i} \\\ &+ \beta_{\mathrm{walk}}(\mathrm{walk}_{i} - 1) \\\ &+ \beta_{\mathrm{floor}} \max {(\mathrm{floor}_{i} - 2, 0)} \\\ &+ \beta_{\mathrm{isTop}} \mathrm{isTop}_{i} \\\ &+ \beta_{\mathrm{isGround}} \mathrm{isGround}_{i} \\\ &+ \beta_{\mathrm{isUnderground}} \mathrm{isUnderground}_{i} \\\ a_{sta[i]} & \sim N(a_{all}, \sigma_{a_{all}}) \\\ b_{sta[i]} & \sim N(b_{all}, \sigma_{b_{all}}) \\\ \end{align*} $$
このとき、物件の対数家賃の相場は$\mu_{i}$万円であると考えます。ただし、式中の変数は以下のとおりです。
- $y_{i}$: 家賃+管理費（万円）
- $\mathrm{area}_{i} (20 \leq \mathrm{area}_{i} \leq 100)$: 面積（m2）
- $\mathrm{age}_{i} (= 0, 1, \dots, 40)$: 築年数（新築は0年とする）
- $\mathrm{walk}_{i} (= 1, 2, \dots, 20)$: 最寄り駅からの徒歩分数
- $\mathrm{floor}_{i} (= -1, 1, 2, \dots, 15)$: 物件の階数
- $\mathrm{isTop}_{i} (= 0, 1)$: 最上階なら1, そうではないなら0
- $\mathrm{isGround}_{i} (= 0, 1)$: 1階なら1, そうではないなら0
- $\mathrm{isUnderground}_{i} (= 0, 1)$: 地下1階なら1, そうではないなら0
式中の$\nu$はt分布の自由度です。StanとRでベイズ統計モデリングによれば、そこまで裾が重くない分布の場合は自由度6～8を設定するとよいそうです。なので自由度を6～8くらいの定数と決め打ちしてもいいと思いますが、自由度もパラメータとして推定することにします。
実装
以下のStanコードを”model.stan”というファイル名で保存します。
data {
int N; // 物件の数
vector[N] Y; // 家賃+管理費
vector[N] AREA; // 面積
int S; // 最寄り駅の数
array[N] int<lower=1, upper=S> STATION; // 最寄り駅index
vector[N] AGE; // 築年数（0 - 40）
vector[N] WALK; // 徒歩分数（1 - 20）
vector[N] FLOOR; // 階数（-1, 1 - 15）
vector[N] IS_TOP; // 最上階かどうか（0/1）
vector[N] IS_GROUND; // 1階かどうか（0/1）
vector[N] IS_UNDERGROUND; // 地下1階かどうか（0/1）
}
transformed data {
vector[N] FLOOR2;
for (i in 1:N) {
if (FLOOR[i] <= 1) {
FLOOR2[i] = 2;
} else {
FLOOR2[i] = FLOOR[i];
}
}
}
parameters {
real a_all;
real b_all;
vector[S] a;
vector[S] b;
real<upper=0> age;
real<upper=0> walk;
real<lower=0> floor_num;
real<lower=0> is_top;
real<upper=0> is_ground;
real<upper=0> is_underground;
real<lower=0> sigma_a;
real<lower=0> sigma_b;
real<lower=0> sigma;
real<lower=1> nu;
}
model {
a ~ normal(a_all, sigma_a);
b ~ normal(b_all, sigma_b);
log(Y) ~ student_t(
nu,
a[STATION] + b[STATION] .* log(AREA) +
age*AGE +
walk*(WALK - 1)+
floor_num*(FLOOR2 - 2) +
is_top*IS_TOP +
is_ground*IS_GROUND +
is_underground*IS_UNDERGROUND,
sigma
);
}
次にこちらのRコードでStanコードをキックします。さきほど載せたデータフレームをdf_uniqueという変数で持っている前提です。
library(tidyverse)
library(cmdstanr)
library(bayesplot)
library(tidybayes)
df <- df_unique |>
filter(rent_admin <= 100) |>
filter(area <= 100 & area >= 20) |>
filter(age <= 40) |>
filter(walk <= 20) |>
filter(story_under <= 1L & story_above <= 15L) |>
filter(floor >= -1L & floor <= 15L) |>
filter(layout %in% c("1K", "1DK", "1LDK", "2K", "2DK", "2LDK")) |>
# Stanに入れるために、駅名をintegerに変換する
mutate(station_index=as.integer(as.factor(station))) |>
mutate(
# 平屋や2階建ての2階の場合は「最上階」とはみなさないことにする（その方が直感的に自然なので）
is_top=as.integer(floor == story_above & story_above >= 3),
is_ground=as.integer(floor == 1L),
is_underground=as.integer(floor <= -1L)
)
mod <- cmdstanr::cmdstan_model("model.stan")
fit <- mod$sample(
data=list(
N=nrow(df),
Y=df$rent_admin,
AREA=df$area,
S=length(unique(df$station_index)),
STATION=df$station_index,
AGE=df$age,
WALK=df$walk,
FLOOR=df$floor,
IS_TOP=df$is_top,
IS_GROUND=df$is_ground,
IS_UNDERGROUND=df$is_underground
),
chains=4, parallel_chains=4, iter_warmup=1000, iter_sampling=1000, thin=1,
seed=1234, refresh=10
)
warmupを入れて合計2000回のiterationで約4日かかりました。モデルがすこしずつ複雑になってきたので時間がかかりますね…。
結果
パラメータ
fit$print(
c("a_all", "b_all", "age", "walk", "floor_num", "is_top", "is_ground", "is_underground",
"sigma_a", "sigma_b", "sigma", "nu"),
max_rows=12, digits=3
)
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#> a_all -0.141 -0.141 0.016 0.016 -0.167 -0.114 1.000 4505 3010
#> b_all 0.838 0.838 0.006 0.006 0.829 0.847 1.002 5600 3239
#> age -0.012 -0.012 0.000 0.000 -0.012 -0.012 1.001 3951 2624
#> walk -0.008 -0.008 0.000 0.000 -0.008 -0.008 1.001 7615 3078
#> floor_num 0.009 0.009 0.000 0.000 0.009 0.009 1.000 6234 2861
#> is_top 0.000 0.000 0.000 0.000 0.000 0.000 1.001 3665 2064
#> is_ground -0.049 -0.049 0.000 0.000 -0.049 -0.048 1.000 7561 2916
#> is_underground -0.066 -0.066 0.002 0.002 -0.070 -0.062 1.002 7136 2928
#> sigma_a 0.350 0.350 0.012 0.012 0.331 0.370 1.000 4787 2908
#> sigma_b 0.124 0.124 0.004 0.004 0.117 0.131 1.001 6393 2580
#> sigma 0.093 0.093 0.000 0.000 0.093 0.094 1.000 6458 3222
#> nu 7.562 7.561 0.076 0.074 7.440 7.689 1.000 7036 3119
自由度$\nu$の推定値は7.56でした。そこまで裾が重くはないですが、正規分布とはいえないくらいには裾が重い分布ということを示します。
築年数効果（age）、徒歩分数効果（walk）、階数効果（floor_num）、最上階効果（is_top）、1階効果（is_ground）、地下1階効果（is_underground）の各パラメータより、以下のことが分かりました。
- 築年数が1年増えるごとに家賃相場は1.2%下がる
- 徒歩1分増えるごとに0.8%下がる
- 2階から上に1階高くなるごとに0.9%上がる
- 1階と地下1階は2階から見てそれぞれ4.6%、6.4%下がる
- 最上階でも家賃相場は変わらない
これらは正規分布を用いた場合と大きく変わりませんでした。
最寄り駅ごとの家賃相場
正規分布ではなくt分布を設定したことで、最寄り駅ごとの家賃相場をよりうまく推定することができたように思われます。その一例を挙げます。
25m2、新築、駅から徒歩5分、3階の物件を仮定して、小田急線の最寄り駅別の家賃相場を求めてみます2。
左のプロットが今回推定したt分布のモデル、右のプロットは前回の記事で推定した正規分布のモデルです。横軸の単位は万円で、黒い点は事後分布の中央値、横の棒は95%信用区間です。
今回のt分布のモデルの方がうまく推定できたと思われるところは二つあります。
まず、成城学園前と祖師ヶ谷大蔵を比べると、正規分布のモデルの方はほぼ差がありませんが、t分布の方はそれよりは差が開いています。成城学園前は急行が止まり、祖師ヶ谷大蔵は各駅のみの停車駅です。急行停車駅とその一つ隣の各駅停車駅では、ふつう前者の方が相場が高くなるため、t分布の方の結果は納得できます。
また、t分布の方は95%有意とはいえないものの豪徳寺＞梅ヶ丘、千歳船橋＞祖師ヶ谷大蔵ですが、正規分布の方は豪徳寺＜梅ヶ丘、祖師ヶ谷大蔵＜千歳船橋です。この4駅はすべて各駅のみ止まるため、新宿に近い方が家賃が高くなるt分布の結果の方がドメイン知識に合っていますね。
とはいえ、このようなきれいなストーリーはたまたまかもしれません。ちゃんと判断するために、汎化誤差の小さいモデルはどちらなのかをみてみます。
汎化誤差が小さいモデルはWAICが小さいモデルです。WAICを算出するには、まずStanのコードのmodelブロックの下に以下のgenerated quantitiesブロックを追加して対数尤度を計算しておきます。
generated quantities {
vector[N] log_lik; // 対数尤度
for (i in 1:N) {
log_lik[i] = student_t_lpdf(
log(Y[i]) |
nu,
a[STATION[i]] + b[STATION[i]] .* log(AREA[i]) +
age*AGE[i] +
walk*(WALK[i] - 1)+
floor_num*(FLOOR2[i] - 2) +
is_top*IS_TOP[i] +
is_ground*IS_GROUND[i] +
is_underground*IS_UNDERGROUND[i],
sigma
);
}
}
次にcmdstanrのCmdStanFitオブジェクトを以下の関数の引数で渡すことでWAICを計算できます。
この関数のコードはcmdstanrで推定後にWAICを計算 - to be continued…をそのまま参考にしました。
# fitはcmdstanrオブジェクト
# generated quantitiesブロックで`log_lik`という変数名で対数尤度を計算している前提
waic <- function(fit) {
log_lik <- fit$draws("log_lik", format="df") |>
select(contains("log_lik"))
lppd <- sum(log(colMeans(exp(log_lik))))
p_waic <- sum(apply(log_lik, 2, var))
waic <- -2*lppd + 2*p_waic
return(list(waic=waic, lppd=lppd, p_waic=p_waic))
}
waic(fit)
今回使ったデータのうち、渋谷区の物件のデータのみ（約41000件）を用いてパラメータを推定してWAICを求めたところ、以下のとおりでした。
- t分布（自由度もパラメータとして推定）: -74652.41
- ちなみにこのとき自由度は7.88（90%CI: 7.46-8.33）だった
- 正規分布: -73071.94
このことからも正規分布よりはt分布の方が汎化性能がよいことが分かります。なお、渋谷区のデータのみでWAICを計算したのは、全てのデータで対数尤度を計算するとメモリに乗り切らなくなるためです。
おわりに
t分布の導入により外れ値にロバストなモデルを作ることができました。現実のデータは裾が厚い分布に従うことがあるのでt分布は役立ちますね。