はじめに

賃貸マンションの家賃相場は、おおむね、最寄り駅、面積、築年数、階数と、駅からの徒歩分数で決まります。これらの要因から家賃相場を推定する階層ベイズモデルを構築しました。なぜこんなことをしているのかというと、部屋探しの過程で家賃が決まるメカニズムに興味を持ったことがきっかけです。

家賃データには誤入力や並外れた高額物件などによる外れ値が含まれます。誤差項に正規分布を仮定した通常のモデルでは、これらの外れ値に推定結果が引きずられてしまいます。そのためこの記事では、誤差項にt分布を用いることで、この問題に対処したロバストなモデルを作りました。

2024年12月にSUUMOに掲載されていた東京23区の賃貸マンションの家賃データ（62万件）を用いて家賃相場を推定したところ、次の内容が分かりました。

• t分布の自由度は7.5程度
  • 正規分布ではなく裾の重さを考慮することが妥当
• 家賃に与える影響
  • 築年数1年あたり1.2%下がる
  • 徒歩1分あたり0.8%下がる
  • 階数が1階上がると0.9%上がる
  • 1階と地下1階は2階からみてそれぞれ4.6%、6.4%下がる
  • 最上階でも家賃相場は変わらない
• t分布モデルでは、急行停車駅が各駅停車駅よりも家賃が高い現象や、都心に近い駅ほど家賃が高い現象をとらえることができた
• t分布モデルは正規分布モデルよりもWAICでの汎化性能が高かった

データ

SUUMOから2024年12月にPythonでスクレイピングした、東京23区の賃貸マンションの家賃相場データ（約62万件）を用います。

ただし、対象はマンションのみ、間取りは1K, 1DK, 1LDK, 2K, 2DK, 2LDK, 家賃+管理費が100万円以下、面積は20m2～100m2、築40年以内、最寄り駅徒歩20分以内、地上15階以下です¹。

このようなデータフレームです。

外れ値の存在

x軸に面積の対数、y軸に家賃+管理費（以下家賃と呼びます）の対数をとり、築年数で色分けしてプロットしてみました。

面積の対数と家賃の対数はほぼ線形に並んでいるため、面積の対数と家賃の対数を線形回帰することは妥当だと思われますが、一部外れ値のような点があります。例えばexp(3.25)=25m2, exp(4.5)=90万円は明らかな誤入力ですね（家賃を1桁多く入力しているようでした）。

線形のライン状に密集しているエリアの少しだけ上と下にも点があります。これらは他の物件とは並外れた特徴がある物件かもしれません。このような物件が一定数あるということは、裾が正規分布より厚い分布に従っているということです。

線形回帰の誤差項を正規分布より裾が厚いt分布やコーシー分布とすることで、どちらのケースにも対応して外れ値に引っ張られないロバストな回帰を行うことができます。こちらの記事が分かりやすいです。

尤度関数におけるガウス分布とスチューデントのt分布の比較 - suzuzusu日記

散布図を見る限りはそこまで裾が重そうな分布ではないのでコーシー分布とするとやりすぎかと思ったのと、自由度が1のt分布はコーシー分布でありt分布はコーシー分布を包含するため、この記事ではt分布を適用してみます。

モデル

家賃相場は最寄り駅、面積、築年数、最寄り駅からの徒歩分数、階数で決まるとします。これはEDAから妥当な仮定であることが分かっているのと、部屋の方角などの他の特徴量を取得することはスクレイピングの時間的制約から難しいため、これらの特徴量だけを用います。ちなみに、このように家賃などの不動産価格を面積などの属性の関数として表すアプローチをヘドニック法といいます（ヘドニック法については清水・唐渡 (2007) が分かりやすかったです）。

最寄り駅によって家賃の水準と面積に対する家賃の弾力性が異なることを階層ベイズで表現します。これにより、物件数が少ない駅でも、東京23区の全体の平均の傾向を借用できる（「縮約」という）ことからパラメータの推定が安定します。

築年数と徒歩分数、階数の影響は共通とします。駅近の方が築古でも家賃が下がりにくいと一般に言われますが、最も効くのは最寄り駅別の家賃水準の違いなので、こちらの影響は共通とします。もちろん高度化の余地はあります。

散布図で見たとおり、家賃と面積は両対数線形の関係とします。また、築年数、徒歩分数、階数の効果は家賃に対して乗算で効くというドメイン的に自然な仮定を置きます（例えば、築1年増えるごとに1%下がるようなイメージです）。これによって、対数を取ると線形モデルとなり、扱いやすくなります。線形モデルの誤差項は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くらいの定数と決め打ちしてもいいと思いますが、自由度もパラメータとして推定することにします。

なお、モデル比較用に別途、t分布のところを正規分布としたモデルも推定しました。

実装

先のモデルを以下の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
  );
}

cmdstanrのcmdstan_modelのsampleを用いて、chains=4, iter_warmup=1000, iter_sampling=1000, thin=1でサンプリングしました。4chain並列で4日かかりました。環境はR=4.4.2, CmdStan=2.35.0, cmdstanr=0.8.1です。

結果

パラメータ

fit$print(
  c("age", "walk", "floor_num", "is_top", "is_ground", "is_underground", "sigma", "nu"),
  max_rows=12, digits=3
)
#>        variable   mean median    sd   mad     q5    q95  rhat ess_bulk ess_tail
#>  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           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%下がる
• 最上階でも家賃相場は変わらない

これらは比較用の正規分布モデルと大きく変わりませんでした。

階数のところですが、例えば4階建てのマンションで2階が家賃10万円なら、3階は10.09万円、4階は10.18万円、1階は9.54万円であり、直感的な結果なのではないでしょうか。

最寄り駅ごとの家賃相場

正規分布ではなくt分布を設定したことで、最寄り駅ごとの家賃相場をドメイン知識に合った形でうまく推定することができました。その一例を挙げます。

25m2、新築、駅から徒歩5分、3階の物件を仮定して、小田急線の最寄り駅別の家賃相場を求めてみます³。

左のプロットが今回推定したt分布のモデル、右のプロットは比較用に別途推定した正規分布のモデルです。横軸の単位は万円で、黒い点は事後分布の中央値、横の棒は95%信用区間です。

今回のt分布のモデルの方がうまく推定できたと思われるところは二つあります。

• 成城学園前と祖師ヶ谷大蔵
  • 成城学園前は急行が止まり、祖師ヶ谷大蔵は各駅のみの停車駅
  • 正規分布のモデルではほぼ差がないが、t分布モデルでは成城学園前＞祖師ヶ谷大蔵
  • 祖師ヶ谷大蔵は成城学園前より一駅分都心に近いが、成城学園前は急行停車駅であることに加えて高級住宅街であるため、t分布の方の結果が納得できる
• 梅ヶ丘と豪徳寺、千歳船橋と祖師ヶ谷大蔵
  • 4駅とも各駅停車駅
  • 正規分布モデルでは梅ヶ丘＜豪徳寺、千歳船橋＜祖師ヶ谷大蔵。t分布モデルでは、95%有意ではないものの、梅ヶ丘＞豪徳寺、千歳船橋＞祖師ヶ谷大蔵
  • 都心に近い方が家賃が高くなるt分布モデルの方がドメイン知識に整合的

WAICによる汎化誤差の比較

このようなきれいなストーリーが偶然ではないことを知りたいですね。汎化誤差が小さいモデルを選ぶために、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
    );
  }
}

これによりlog_likで対数尤度を持っているCmdStanFitモデルを作っておけば、loo::waic()でWAICが計算できます。

今回使ったデータのうち、渋谷区の物件のデータ（約41000件）を用いてパラメータを推定してWAICを求めたところ、以下のとおりでした。

• t分布（自由度もパラメータとして推定）: -74652.41
• 正規分布: -73071.94

このことからも正規分布よりはt分布の方が汎化性能がよいことが分かります。

おわりに

t分布の導入により、外れ値にロバストなモデルを作ることができました。現実のデータは裾が厚い分布に従うことがあるのでt分布は役立ちますね。

参考文献

• 清水千弘・唐渡広志 (2007).『不動産市場の計量経済分析』朝倉書店.