CARBayesで2次元のCAR：Taglibro de H：SSブログ

	ブログをはじめるログイン

CARBayesで2次元のCAR　[統計] [編集]

前回のデータを今度はRのCARBayesパッケージで解析してみる。以前ためしたときの記事→R: CARBayesをためす。

以前は、分布別に関数がわかれていたが、CARBayes 3.0ではまとめられていて、familyで指定するようになっている。

library(CARBayes)
qg <- read.csv("Qglauca.csv")

n.site <- nrow(qg)
adj <- matrix(0, nrow = n.site, ncol = n.site)
## corner
adj[1, 2] <- 1
adj[1, 11] <- 1
adj[10, 9] <- 1
adj[10, 20] <- 1
adj[191, 181] <- 1
adj[191, 192] <- 1
adj[200, 190] <- 1
adj[200, 199] <- 1
## upper side
for (i in seq(11, 181, 10)) {
  adj[i, i - 10] <- 1
  adj[i, i + 1] <- 1
  adj[i, i + 10] <- 1
}
## lower side
for (i in seq(20, 190, 10)) {
  adj[i, i - 10] <- 1
  adj[i, i - 1] <- 1
  adj[i, i + 10] <- 1
}
## left side
for (i in 2:9) {
  adj[i, i - 1] <- 1
  adj[i, i + 10] <- 1
  adj[i, i + 1] <- 1
}
## right side
for (i in 192:199) {
  adj[i, i - 1] <- 1
  adj[i, i - 10] <- 1
  adj[i, i + 1] <- 1
}
## inside
for (i in seq(11, 181, 10)) {
  for (j in 1:8) {
    adj[i + j, i + j - 1] <- 1
    adj[i + j, i + j - 10] <- 1
    adj[i + j, i + j + 1] <- 1
    adj[i + j, i + j + 10] <- 1
  }
}

## Intrinsic Conditional Autoregressive model
fit <- iarCAR.re(N ~ 1, W = adj, family = "poisson", data = qg,
                 burnin = 1000, n.sample = 41000, thin = 20)

隣のメッシュを指定するところはもっと効率的にできそうな気もするが。

結果

軌跡を表示

> plot(fit$samples[["beta"]])

要約量の表示

> print(fit)
#################
#### Model fitted
#################
Likelihood model - Poisson (log link function) 
Random effects model - Intrinsic CAR
Regression equation - N ~ 1

############
#### Results
############
Posterior quantiles and DIC

             Median    2.5%   97.5% n.sample % accept
(Intercept) -0.5132 -0.7926 -0.2722     2000       60
tau2         1.2885  0.6706  2.3089     2000      100

DIC =  489.3751       p.d =  53.7228

元のデータと95%信用区間を重ねてプロットしてみる。

library(ggplot2)
ggplot(data.frame(x = 1:n.site,
                  y = qg$N,
                  ymean = fit$fitted.values[, "Mean"],
                  ymax = fit$fitted.values[, "97.5%"],
                  ymin = fit$fitted.values[, "2.5%"])) +
    geom_point(aes(x = x, y = y)) +
    geom_pointrange(aes(x = x, y = ymean, ymax = ymax, ymin = ymin),
                  colour = "red", alpha = 0.5, shape = 3) +
    xlab("Index") + ylab("N") + ylim(0, 9)

結果は、OpenBUGSのcar.normalにちかいようだ。