Jones, J. E.; Kroll, A. J.; Giovanini, J.; Duke, S. D.; Ellis, T. M.; Betts, M. G. Avian Species Richness in Relation to Intensive Forest Management Practices in Early Seral Tree Plantations. PLoS ONE 2012, 7, e43290.
지난 번에 발표한 논문 실행 코드를 올립니다. 제가 작성한 부분은 원저자가 만들어 놓은 BUGS 코드를 R에서 돌려서 결과 보는 부분만 입니다. 실제로 돌리려면 PLOS 가서 BUGS 코드 다운 받고, 텍스트 파일로 만들어야 함.
Jones, J. E.; Kroll, A. J.; Giovanini, J.; Duke, S. D.; Ellis, T. M.; Betts, M. G. Avian Species Richness in Relation to Intensive Forest Management Practices in Early Seral Tree Plantations. PLoS ONE 2012, 7, e43290.
Appendix S4: WinBUGS code for hierarchical community model and R code for calculating Average Predictive Comparisons, Oregon Coast Range, USA, 2008-2009.
데이터는 올라와 있지 않아 가상 데이터를 생성해서 돌아가도록 했습니다. 지난 번 토의 내용에 연결된 내용 조금 적겠습니다.
무엇을 말 할 수 있는가?
-> mu.psi 는 평균 관측확률, 새를 얼마나 수월하게 관측할 수 있었는지 정보를 줌. Z는 사이트에 새가 존재했는지 아니었는지 (P or A). 종 별, 길드 별로 볼 수 있다.
관측 확률은 어떻게 정보를 주는가?
-> b0, b1, ... , b4 가 관측확률을 회귀하는 식 계수. 이 계수에 사전에 생각하고 있던 값을 넣을 수 있다. 자료가 부족할 때는 사전에 넣은 계수에 좌지우지되지만 자료가 많아짐에 따라 자료에 적합된 계수가 도출됨 (사후 분포를 얻음).
# R code for hierarchical community model and average predictive comparisons of species richness.
# Ported to the JAGS-rjags environment by Alan B. Seo, 21 Feb 2013
require(rjags) # This loads in the library that uses the 'JAGS' Bayesian analysis package
setwd("~/Dropbox/BayesianBird/")
model.loc <- "AvianSpeciesRichness_bugs.R" # The model specification file
data.loc <- "AvianSpeciesRichness.csv" # Which we don't have
# --- POSSIBLE VALUES TO CHANGE ---
num.chains <- 1 # Number of MCMC chains to run
burnin <- 10000 # The number of iterations to perform at the beginning (the 'burn-in')
### ONLY IN ADAPTATION PHASE BUGS ADJUST PROPOSAL FUNCTIONS!!! After burn-in, nothing changes )chain properties)
num.iters <- 10000 # The number of iterations to perform after the 'burn-in'
thin <- 1 # A thinning parameter (should usually be 1 unless you're running a lot of samples)
# ---------------------------------
# --- TO COMPLETE ---
# Change the next line to a character vector of the name of the parameters that you would like to save (as they appear in they appear in the JAGs model specification)
vars.to.save <- c( "mu.psi", "a0", "b0", "Z")
# -------------------
# Create a virtual data of bird observation.
data = array(rbinom(n=50, size=1, prob=0.05), dim=c(1,50,1)) # One species, one site and 50 observation events
inits <- list("mu.a0" = 0.1, "mu.a1" = 0.1, "mu.a2" = 0.1, "mu.a3" = 0.1, "mu.a4" = 0.1, "mu.a5" = 0.1, "mu.b0" = 0.1, "mu.b1" = 0.1, "mu.b2" = 0.1, "mu.b3" = 0.1, "mu.b4" = 0.1, "b0"=0) # Initial mu values
avian.list <- list(
n = 1, J=1, K=50, leafgleaners = rep(1, times=50), X = data,
ELEVM = c(100), CONIFER=c(0.5), BROADsansDEC=c(0.3), DECsansHWD=c(0.2), HWD=c(0.6)
) # Number of timeperiods to pass to JAGS as the parameter 'total.T'
# Create and initialise the model
statespace.model <- jags.model(file=model.loc, n.chains = num.chains, n.adapt = burnin, data=avian.list, quiet=F, inits= inits)
# Get samples from the posterior
model.output <- coda.samples(statespace.model, vars.to.save, n.iter = num.iters, thin = 1)
# Posterior description
summary(model.output)
plot(model.output)