R: MCMCpackによる、ランダム効果つきポアソン回帰 [統計]
MCMCpackをあらためてながめていたら、実はいろいろできることに気がついたので、ランダム効果のあるポアソン回帰など。
久保みどり本 10.5節のデータでやってみます。MCMCpackでは、MCMChpoisson関数を使用します。
> library(MCMCpack) > > data <- read.csv(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/hbm/nested/d1.csv")) > > fit <- MCMChpoisson(fixed = y ~ f, random = ~1, group = "pot", data = data, + burnin = 5000, mcmc = 1000, thin = 1, verbose = 1, + r = 1, R = matrix(0.001), nu = 0.001, delta = 0.001) Running the Gibbs sampler. It may be long, keep cool :) **********:10.0%, mean accept. rate=0.466 **********:20.0%, mean accept. rate=0.449 **********:30.0%, mean accept. rate=0.459 **********:40.0%, mean accept. rate=0.441 **********:50.0%, mean accept. rate=0.441 **********:60.0%, mean accept. rate=0.464 **********:70.0%, mean accept. rate=0.448 **********:80.0%, mean accept. rate=0.441 **********:90.0%, mean accept. rate=0.441 **********:100.0%, mean accept. rate=0.434 > summary(fit$mcmc) Iterations = 5001:6000 Thinning interval = 1 Number of chains = 1 Sample size per chain = 1000 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE beta.(Intercept) 1.3700 0.4546 0.014374 0.01556 beta.fT -0.8474 0.6311 0.019956 0.02545 b.(Intercept).A 0.3681 0.5190 0.016413 0.01861 b.(Intercept).B -0.3282 0.5149 0.016283 0.01935 b.(Intercept).C -0.5722 0.5318 0.016818 0.02222 b.(Intercept).D 0.6869 0.5086 0.016084 0.01743 b.(Intercept).E -0.1717 0.5125 0.016205 0.02060 b.(Intercept).F -0.8380 0.5577 0.017635 0.02971 b.(Intercept).G 0.6006 0.5349 0.016914 0.02593 b.(Intercept).H -0.8379 0.5682 0.017968 0.03089 b.(Intercept).I 1.2574 0.5407 0.017099 0.03229 b.(Intercept).J -0.1821 0.5414 0.017120 0.03030 VCV.(Intercept).(Intercept) 0.8333 0.6889 0.021784 0.04405 sigma2 1.0482 0.2343 0.007408 0.02293 Deviance 368.6206 13.1145 0.414715 0.92102 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% beta.(Intercept) 0.4753 1.11171 1.3678 1.652583 2.2754 beta.fT -2.0790 -1.22934 -0.8428 -0.465608 0.3842 b.(Intercept).A -0.5691 0.04847 0.3528 0.681482 1.4369 b.(Intercept).B -1.3785 -0.63960 -0.3133 -0.003193 0.6271 b.(Intercept).C -1.6971 -0.91092 -0.5602 -0.246593 0.4784 b.(Intercept).D -0.2789 0.35475 0.6707 0.999336 1.7293 b.(Intercept).E -1.1833 -0.50586 -0.1636 0.148726 0.8077 b.(Intercept).F -1.9300 -1.20657 -0.8252 -0.436459 0.1785 b.(Intercept).G -0.3752 0.25301 0.5624 0.932266 1.6750 b.(Intercept).H -2.1003 -1.15601 -0.7928 -0.487770 0.1638 b.(Intercept).I 0.2658 0.90008 1.2233 1.576884 2.4288 b.(Intercept).J -1.2328 -0.50982 -0.1802 0.151011 0.8465 VCV.(Intercept).(Intercept) 0.1878 0.40011 0.6400 1.012997 2.7570 sigma2 0.6830 0.88140 1.0157 1.186421 1.5983 Deviance 343.9016 359.29448 368.2692 376.830900 395.8377
ランダム効果の事前分布は逆Wishart分布で、rがshapeパラメーター、Rがscale行列ということです。マニュアルによれば、無情報事前分布ならr=ランダム効果の数、だそうで、ここでは r=1, r = 1, R = matrix(0.001) としました。この事前分布は以下のような確率分布になります。
> x <- seq(0.01, 1, length = 100) > p <- sapply(x, function(x) diwish(W = matrix(x), v = 1, S = 0.001)) > plot(p ~ x, type = "l", las = 1)
コメント 0