JAGSで1次元のIntrinsic CARモデル

beroberoさんの「[Stan] 久保緑本11章のCAR model（空間構造のあるベイズモデル） 訂正版」を参考に、JAGSで1次元のIntrinsic CARモデルを実装してみました。

状態空間モデルと同じ形式になるというのは目からウロコでした。

JAGSのモデル

var
  N.site,
  Y[N.site],
  r[N.site],
  sigma,
  tau;
model {
  for (i in 1:N.site) {
    Y[i] ~ dpois(exp(r[i]));
  }
  for (i in 2:N.site) {
    r[i] ~ dnorm(r[i - 1], tau);
  }
  r[1] ~ dnorm(0, 1.0E-4);
  tau <- 1 / (sigma * sigma);
  sigma ~ dunif(0, 1.0E+4);
}

Rコード

Y <- c( 0,  3,  2,  5,  6, 16,  8, 14, 11, 10,
       17, 19, 14, 19, 19, 18, 15, 13, 13,  9,
       11, 15, 18, 12, 11, 17, 14, 16, 15,  9,
        6, 15, 10, 11, 14,  7, 14, 14, 13, 17,
        8,  7, 10,  4,  5,  5,  7,  4,  3,  1)
N.site <- length(Y)

library(rjags)
model <- jags.model("jags_car_bug.txt",
                    data = list(Y = Y, N.site = N.site),
                    n.chains = 3, n.adapt = 10000)
fit <- coda.samples(model,
                    variable.names = c("r", "sigma"),
                    n.iter = 10000, thin = 10)

Y.hat <- exp(summary(fit[, grep("r", varnames(fit))])$statistics[, "Mean"])
Y.qnt <- exp(summary(fit[, grep("r", varnames(fit))])$quantiles[, c("2.5%", "97.5%")])

library(ggplot2)
df <- data.frame(site = seq_along(Y),
                 Y = Y,
                 Y.hat = Y.hat,
                 Y.975 = Y.qnt[, 2],
                 Y.025 = Y.qnt[, 1])
ggplot(df) +
  geom_point(aes(x = site, y = Y), size = 3) +
  geom_line(aes(x = site, y = Y.hat), size = 1) +
  geom_ribbon(aes(x = site, ymin = Y.025, ymax = Y.975), alpha = 0.3)

結果です。

ただ、2次元になると尤度をいじる必要がでてくるので（Stanなら2重に"~"構文をつかったり、increment_log_probをつかえます）、一筋縄ではいかない模様です。