Stan:階層モデル②
過分散
今回は緑本(データ解析のための統計モデリング入門)で紹介されている階層モデルの導入部分を実践する。
データとしては以下を用意する。
x [1] 1 8 9 2 4 3 2 9 5 8
緑本では植物の種子数だったが、これは別に何でもよい。
今回は10点満点のテストの生徒の点数などと考える。
これデータの平均と分散を最尤法で求めると以下となる。
n <- length(x) opt <- optimize(l, interval = c(0, 1), maximum = TRUE) p <- opt$maximum p # mean [1] 0.5100017 n*p*(1-p) # var [1] 2.499
optimize
は最適解を計算してくれる関数。
しかし、実際の平均と分散を求めると以下となる。
mean(x) [1] 5.1 var(x) [1] 9.877778
最尤法で求めた分散と大きく乖離していることが分かる。
(これを「過分散」と言う。)
データ全体は二項分布に従っているが、個人差の影響が大きすぎるためこのようなことが起きてしまっている。
階層モデル
では、Stanで個人差を考慮したモデルを作る。
詳しいことはTJOさんのブログが詳しい。
tjo.hatenablog.com
Stanコードは以下の通りである。
data { int<lower=0> N; int<lower=0> x[N]; } parameters { real beta; real r[N]; real<lower=0> s; } transformed parameters { real q[N]; for (i in 1:N) q[i] = inv_logit(beta+r[i]); } model { for (i in 1:N) x[i] ~ binomial(N,q[i]); for (i in 1:N) r[i]~normal(0,s); }
推定結果は以下のとおりである。
model <- stan_model("hierarchial.stan") fit <- sampling(model,data=list(N=length(x),x=x)) fit Inference for Stan model: hierarchial. 4 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta 0.07 0.02 0.63 -1.16 -0.32 0.06 0.45 1.37 889 1.01 r[1] -1.84 0.02 1.01 -4.05 -2.44 -1.77 -1.15 -0.07 1923 1.00 r[2] 1.16 0.02 0.93 -0.57 0.55 1.12 1.74 3.16 1486 1.00 r[3] 1.74 0.02 1.00 0.01 1.06 1.68 2.34 3.93 2059 1.00 r[4] -1.29 0.02 0.91 -3.23 -1.83 -1.24 -0.69 0.38 1592 1.00 r[5] -0.43 0.02 0.82 -2.08 -0.93 -0.41 0.10 1.18 1313 1.00 r[6] -0.81 0.02 0.86 -2.64 -1.35 -0.77 -0.24 0.80 1356 1.01 r[7] -1.28 0.02 0.90 -3.24 -1.81 -1.24 -0.67 0.33 1708 1.00 r[8] 1.72 0.02 1.00 -0.01 1.05 1.63 2.32 3.90 1839 1.00 r[9] -0.06 0.02 0.82 -1.70 -0.59 -0.06 0.47 1.53 1284 1.00 r[10] 1.15 0.02 0.89 -0.52 0.55 1.11 1.71 2.96 1619 1.00 s 1.75 0.01 0.64 0.83 1.30 1.64 2.08 3.28 1835 1.00 q[1] 0.18 0.00 0.11 0.02 0.09 0.16 0.24 0.44 4000 1.00 q[2] 0.75 0.00 0.13 0.48 0.67 0.76 0.85 0.94 4000 1.00 q[3] 0.83 0.00 0.10 0.58 0.77 0.85 0.91 0.98 4000 1.00 q[4] 0.25 0.00 0.12 0.05 0.16 0.24 0.33 0.51 4000 1.00 q[5] 0.42 0.00 0.14 0.17 0.32 0.41 0.52 0.69 4000 1.00 q[6] 0.34 0.00 0.13 0.11 0.24 0.33 0.43 0.63 4000 1.00 q[7] 0.25 0.00 0.12 0.06 0.16 0.24 0.33 0.53 4000 1.00 q[8] 0.83 0.00 0.11 0.57 0.77 0.85 0.91 0.98 4000 1.00 q[9] 0.50 0.00 0.14 0.24 0.40 0.50 0.60 0.77 4000 1.00 q[10] 0.75 0.00 0.12 0.48 0.67 0.76 0.84 0.95 4000 1.00 lp__ -63.37 0.08 2.80 -69.92 -65.00 -63.00 -61.33 -58.96 1340 1.00
これだけではよくわからない。
個人の二項分布の事象の生起確率は以下より、と分かる。
これに試行回数を掛けると、各個人の確率分布の平均が分かる。
これより、各生徒の二項分布の平均は
「1.8, 7.5, 8.3, 2.5, 4.2, 3.4, 2.5, 8.3, 5.0, 7.5」だと分かる。
実際の点数は
x [1] 1 8 9 2 4 3 2 9 5 8
だったので、割といい線いっていることが分かる。