水槽ごとにカエルの卵の生存率を推定する。 階層ベイズモデルにおける縮退(Shrinkage)効果を確認する。
Richard McElreath, "Statistical Rethinking", Chapter12.1. Example: Multilevel tadpoles
今回は水槽ごとのdensity(卵の数=水槽の大きさ)とsurv(生存数)の二つの変数を使用する.
各水槽iにはDensity[i]個の卵があり, そのうちSurv[i]個が生き残る. このとき, 生存する卵は二項分布Binomial(Density[i], p[i])に従うと考える. ここで, p[i]は水槽iにおける生存確率である. 生存確率p[i]は水槽によって異なる係数αpond[i] をロジスティック関数${\rm inv_logit}$で(0, 1)の区間に変換したものとする. αpond[i]は平均α, 標準偏差σの正規分布から生成されるものとする. また, αの事前分布は平均0, 標準偏差1の正規分布, σの事前分布は平均0, 標準偏差1の半コーシー分布とする.
$$
\begin{align}
Surv[i] & \sim {\rm Binomial}(Density[i], p[i]) &i&=1,2,...,I\\
p[i] & = {\rm inv\_logit}(\alpha_{pond}[i]) &i&=1,2,...,I\\
\alpha_{pond}[i] & \sim {\rm Normal}(\alpha, \sigma) &i&=1,2,...,I\\
\alpha & \sim {\rm Normal}(0, 1)\\
\sigma & \sim {\rm HalfCauchy}(0, 1)\\
\end{align}
$$
ms <-rstan::extract(fit)propotion_survival <- reedfrogs$surv / reedfrogs$density
average_propotion <- sum(reedfrogs$surv)/sum(reedfrogs$density)
plot_df <- data.frame(pond_index=1:48, empirical=propotion_survival, estimated=colMeans(ms$p))
plot_df.long <- melt(plot_df, id.vars="pond_index", variable.name="kind", value.name="propotion")ggplot(plot_df.long, aes(x=pond_index, y=propotion, colour=kind)) +
geom_point() +
geom_hline(yintercept=average_propotion, linetype="dashed") +
xlab("Pond Index") +
ylab("Propotion Survival")
colMeans(ms$p)## [1] 0.8702066 0.9323779 0.7144771 0.9321336 0.8689858 0.8706553 0.9307631
## [8] 0.8665497 0.4588228 0.8673171 0.7079410 0.6241969 0.7125905 0.5403379
## [15] 0.8649620 0.8713687 0.9351615 0.9003667 0.8686281 0.9623659 0.9034832
## [22] 0.9029950 0.9022850 0.8296456 0.2732475 0.5365202 0.2000485 0.3881909
## [29] 0.5388969 0.7968774 0.3506360 0.4249060 0.9508013 0.9263799 0.9259601
## [36] 0.8762601 0.8761870 0.9701417 0.9307551 0.9025013 0.1494744 0.3638514
## [43] 0.3924659 0.4178878 0.6398540 0.3600154 0.8790053 0.4968829ggplot()+theme_set(theme_bw(base_size=18,base_family="HiraKakuProN-W3"))
print(fit)## Inference for Stan model: Shrinkage.
## 4 chains, each with iter=500; warmup=250; thin=1;
## post-warmup draws per chain=250, total post-warmup draws=1000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## alpha_pond[1] 2.14 0.03 0.85 0.70 1.55 2.07 2.65
## alpha_pond[2] 3.08 0.06 1.16 1.24 2.29 2.95 3.76
## alpha_pond[3] 1.00 0.02 0.63 -0.24 0.59 0.98 1.43
## alpha_pond[4] 3.05 0.03 1.10 1.31 2.24 2.98 3.77
## alpha_pond[5] 2.14 0.03 0.88 0.64 1.50 2.10 2.67
## alpha_pond[6] 2.15 0.03 0.85 0.69 1.54 2.05 2.68
## alpha_pond[7] 3.10 0.04 1.19 0.96 2.29 2.95 3.82
## alpha_pond[8] 2.13 0.03 0.89 0.56 1.55 2.06 2.67
## alpha_pond[9] -0.18 0.02 0.61 -1.39 -0.59 -0.17 0.25
## alpha_pond[10] 2.14 0.03 0.91 0.53 1.55 2.04 2.67
## alpha_pond[11] 0.97 0.02 0.64 -0.26 0.52 0.96 1.37
## alpha_pond[12] 0.55 0.02 0.62 -0.64 0.15 0.54 0.92
## alpha_pond[13] 1.00 0.02 0.68 -0.25 0.51 0.98 1.48
## alpha_pond[14] 0.18 0.02 0.62 -1.04 -0.21 0.18 0.58
## alpha_pond[15] 2.13 0.03 0.92 0.48 1.52 2.09 2.70
## alpha_pond[16] 2.20 0.03 0.94 0.61 1.51 2.10 2.80
## alpha_pond[17] 2.93 0.03 0.83 1.51 2.34 2.89 3.42
## alpha_pond[18] 2.38 0.02 0.71 1.20 1.87 2.31 2.79
## alpha_pond[19] 2.01 0.02 0.59 1.00 1.60 1.98 2.34
## alpha_pond[20] 3.66 0.03 1.04 1.97 2.90 3.55 4.32
## alpha_pond[21] 2.41 0.02 0.68 1.24 1.93 2.35 2.82
## alpha_pond[22] 2.40 0.02 0.67 1.20 1.94 2.34 2.80
## alpha_pond[23] 2.38 0.02 0.64 1.24 1.92 2.33 2.79
## alpha_pond[24] 1.68 0.02 0.57 0.67 1.29 1.65 2.06
## alpha_pond[25] -1.02 0.01 0.44 -1.89 -1.32 -0.99 -0.71
## alpha_pond[26] 0.15 0.01 0.39 -0.58 -0.12 0.15 0.42
## alpha_pond[27] -1.45 0.02 0.48 -2.43 -1.75 -1.41 -1.13
## alpha_pond[28] -0.47 0.01 0.41 -1.34 -0.72 -0.48 -0.19
## alpha_pond[29] 0.16 0.01 0.39 -0.61 -0.10 0.16 0.41
## alpha_pond[30] 1.44 0.02 0.51 0.52 1.09 1.41 1.78
## alpha_pond[31] -0.64 0.01 0.42 -1.46 -0.92 -0.62 -0.36
## alpha_pond[32] -0.31 0.01 0.38 -1.05 -0.58 -0.30 -0.05
## alpha_pond[33] 3.16 0.02 0.71 2.00 2.65 3.11 3.62
## alpha_pond[34] 2.70 0.02 0.65 1.56 2.25 2.66 3.08
## alpha_pond[35] 2.68 0.02 0.63 1.55 2.23 2.67 3.09
## alpha_pond[36] 2.04 0.02 0.49 1.12 1.71 2.03 2.36
## alpha_pond[37] 2.06 0.02 0.53 1.11 1.68 2.04 2.40
## alpha_pond[38] 3.87 0.03 1.01 2.23 3.15 3.75 4.54
## alpha_pond[39] 2.76 0.02 0.64 1.66 2.30 2.72 3.15
## alpha_pond[40] 2.34 0.02 0.56 1.34 1.95 2.31 2.71
## alpha_pond[41] -1.82 0.02 0.51 -2.98 -2.13 -1.78 -1.48
## alpha_pond[42] -0.58 0.01 0.35 -1.31 -0.81 -0.57 -0.33
## alpha_pond[43] -0.45 0.01 0.34 -1.16 -0.68 -0.43 -0.23
## alpha_pond[44] -0.34 0.01 0.37 -1.05 -0.60 -0.34 -0.10
## alpha_pond[45] 0.59 0.01 0.34 -0.06 0.37 0.59 0.82
## alpha_pond[46] -0.59 0.01 0.34 -1.29 -0.81 -0.59 -0.36
## alpha_pond[47] 2.09 0.02 0.54 1.13 1.69 2.07 2.43
## alpha_pond[48] -0.01 0.01 0.33 -0.65 -0.22 -0.02 0.20
## alpha 1.30 0.01 0.26 0.82 1.12 1.29 1.47
## sigma 1.62 0.01 0.21 1.25 1.48 1.60 1.75
## p[1] 0.87 0.00 0.09 0.67 0.82 0.89 0.93
## p[2] 0.93 0.00 0.06 0.78 0.91 0.95 0.98
## p[3] 0.71 0.00 0.12 0.44 0.64 0.73 0.81
## p[4] 0.93 0.00 0.06 0.79 0.90 0.95 0.98
## p[5] 0.87 0.00 0.09 0.65 0.82 0.89 0.93
## p[6] 0.87 0.00 0.09 0.67 0.82 0.89 0.94
## p[7] 0.93 0.00 0.07 0.72 0.91 0.95 0.98
## p[8] 0.87 0.00 0.09 0.64 0.83 0.89 0.94
## p[9] 0.46 0.00 0.14 0.20 0.36 0.46 0.56
## p[10] 0.87 0.00 0.09 0.63 0.82 0.89 0.94
## p[11] 0.71 0.00 0.12 0.44 0.63 0.72 0.80
## p[12] 0.62 0.00 0.13 0.35 0.54 0.63 0.72
## p[13] 0.71 0.00 0.13 0.44 0.63 0.73 0.81
## p[14] 0.54 0.00 0.14 0.26 0.45 0.54 0.64
## p[15] 0.86 0.00 0.10 0.62 0.82 0.89 0.94
## p[16] 0.87 0.00 0.09 0.65 0.82 0.89 0.94
## p[17] 0.94 0.00 0.05 0.82 0.91 0.95 0.97
## p[18] 0.90 0.00 0.06 0.77 0.87 0.91 0.94
## p[19] 0.87 0.00 0.06 0.73 0.83 0.88 0.91
## p[20] 0.96 0.00 0.03 0.88 0.95 0.97 0.99
## p[21] 0.90 0.00 0.05 0.78 0.87 0.91 0.94
## p[22] 0.90 0.00 0.05 0.77 0.87 0.91 0.94
## p[23] 0.90 0.00 0.05 0.77 0.87 0.91 0.94
## p[24] 0.83 0.00 0.07 0.66 0.78 0.84 0.89
## p[25] 0.27 0.00 0.08 0.13 0.21 0.27 0.33
## p[26] 0.54 0.00 0.09 0.36 0.47 0.54 0.60
## p[27] 0.20 0.00 0.07 0.08 0.15 0.20 0.24
## p[28] 0.39 0.00 0.09 0.21 0.33 0.38 0.45
## p[29] 0.54 0.00 0.09 0.35 0.48 0.54 0.60
## p[30] 0.80 0.00 0.08 0.63 0.75 0.80 0.86
## p[31] 0.35 0.00 0.09 0.19 0.29 0.35 0.41
## p[32] 0.42 0.00 0.09 0.26 0.36 0.43 0.49
## p[33] 0.95 0.00 0.03 0.88 0.93 0.96 0.97
## p[34] 0.93 0.00 0.04 0.83 0.90 0.93 0.96
## p[35] 0.93 0.00 0.04 0.83 0.90 0.94 0.96
## p[36] 0.88 0.00 0.05 0.75 0.85 0.88 0.91
## p[37] 0.88 0.00 0.05 0.75 0.84 0.88 0.92
## p[38] 0.97 0.00 0.03 0.90 0.96 0.98 0.99
## p[39] 0.93 0.00 0.04 0.84 0.91 0.94 0.96
## p[40] 0.90 0.00 0.05 0.79 0.88 0.91 0.94
## p[41] 0.15 0.00 0.06 0.05 0.11 0.14 0.18
## p[42] 0.36 0.00 0.08 0.21 0.31 0.36 0.42
## p[43] 0.39 0.00 0.08 0.24 0.34 0.39 0.44
## p[44] 0.42 0.00 0.09 0.26 0.35 0.42 0.48
## p[45] 0.64 0.00 0.08 0.49 0.59 0.64 0.69
## p[46] 0.36 0.00 0.08 0.22 0.31 0.36 0.41
## p[47] 0.88 0.00 0.05 0.76 0.84 0.89 0.92
## p[48] 0.50 0.00 0.08 0.34 0.44 0.50 0.55
## lp__ -533.70 0.29 5.52 -544.94 -537.09 -533.26 -529.74
## 97.5% n_eff Rhat
## alpha_pond[1] 3.98 1000 1.00
## alpha_pond[2] 5.56 441 1.00
## alpha_pond[3] 2.21 1000 1.00
## alpha_pond[4] 5.20 1000 1.00
## alpha_pond[5] 4.07 1000 1.00
## alpha_pond[6] 3.97 1000 1.00
## alpha_pond[7] 5.72 1000 1.00
## alpha_pond[8] 4.12 1000 1.00
## alpha_pond[9] 0.97 1000 1.00
## alpha_pond[10] 4.11 1000 1.00
## alpha_pond[11] 2.31 1000 1.00
## alpha_pond[12] 1.85 1000 1.00
## alpha_pond[13] 2.36 1000 1.00
## alpha_pond[14] 1.33 1000 1.00
## alpha_pond[15] 4.04 1000 1.01
## alpha_pond[16] 4.23 1000 1.00
## alpha_pond[17] 4.77 664 1.00
## alpha_pond[18] 3.96 1000 1.00
## alpha_pond[19] 3.25 1000 1.00
## alpha_pond[20] 5.94 1000 1.00
## alpha_pond[21] 3.79 1000 1.00
## alpha_pond[22] 3.82 1000 1.01
## alpha_pond[23] 3.73 1000 1.00
## alpha_pond[24] 2.83 1000 1.00
## alpha_pond[25] -0.17 1000 1.00
## alpha_pond[26] 0.94 1000 1.00
## alpha_pond[27] -0.56 1000 1.00
## alpha_pond[28] 0.31 1000 1.00
## alpha_pond[29] 0.97 1000 1.00
## alpha_pond[30] 2.48 1000 1.00
## alpha_pond[31] 0.16 1000 1.00
## alpha_pond[32] 0.40 1000 1.00
## alpha_pond[33] 4.67 1000 1.00
## alpha_pond[34] 4.10 1000 1.00
## alpha_pond[35] 4.01 1000 1.00
## alpha_pond[36] 3.05 1000 1.00
## alpha_pond[37] 3.15 1000 1.00
## alpha_pond[38] 6.09 1000 1.00
## alpha_pond[39] 4.15 1000 1.00
## alpha_pond[40] 3.55 1000 1.00
## alpha_pond[41] -0.95 820 1.00
## alpha_pond[42] 0.10 1000 1.00
## alpha_pond[43] 0.23 1000 1.00
## alpha_pond[44] 0.38 1000 1.00
## alpha_pond[45] 1.25 1000 1.00
## alpha_pond[46] 0.08 1000 1.00
## alpha_pond[47] 3.22 1000 1.00
## alpha_pond[48] 0.64 1000 1.00
## alpha 1.81 1000 1.00
## sigma 2.11 403 1.01
## p[1] 0.98 1000 1.00
## p[2] 1.00 1000 1.00
## p[3] 0.90 1000 1.00
## p[4] 0.99 1000 1.00
## p[5] 0.98 1000 1.00
## p[6] 0.98 1000 1.00
## p[7] 1.00 1000 1.00
## p[8] 0.98 1000 1.00
## p[9] 0.72 1000 1.00
## p[10] 0.98 1000 1.00
## p[11] 0.91 1000 1.00
## p[12] 0.86 1000 1.00
## p[13] 0.91 1000 1.00
## p[14] 0.79 1000 1.00
## p[15] 0.98 1000 1.00
## p[16] 0.99 1000 1.00
## p[17] 0.99 1000 1.00
## p[18] 0.98 1000 1.00
## p[19] 0.96 1000 1.00
## p[20] 1.00 1000 1.00
## p[21] 0.98 1000 1.00
## p[22] 0.98 1000 1.00
## p[23] 0.98 1000 1.00
## p[24] 0.94 1000 1.00
## p[25] 0.46 1000 1.00
## p[26] 0.72 1000 1.00
## p[27] 0.36 1000 1.00
## p[28] 0.58 1000 1.00
## p[29] 0.72 1000 1.00
## p[30] 0.92 1000 1.00
## p[31] 0.54 1000 1.00
## p[32] 0.60 1000 1.00
## p[33] 0.99 1000 1.00
## p[34] 0.98 1000 1.00
## p[35] 0.98 1000 1.00
## p[36] 0.95 1000 1.00
## p[37] 0.96 1000 1.00
## p[38] 1.00 1000 1.00
## p[39] 0.98 1000 1.00
## p[40] 0.97 1000 1.00
## p[41] 0.28 1000 1.00
## p[42] 0.53 1000 1.00
## p[43] 0.56 1000 1.00
## p[44] 0.59 1000 1.00
## p[45] 0.78 1000 1.00
## p[46] 0.52 1000 1.00
## p[47] 0.96 1000 1.00
## p[48] 0.65 1000 1.00
## lp__ -524.08 355 1.00
##
## Samples were drawn using NUTS(diag_e) at Sun Jul 29 22:37:40 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).traceplot(fit, pars=c("sigma"))
plot(fit, pars=c("alpha", "sigma"))## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)
ms <- extract(fit, permuted=TRUE)#library(shinystan)
#sso <- launch_shinystan(fit)