########### DESCRIPTION # Bayesian occupancy analysis of Beetle and Burn Data # published in Tingley et al., 2023, PLOS ONE ########### ## Load Libraries library(R2jags) library(rjags) library(foreign) library(dclone) library(MCMCvis) library(runjags) ## Load data folder<-getwd() load("Occ_data_beetles.Rdata") load("Occ_data_birds.Rdata") attach(export.birds) attach(export.beetles) ## Limit Beetle data to matching sites in woodpecker data sites <- sites[sites$Trap_code %in% events$Station_visit_ID, ] ################ Beetle activity metrics # creating empty frames bar_s <- array(dim = c(dim(sites)[1], 6)) bar_n <- array(dim = c(dim(sites)[1], 6)) spec.tree <- array(dim = c(dim(sites)[1], 6)) # filling in frames site by site for(i in 1:dim(sites)[1]) { data.j <- beetles[beetles$Trap_code == sites$Trap_code[i], ] bar_s[i, 1:dim(data.j)[1]] <- data.j$Beetle_activity_rating_south bar_n[i, 1:dim(data.j)[1]] <- data.j$Beetle_activity_rating_north spec.tree[i, 1:dim(data.j)[1]] <- data.j$Species } # creating a composite activity index for each tree at each site activity <- (bar_s + bar_n) # index of the proportion of pines in measured samples at each site prop.pine <- apply(array(data = spec.tree %in% c("Jeffrey Pine", "Lodgepole Pine", "Ponderose Pine", "Sugar pine", "Yellow Pine")*1, dim = dim(spec.tree)), 1, sum) / (6 - apply(is.na(bar_s), 1, sum)) ######### Format woodpeckers into site x year array X.bbwo <- X.bbwo[, c(1:3, 6:8)] years <- unique(events$Year) bbwo.3d <- array(dim = c(length(sites$Trap_line_code), dim(X.bbwo)[2], length(unique(events$Year)))) # day of year array jday.2d <- array(dim = c(length(sites$Trap_line_code),length(unique(events$Year)))) # snag density array snag.2d <- array(dim = c(length(sites$Trap_line_code),length(unique(events$Year)))) # years since fire array age.2d <- array(dim = c(length(sites$Trap_line_code),length(unique(events$Year)))) # fill in arrays for(t in 1:dim(bbwo.3d)[3]) { x.tmp <- X.bbwo[events$Year == years[t], ] jday.tmp <- jday[events$Year == years[t]] snag.tmp <- snag[events$Year == years[t]] age.tmp <- fire.age[events$Year == years[t]] ev.tmp <- events[events$Year == years[t], ] bbwo.3d[,,t] <- x.tmp[match(sites$Trap_line_code, ev.tmp$Station_code), ] jday.2d[, t] <- jday.tmp[match(sites$Trap_line_code, ev.tmp$Station_code)] snag.2d[, t] <- snag.tmp[match(sites$Trap_line_code, ev.tmp$Station_code)] age.2d[, t] <- age.tmp[match(sites$Trap_line_code, ev.tmp$Station_code)] } age.2d <- matrix(data = -9:0, nrow = dim(age.2d)[1], ncol = dim(age.2d)[2], byrow = T) + age.2d[, 10] age.2d[age.2d < 0] <- 0 # can't have years since fire < 0 # fill in missing predictor values (these are surveys that did not occur, so they take the mean) jday.2d[is.na(jday.2d)] <- mean(jday) snag.2d[is.na(snag.2d)] <- mean(snag) ##### Restrict woodpecker data to just sites that were visited in 2018 which.bbwo <- match(sites$Trap_code, events$Station_visit_ID) ############################################### JAGS data ##################### jags.data <- list(y = bbwo.3d, jday = array(scale(as.vector(jday.2d)), dim = dim(jday.2d)), ef = ef[c(1:3, 6:8)], type = itype[c(1:3, 6:8)], snag = array(scale(as.vector(snag.2d)), dim = dim(snag.2d)), age = array(scale(as.vector(age.2d)), dim = dim(age.2d)), elev = scale(elev)[which.bbwo, 1], lat = scale(lat)[which.bbwo, 1], fire = as.numeric(factor(sites$Fire)), n.fires = length(unique(sites$Fire)), n.points = length(which.bbwo), n.surveys = dim(bbwo.3d)[2], n.trees = 6 - apply(is.na(bar_s), 1, sum), n.years = length(years), activ = apply(activity, 1, sum, na.rm = T), prop.pine = prop.pine) # build a naive z-matrix as initial values to keep JAGS happy z.naive <- apply(jags.data$y, MARGIN = c(1, 3), max, na.rm=T) z.naive[(z.naive) == "-Inf"] <- 1 inits <- function() {list(z = z.naive)} #### This is the main model code, written in the JAGS modeling language jags.mod <- function() { #PRIORS #priors for intercepts psi0 ~ dunif(0, 1) p0 ~ dunif(0, 1) a.0 <- log(p0 / (1 - p0)) mu.b0 <- log(psi0 / (1 - psi0)) tau.b0 ~ dgamma(0.1, 0.1) #priors for logit-linear model coefficients a.type ~ dnorm(0,0.1) b.elev ~ dnorm(0,0.1) b.lat ~ dnorm(0,0.1) b.snag ~ dnorm(0,0.1) b.activ ~ dnorm(0, 0.1) b.activage ~ dnorm(0, 0.1) phi ~ dnorm(0, 0.1) g.0 ~ dnorm(0,0.1) g.yr ~ dnorm(0,0.1) g.pine ~ dnorm(0,0.1) g.yrpine ~ dnorm(0,0.1) # Random effect for each fire for(l in 1:n.fires) { b.0[l] ~ dnorm(mu.b0, tau.b0) } # Main model for (j in 1:n.points) { logit(intensity[j, 1]) <- g.0 + g.yr * age[j, 1] + g.pine * prop.pine[j] + g.yrpine * age[j, 1] * prop.pine[j] activ[j] ~ dbinom(intensity[j, 10], n.trees[j] * 8) logit(psi[j, 1]) <- b.0[fire[j]] + elev[j] * b.elev + lat[j] * b.lat + snag[j, 1] * b.snag + intensity[j, 1] * b.activ + intensity[j, 1] * age[j, 1] * b.activage z[j, 1] ~ dbern(psi[j, 1]) for(k in 1:n.surveys) { logit(p[j, k, 1]) <- a.0 + a.type * type[k] mu.p[j, k, 1] <- z[j, 1] * p[j, k, 1] y[j, k, 1] ~ dbern(mu.p[j, k, 1]) } for(t in 2:n.years) { logit(intensity[j, t]) <- g.0 + g.yr * age[j, t] + g.pine * prop.pine[j] + g.yrpine * age[j, t] * prop.pine[j] logit(psi[j, t]) <- b.0[fire[j]] + elev[j] * b.elev + lat[j] * b.lat + snag[j, t] * b.snag + intensity[j, t] * b.activ + intensity[j, t] * age[j, t] * b.activage + phi * z[j, t - 1] z[j, t] ~ dbern(psi[j, t]) for(k in 1:n.surveys) { logit(p[j, k, t]) <- a.0 + a.type * type[k] mu.p[j, k, t] <- z[j, t] * p[j, k, t] y[j, k, t] ~ dbern(mu.p[j, k, t]) } } } } ########## RUN JAGS ############### params = c("a.0", "a.ef", "a.jday", "a.type", "mu.b0", "b.elev", "b.lat", "b.snag", "b.activ", "b.activage", "phi", "g.0", "g.yr", "g.pine", "g.yrpine") setwd(folder) nc <- 3 n.adapt <- 10000 n.burn <- 50000 n.iter <- 50000 thin <- 50 ## Parallelize ## start.time<-Sys.time() cl <- makePSOCKcluster(nc) # Make your socket cluster and name it 'cl'. nc is the number of cores you want your cluster to have--generally the number of chains you want to run, hence the name. tmp <- clusterEvalQ(cl, library(dclone)) # Check that `dclone` is loaded on each of cl's workers. dclone includes JAGS functionality parLoadModule(cl, "glm") # load the JAGS module 'glm' on each worker parListModules(cl) # make sure previous line worked. beetle_pars <- jags.parfit(cl, jags.data, params = params, model = jags.mod, inits = inits, n.chains = nc, n.adapt = n.adapt, n.update = n.burn, thin = thin, n.iter = n.iter) stopCluster(cl) # close out the cluster. end.time=Sys.time() elapsed.time <- difftime(end.time, start.time, units='mins') #2 minutes for 5K iterations elapsed.time ## save setwd(folder) save(beetle_pars, file = paste('pars_multiyr_', Sys.Date(), '.Rdata', sep='')) ###### view results summary(beetle_pars) MCMCplot(beetle_pars, ref_ovl = T) MCMCsummary(beetle_pars) MCMCtrace(beetle_pars)