R: glmmAKによる累積ロジットモデルのベイズ推定

glmmAKで累積ロジットモデルのベイズ推定をしてみる。累積ロジットモデルのベイズ推定(2)で使ったデータを使用。時間によるランダム効果は無視した。

# inverse logit
invlogit <- function(x) exp(x)/(1+exp(x))

set.seed(12345)
n.rank <- 5
plot <- 1:20
n.plot <- length(plot)
time <- 1:5
n.time <- length(time)
n <- n.plot * n.time
data <- expand.grid(plot = plot, time = time)
slope <- rgamma(n.plot, shape = 15^2/5^2, rate = 15/5^2)
data$slope <- rep(slope, n.time)
ranef.plot <- rnorm(n.plot, 0, 1)
ranef.time <- rnorm(n.time, 0, 1)
logit.c1 <- 7.5 - 0.5 * slope +
            ranef.plot[data$plot] + ranef.time[data$time]
data$c1 <- rbinom(n, n.rank - 1, invlogit(logit.c1)) + 1

library(glmmAK)

wd <- getwd()
fit <- cumlogitRE(y = data$c1 - 1,
                  x = data.frame(slope = data$slope),
                  cluster = data$plot,
                  intcpt.random = TRUE,
                  hierar.center = FALSE,
                  drandom = "normal",
                  C = max(data$c1 - 1),
                  logit.order = "decreasing",
                  prior.fixed = list(mean = 0, var = 10000),
                  prior.random = list(Ddistrib = "gamma",
                                      Dshape = 0.001,
                                      DinvScale = 0.001),
                  nsimul = list(niter = 2000, nthin = 20,
                                  nburn = 1000, nwrite = 20),
                  dir = wd)
post <- glmmAK.files2coda(dir = wd,
                          drandom = "normal", 
                          skip = 0)

summary(post$betaF)

結果

Iterations = 1001:2000
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):1 10.8667 1.75427 0.055475       0.067901
(Intercept):2  8.8216 1.57604 0.049839       0.059188
(Intercept):3  7.7310 1.51383 0.047871       0.056044
(Intercept):4  5.5203 1.36882 0.043286       0.047840
slope         -0.5514 0.09632 0.003046       0.003524

2. Quantiles for each variable:

                 2.5%     25%    50%     75%   97.5%
(Intercept):1  7.9021  9.6676 10.652 12.0425 14.5893
(Intercept):2  6.1378  7.7112  8.657  9.8661 11.9907
(Intercept):3  5.2030  6.6658  7.627  8.7007 10.8928
(Intercept):4  3.0617  4.6209  5.363  6.3975  8.3975
slope         -0.7555 -0.6134 -0.543 -0.4817 -0.3815

グラフ

plot(post$betaF)

タグ：累積ロジットモデル glmmAK MCMC R