R: hSDMパッケージをためす [統計]
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)
コメント 0