階層ベイズを用いた地域ごとの家賃の変化の推測

目的

RStanによる階層ベイズを用いて, 日本の各地域ごとの家賃の変化を推測する.

今回は, e-statの日本の様々な町ごとの1991年と2000年の家賃のデータを用いて, 各地域の家賃の変化の推測を行う. その際, 日本全体で大まかな変化の傾向があることを仮定して, 階層ベイズモデルを適用する.

具体的には, データに含まれる67箇所の都市を6つの地域(北海道&東北, 関東, 中部, 関西, 中国&四国, 九州)に分類して, 次のようなモデルにより線形回帰を行う.

$Y[n] \sim {\rm Normal}(a[area(n)] + b[area(n)] * X[n], \sigma_Y[area(n)]) \\ a[k] \sim {\rm Normal}(a_{ave}, \sigma_a) (k = 1, 2, ...,6) \\ b[k] \sim {\rm Normal}(b_{ave}, \sigma_b) (k = 1, 2, ...,6) \\ \sigma_Y[k] \sim \Gamma(1.0, beta0) (k = 1, 2, ...,6)$

(ただし $area(n)$ はサンプル点 $n$ の所属する地域を表す)

このモデルでは, 地域 $k$ のX(1991年の家賃)とY(2000年の家賃)の関係を示す線形モデルのパラメータ $a[k]$ が, 平均 $a_{ave}$ , 標準偏差 $\sigma_a$ の正規分布から生成されていると考える( $b[k]$ も同様に平均 $b_{ave}$ , 標準偏差 $\sigma_b$ の正規分布から生成されているとする). また, 標準偏差 $\sigma_Y[k]$ は, $k$ ごとに定数1.0と $beta0$ をパラメータとしたガンマ分布から生成しているとする. こうした階層的なモデルを用いることで, それぞれの地域ごとの少ないサンプルで線形回帰した場合に起こる過学習を防げる可能性がある. また. 他の地方の情報が, 事前分布を通じて, 別地方の推定に活用できる.

上記のモデルをStanで定義し, 各パラメータの事後分布のサンプリングを行うコードは次の通りである.

library(rstan)

options(mc.cores = parallel::detectCores())
d <- read.csv("yachin.csv")
area <- d$X[2:68]
town <- d$X.1[2:68]
rent1991 <- d$X.2[2:68] / 1000
rent2000 <- d$X.3[2:68] / 1000
df <- data.frame(town, rent1991, rent2000, area, stringsAsFactors=FALSE)

data <- list(N=nrow(df), K=6, X=df$rent1991, Y=df$rent2000, KID=df$area)
fit <- stan(file="kaisou.stan", data=data, seed=123)

Stan(kaisou.stan) ... StanとRでベイズ統計モデリング(松浦健太郎, 共立出版) p137のコードを元に, 標準偏差 $\sigma_Y[k]$ も地域ごとにガンマ分布から生成するよう変更

data {
    int N;
    int K;
    real X[N];
    real Y[N];
    int<lower=1, upper=K> KID[N];
}

parameters {
    real a0;
    real b0;
    real<lower=0> beta0;
    real a[K];
    real b[K];
    real<lower=0> s_a;
    real<lower=0> s_b;
    real<lower=0> s_Y[K];
}

model {
    for (k in 1:K) {
        a[k] ~ normal(a0, s_a);
        b[k] ~ normal(b0, s_b);
        s_Y[k] ~ gamma(1.0, beta0);
    }

    for (n in 1:N) {
        Y[n] ~ normal(a[KID[n]] + b[KID[n]] * X[n], s_Y[KID[n]]);
    }
}

generated quantities {
    vector[N] log_lik;

    for (n in 1:N) {
        log_lik[n] = normal_lpdf(Y[n] | a[KID[n]] + b[KID[n]] * X[n], s_Y[KID[n]]);
    }
}

なお, 元データの家賃は民営借家一ヶ月( $3.3m^2$ あたり)の値であり, 前処理の段階で扱いやすくするために1/1000倍している.

これとは別に, 各地域ごとに別々に線形回帰を行った. 地域ごとに行った回帰により予測された線を黒, 階層ベイズを用いた回帰により予測された線を赤でプロットすると, 次の図ようになった.

f:id:gakushukun1:20190604205214p:plain

(グラフ中の番号は , 1...北海道&東北, 2...関東, 3...中部, 4...関西, 5...中国&四国, 6...九州に該当, また点は各地方に属している都市のデータを表す)

グラフを見ると, 2(関東)のように２つのモデルで予測された線がほぼ同一になるものもあれば, 5(中国&四国)のようにまったく異なる線になる場合があると分かる. 特に5の場合, 家賃が同じ価格帯の都市のデータばかりであったために, 地域単体の回帰では, 日本全体の傾向(傾きが1.00に近い=1991年から2000年にかけて家賃が横ばい)から外れたものになっていると考えられる.

地域ごとのデータで行った回帰と, 階層ベイズによる回帰の汎化性能を比べるため, それぞれの場合での予測分布に対するWAICとLOOの値を次に示す. 値が小さい方が汎化性能が優れているといえる.

	WAIC(個々)	WAIC(階層)	LOO(個々)	LOO(階層)
1.北海道&東北	11.9	10.1	12.2	10.2
2.関東	25.9	24.1	26.1	24.2
3.中部	29.7	26.7	30.6	27.1
4.関西	22.4	20.9	22.8	21.5
5.中国&四国	12.4	14.3	12.8	14.6
6.九州	15.7	12.8	15.9	13.0