Stan で階層ベイズモデルの構築およびパラメーターの推測

階層モデル(単回帰)

サンプルデータの作成

説明変数だけでは説明がつかないような群間差や個体差を説明するには、階層ベイズモデルが使われる。階層ベイズモデルのパラメーターを Stan で推定するため、ここで乱数を使って、次の図にあるような構造からなる 3 つのグループからなるサンプルデータを作成する。

単回帰階層ベイズモデル
library(ggplot2)
library(ggsci)

group.1.x <- runif(100, 10, 60)
group.2.x <- runif(100, 20, 80)
group.3.x <- runif(100, 40, 80)
group.1.y <- rnorm(100, group.1.x * rnorm(100, 2, 0.2) + rnorm(100, 30, 2), 10)
group.2.y <- rnorm(100, group.2.x * rnorm(100, 3, 0.2) + rnorm(100, 50, 2), 10)
group.3.y <- rnorm(100, group.3.x * rnorm(100, 4, 0.2) + rnorm(100, 20, 2), 10)

d <- data.frame(x = c(group.1.x, group.2.x, group.3.x),
                y = c(group.1.y, group.2.y, group.3.y),
                group = rep(c(1, 2, 3), each = 100))

g <- ggplot(d, aes(x = x, y = y, color = factor(group)))
g <- g + geom_point() + scale_color_nejm()
print(g)
単回帰階層ベイズモデル(サンプルデータの図示)

Stan モデル

次に、Stan コードでモデルを記述していく。このサンプルデータでは 3 つの実験群があるので、実験群間の差を表すためのパラメーターを導入する。まず、群間効果の平均値を *_ave を導入し、その平均値との差 *_diff を各実験群の群差とおく。これらの変数を導入して、群間効果を考慮したモデルは次のようになる。また、解析後に、予測結果も合わせて図示するため、generated quantities ブロックで、推定されたパラメーターを使って再サンプリングを行うようにする。

data {
    int G;
    int N;
    real x[N];
    real y[N];
    int<lower=1,upper=G> group[N];
}
parameters {
    // group averages of parameters beta0 and beta1
    real beta0_ave;
    real beta1_ave;

    // differences from group averages for each group
    real beta0_diff[G];
    real beta1_diff[G];

    // standard deviation for each parameters
    real<lower=0> sigma_beta0;
    real<lower=0> sigma_beta1;
    real<lower=0> sigma_y;
}
transformed parameters {
    real beta0[G];
    real beta1[G];
    for (g in 1:G) {
        beta0[g] = beta0_ave + beta0_diff[g];
        beta1[g] = beta1_ave + beta1_diff[g];
    }
}
model {
    for (g in 1:G) {
        beta0_diff[g] ~ normal(0, sigma_beta0);
        beta1_diff[g] ~ normal(0, sigma_beta1);
    }
    for (n in 1:N) {
        y[n] ~ normal(beta0[group[n]] + beta1[group[n]] * x[n], sigma_y);
    }
}
generated quantities {
    real yhat[N];
    for (n in 1:N) {
        yhat[n] =  normal_rng(beta0[group[n]] + beta1[group[n]] * x[n], sigma_y);
    }
}

パラメーター推定

次に rstan パッケージを利用して、パラメーター推定を行う。

library(rstan)

d <- list(x = d$x, y = d$y, group = d$group,
          N = length(d$x), G = length(unique(d$group)))

fit <- stan(file = 'lb.stan', data = d)

fit オブジェクトには様々な情報が含まれているので、そのまま表示すると分かりづらくなるので、ここでは 3 つの実験群のパラメーターである β0 および β1 の 95% 信頼区間およびその平均値を求めて表示させる。

ms <- rstan::extract(fit, pars = c('beta0', 'beta1'))

apply(ms$beta0, 2, quantile, probs = c(0.025, 0.500, 0.975))
##             [,1]     [,2]     [,3]
##   2.5%  28.50130 40.13231 15.21529
##   50%   36.35821 49.97955 31.75733
##   97.5% 43.88187 59.64620 46.71032

apply(ms$beta1, 2, quantile, probs = c(0.025, 0.500, 0.975))
##             [,1]     [,2]     [,3]
##   2.5%  1.668467 2.754701 3.546323
##   50%   1.845419 2.943658 3.790877
##   97.5% 2.036375 3.132323 4.059905

推定結果と乱数を発生させた時の設定値を見比べると、3 つの実験群における β0 の推定値は 36.4, 50.0、31.8 であり、真の値は 30, 50, 20 である。また、3 つの実験群における β0 の推定値は 1.8、2.9、3.8 であり、真の値は 2, 3, 4 である。この結果をみると、β0 の推定は大きく外れているが、β1 の推定値は真の値と非常に近いことがわかる。

stan_hist(fit, pars = c('beta0', 'beta1'))
単回帰階層ベイズモデルのパラメーター推定結果(ヒストグラム)

次に Stan コードの generated quantities のブロックで再サンプリングした標本を用いて、実験群ごとに予測区間を描く。

ms <- rstan::extract(fit, pars = 'yhat')

df.pred <- data.frame(x = d$x, group = factor(d$group),
                      lower  = apply(ms$yhat, 2, quantile, prob = 0.025),
                      median = apply(ms$yhat, 2, quantile, prob = 0.500),
                      upper  = apply(ms$yhat, 2, quantile, prob = 0.925))

g <- ggplot(df.pred, aes(x = x))
g <- g + geom_ribbon(aes(ymin = lower, ymax = upper, fill = group), alpha = 0.4)
g <- g + geom_line(aes(y = median, color = group))
g <- g + geom_point(data = data.frame(x = d$x, y = d$y, group = factor(d$group)),
                    aes(x = x, y = y, colour = group))
g <- g + scale_color_nejm() + scale_fill_nejm()
g <- g + xlab('x') + ylab('y')
print(g)
Stan で推定した階層ベイズモデルの予測区間を ggplot で描く方法