SSブログ

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)"])

Rplot.png


タグ:glmm MCMC
nice!(1)  コメント(3)  トラックバック(0) 
共通テーマ:日記・雑感

nice! 1

コメント 3

fisher

こんばんは、魚類の研究をしている大学院生です。
いつもブログで勉強させて頂いています。

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) 

hiroki

なるほど、現在のバージョンではエラーになりますね。

検索してみると同様の報告がありました。
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) 

fisher

長らく海外調査に行っており、お礼のお返事が遅くなってしまいました。
大変申し訳ございませんでした。

丁寧な解説をありがとうございました。
なるほど、バージョンによって変わってしまうのですね…。


by fisher (2017-07-09 20:24) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0