Taglibro de H

R Package for Hierarchical Bayesian species distribution models now available | Adam M. Wilsonという記事で紹介されていたhSDMパッケージをためしてみた。

例題として、久保本11章の空間構造のある階層ベイズモデルをやってみる。

library(hSDM)

load(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/Y.RData"))

n.site <- length(Y)
adj <- c(2, sapply(2:(n.site - 1), function(i) c(i - 1, i + 1)), n.site - 1)
n.adj <- c(1, rep(2, n.site - 2), 1)
data <- data.frame(Y = Y, site = 1:n.site)

post.hSDM <- hSDM.poisson.iCAR(counts = data$Y,
                               suitability = ~1,
                               cells = data$site,
                               n.neighbors = n.adj,
                               neighbors = adj,
                               data = data,
                               burnin = 3000,
                               mcmc = 10000,
                               thin = 5,
                               beta.start = 1,
                               Vrho.start = 1,
                               priorVrho = "Uniform")

summary(post.hSDM$mcmc)
plot(post.hSDM$mcmc)

## 
plot(1:n.site, data$Y, type = "p", las = 1,
     ylim = c(0, 20), xlab = "site", ylab = "Y")
lines(1:n.site, post.hSDM$prob.p.pred)

> summary(post.hSDM$mcmc)

Iterations = 3001:12996
Thinning interval = 5 
Number of chains = 1 
Sample size per chain = 2000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                      Mean      SD  Naive SE Time-series SE
beta.(Intercept)   2.26934 0.04999 0.0011178       0.002623
Vrho               0.05682 0.02606 0.0005827       0.001690
Deviance         247.66520 7.31907 0.1636594       0.322132

2. Quantiles for each variable:

                      2.5%       25%       50%       75%    97.5%
beta.(Intercept)   2.17036   2.23494   2.26986   2.30431   2.3686
Vrho               0.02066   0.03906   0.05149   0.06928   0.1211
Deviance         234.75545 242.51875 247.44579 252.13904 263.3654

> plot(post.hSDM$mcmc)

> plot(1:n.site, data$Y, type = "p", las = 1,
+      ylim = c(0, 20), xlab = "site", ylab = "Y")
> lines(1:n.site, post.hSDM$prob.p.pred)