R: glmmADMBによるMCMC [統計]
きのうのつづきで、glmmADMBによるベイズ推定をやってみる。
データはきのうとおなじ。glmmadmb()の引数でmcmc = TRUEとする。MCMCの設定はmcmc.opts = mcmcControl(...)にて。
## MCMC fit2 <- glmmadmb(y ~ 1, random = ~ 1|b, family = "nbinom", mcmc = TRUE, mcmc.opts = mcmcControl(mcmc = 2000), data = data)
結果
> summary(fit2$mcmc) Iterations = 1:1000 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 (Intercept) 2.06260 0.1950 0.006165 0.01769 tmpL 0.55821 0.1728 0.005463 0.01272 alpha 4.54474 1.1785 0.037266 0.08812 u.01 -0.73272 0.3211 0.010155 0.02574 u.02 0.01172 0.3008 0.009512 0.02638 u.03 1.14808 0.4117 0.013019 0.03428 u.04 -0.76447 0.3233 0.010222 0.02588 u.05 0.01172 0.3008 0.009512 0.02638 u.06 1.62942 0.4966 0.015705 0.04105 u.07 0.54267 0.3324 0.010512 0.02911 u.08 -1.27532 0.3644 0.011525 0.02874 u.09 -0.52429 0.3092 0.009777 0.02623 u.10 -0.38835 0.3036 0.009602 0.02603 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% (Intercept) 1.69472 1.9251 2.064920 2.1961 2.4271 tmpL 0.31492 0.4337 0.529334 0.6461 0.9999 alpha 2.61871 3.7167 4.413139 5.2493 7.3281 u.01 -1.33372 -0.9592 -0.742789 -0.4888 -0.1149 u.02 -0.54210 -0.2060 0.002575 0.2416 0.5681 u.03 0.39176 0.8536 1.129203 1.4311 1.9364 u.04 -1.36703 -0.9959 -0.781797 -0.5207 -0.1421 u.05 -0.54210 -0.2060 0.002575 0.2416 0.5681 u.06 0.67098 1.2792 1.607220 1.9780 2.5829 u.07 -0.08132 0.2990 0.537961 0.7605 1.1594 u.08 -1.96390 -1.5292 -1.298649 -1.0068 -0.5866 u.09 -1.09830 -0.7524 -0.532121 -0.2875 0.0624 u.10 -0.94048 -0.6162 -0.399869 -0.1627 0.1765
Gewekeの収束診断をしてみる。
> geweke.diag(fit2$mcmc) Fraction in 1st window = 0.1 Fraction in 2nd window = 0.5 (Intercept) tmpL alpha u.01 u.02 u.03 -0.17838 0.09925 -0.52221 0.48390 0.39194 0.24346 u.04 u.05 u.06 u.07 u.08 u.09 0.48440 0.39194 0.20112 0.31917 0.47068 0.42502 u.10 0.42238
切片の軌跡と密度をプロットしてみる。
> plot(fit2$mcmc[, "(Intercept)"])
こんばんは、魚類の研究をしている大学院生です。
いつもブログで勉強させて頂いています。
glmmADMBパッケージでMCMCを施行するとR2admb::read_psv(file_name) でエラー: no PSV file foundというエラーメッセージが出て結果が出力されません。Rhelp等も当たってみましたが、解決しませんでした。
コードはmodel <- glmmadmb(YY~A+(1|Location), data=dat, zeroInflation=TRUE, family="binomial", mcmc = TRUE, mcmc.opts = mcmcControl(mcmc = 50000))です。
YYが目的変数、Aが説明変数、Locationがランダム効果です。
解決策等ご存知でしたら、ご教授頂けますと幸いです。
よろしくお願い申し上げます。
by fisher (2017-06-07 01:35)
なるほど、現在のバージョンではエラーになりますね。
検索してみると同様の報告がありました。
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2016q1/024458.html
glmmadmbを入れ替えてみるとよいとのことでしたが、ADMBのサイトからはglmmadmbがきえているようです。
http://www.admb-project.org/buildbot/glmmadmb/
ということで、解決策はいまのところわかりません。お役に立てずすみません。
by hiroki (2017-06-10 05:45)
長らく海外調査に行っており、お礼のお返事が遅くなってしまいました。
大変申し訳ございませんでした。
丁寧な解説をありがとうございました。
なるほど、バージョンによって変わってしまうのですね…。
by fisher (2017-07-09 20:24)