Dynamic community occupancy model in R and JAGS example

community_occ.R

# Multi-species dynamic occupancy model with R and JAGS
# Written by Max Joseph
# maxwell.b.joseph@colorado.edu
# see http://www.colorado.edu/eeb/gradstudents/joseph/community_occ.html
# for details

# convenience functions
logit <- function(x) {
  log(x/(1 - x))
}

antilogit <- function(x) {
  exp(x)/(1 + exp(x))
}

# initialize parameters
nsite <- 150
nspec <- 6
nyear <- 4
nrep <- 3

# community level hyperparameters
p_beta = 0.7
mubeta <- logit(p_beta)
sdbeta <- 2

p_rho <- 0.8
murho <- logit(p_rho)
sdrho <- 1

# species specific random effects
set.seed(1)  # for reproducability
beta <- rnorm(nspec, mubeta, sdbeta)
set.seed(1008)
rho <- rnorm(nspec, murho, sdrho)

# initial occupancy states
set.seed(237)
rho0 <- runif(nspec, 0, 1)
z0 <- array(dim = c(nsite, nspec))
for (i in 1:nspec) {
  z0[, i] <- rbinom(nsite, 1, rho0[i])
}

# subsequent occupancy
z <- array(dim = c(nsite, nspec, nyear))
lpsi <- array(dim = c(nsite, nspec, nyear))
psi <- array(dim = c(nsite, nspec, nyear))
for (j in 1:nsite) {
  for (i in 1:nspec) {
    for (t in 1:nyear) {
      if (t == 1) {
        lpsi[j, i, t] <- beta[i] + rho[i] * z0[j, i]
        psi[j, i, t] <- antilogit(lpsi[j, i, t])
        z[j, i, t] <- rbinom(1, 1, psi[j, i, t])
      } else {
        lpsi[j, i, t] <- beta[i] + rho[i] * z[j, i, t - 1]
        psi[j, i, t] <- antilogit(lpsi[j, i, t])
        z[j, i, t] <- rbinom(1, 1, psi[j, i, t])
      }
    }
  }
}

# detection probabilities
p_p <- 0.7
mup <- logit(p_p)
sdp <- 1.5
set.seed(222)
lp <- rnorm(nspec, mup, sdp)
p <- antilogit(lp)

# observations
x <- array(dim = c(nsite, nspec, nyear, nrep))
for (j in 1:nsite) {
  for (i in 1:nspec) {
    for (t in 1:nyear) {
      for (k in 1:nrep) {
        x[j, i, t, k] <- rbinom(1, 1, p[i] * z[j, i, t])
      }
    }
  }
}

# model specification
cat("
model{
    #### priors
    # beta hyperparameters
    p_beta ~ dbeta(1, 1)
    mubeta <- log(p_beta / (1 - p_beta))
    sigmabeta ~ dunif(0, 10)
    taubeta <- (1 / (sigmabeta * sigmabeta))
    
    # rho hyperparameters
    p_rho ~ dbeta(1, 1)
    murho <- log(p_rho / (1 - p_rho))
    sigmarho~dunif(0, 10)
    taurho<-1 / (sigmarho * sigmarho)
    
    # p hyperparameters
    p_p ~ dbeta(1, 1)
    mup <- log(p_p / (1 - p_p))
    sigmap ~ dunif(0,10)
    taup <- (1 / (sigmap * sigmap))
    
    #### occupancy model
    # species specific random effects
    for (i in 1:(nspec)) {
      rho0[i] ~ dbeta(1, 1)
      beta[i] ~ dnorm(mubeta, taubeta)
      rho[i] ~ dnorm(murho, taurho)
    }
    
    # occupancy states
    for (j in 1:nsite) {
      for (i in 1:nspec) {
        z0[j, i] ~ dbern(rho0[i])
        logit(psi[j, i, 1]) <- beta[i] + rho[i] * z0[j, i] 
        z[j, i, 1] ~ dbern(psi[j, i, 1]) 
        for (t in 2:nyear) {
          logit(psi[j, i, t]) <- beta[i] + rho[i] * z[j, i, t-1]
          z[j, i, t] ~ dbern(psi[j, i, t])
        }
      }
    }
    
    #### detection model
    for(i in 1:nspec){ 
      lp[i] ~ dnorm(mup, taup)
      p[i] <- (exp(lp[i])) / (1 + exp(lp[i]))
    }
    
    #### observation model
    for (j in 1:nsite){
      for (i in 1:nspec){
        for (t in 1:nyear){
          mu[j, i, t] <- z[j, i, t] * p[i] 
          for (k in 1:nrep){
            x[j, i, t, k] ~ dbern(mu[j, i, t])
          }
        }
      }
    }
  }
  ", fill=TRUE, file="com_occ.txt")

# bundle data
data <- list(x = x, nrep = nrep, nsite = nsite, nspec = nspec, nyear = nyear)

# initial values
zinit <- array(dim = c(nsite, nspec, nyear))
for (j in 1:nsite) {
  for (i in 1:nspec) {
    for (t in 1:nyear) {
      zinit[j, i, t] <- max(x[j, i, t, ])
    }
  }
}

inits <- function() {
  list(p_beta = runif(1, 0, 1), p_rho = runif(1, 0, 1), 
       sigmarho = runif(1, 0, 1), sigmap = runif(1, 0, 10), 
       sigmabeta = runif(1, 0, 10), z = zinit)
}

# parameters to monitor
params <- c("lp", "beta", "rho")

require(rjags)
# build model
ocmod <- jags.model(file = "com_occ.txt", inits = inits, data = data, n.chains = 3)

# specify MCMC settings and start sampling
nburn <- 2000
update(ocmod, n.iter = nburn)
out <- coda.samples(ocmod, n.iter = 7000, variable.names = params)
summary(out)

# check convergence
plot(out)

# compare parameter estimates to true values
require(mcmcplots)
caterplot(out, "beta", style = "plain")
caterpoints(beta)
caterplot(out, "lp", style = "plain")
caterpoints(lp)
caterplot(out, "rho", style = "plain")
caterpoints(rho)