SSブログ

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)

Rplot001.png


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

nice! 1

コメント 0

コメントを書く

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

Facebook コメント

トラックバック 0