#******************************************************************************* # 14.4 Territory occupancy model to estimate survival using BUGS #******************************************************************************* # # Korner-Nievergelt, F., T. Roth, S. von Felten, J. Guelat, B. Almasi and # P. Korner-Nievergelt. 2015. Bayesian Data Analysis in Ecology using Linear # Models with R, BUGS and Stan. Elsevier. # # Last modification: 5 June 2014 #------------------------------------------------------------------------------- library(blmeco) # package to the book #------------------------------------------------------------------------------- # Prepare data #------------------------------------------------------------------------------- # # Example data: Territory occupancy data of nightingales # The data was collected and kindly provided by Valentin Amrhein # # The data contains the territory occupancy data y[i,t,j] # i = territory # t = year (1=2000, 10=2009) # j = visit (eight vistis per year) # if y=1 --> nightingale observed # if y=0 --> no nightingale observed # Load data data(nightingales) y <- nightingales class(y) dim(y) # Prepare some additional data N <- dim(y)[1] # Number of territories T <- dim(y)[2] # Number of study years J <- dim(y)[3] # Number of yearly visits # bundle the data datax <- list(y=y, N=N, T=T, J=J) #------------------------------------------------------------------------------- # Calculate initial values #------------------------------------------------------------------------------- # Calculate x[i,t] --> x[i,t] = 1 if a nightingale was observed at least once during the J visits x <- array(NA, dim = c(N,T)) for(t in 1:T) x[,t] <- as.integer(apply(y[,t,], 1, sum)>0) # Function to calculate initial values inits <- function() { list(omega = runif(1,0,1), phi = runif(1,0,1), r = runif(1,0,1), p = runif(1,0,1), x=x ) } #------------------------------------------------------------------------------- # Run OpenBUGS #------------------------------------------------------------------------------- # Parameters to be monitored parameters <- c("phi", "r", "p", "omega") # Run OpenBUGS library(R2OpenBUGS) bugsdir <- file.path(getwd(), "BUGS") bugsmod <- bugs(datax, inits, parameters, model.file="territoryoccupancy.txt", n.thin=2, n.chains=2, n.burnin=1000, n.iter=5000, debug=TRUE, working.directory = bugsdir) #------------------------------------------------------------------------------- # Predictive model checking #------------------------------------------------------------------------------- nsim <- bugsmod$n.sims nobstot <- integer(nsim) meanobster <- integer(nsim) obsy1 <- integer(nsim) obsy5 <- integer(nsim) for(i in 1:nsim){ x <- array(dim = c(N, T)) z <- array(dim = c(N, T)) yrep <- array(dim = c(N, T, J)) x[,1] <- rbinom(N, 1, bugsmod$sims.matrix[, "omega"][i]) for(j in 1:J) { yrep[,1,j] <- rbinom(N, 1, x[,1]*bugsmod$sims.matrix[, "p"][i]) } for(t in 2:T) { z[,t] <- rbinom(N, 1, x[,t-1]*bugsmod$sims.matrix[, "phi"][i]) x[,t] <- rbinom(N, 1, z[,t] + (1-z[,t])*bugsmod$sims.matrix[, "r"][i]) for(j in 1:J) { yrep[,t,j] <- rbinom(N, 1, x[,t]*bugsmod$sims.matrix[, "p"][i]) } } nobstot[i] <- sum(yrep) nobster <- rep(0,N) for(t in 1:T) { nobster <- nobster + as.integer(apply(yrep[,t,], 1, sum)>0) } meanobster[i] <- mean(nobster) obsy1[i] <- sum(yrep[,,1]) obsy5[i] <- sum(yrep[,,5]) } # Compare summary statistics from simulated data with summary statistics from real data quantile(nobstot, probs = c(0.01, 0.99)) sum(y) quantile(meanobster, probs = c(0.01, 0.99)) nobster <- rep(0,N) for(t in 1:T) { nobster <- nobster + as.integer(apply(y[,t,], 1, sum)>0) } mean(nobster) quantile(obsy1, probs = c(0.01, 0.99)) sum(y[,,1]) quantile(obsy5, probs = c(0.01, 0.99)) sum(y[,,5]) #------------------------------------------------------------------------------- # Inference #------------------------------------------------------------------------------- print(bugsmod$summary[,c("mean", "2.5%", "97.5%", "Rhat")], 3)