# Supplement

# Koons, D.N., F. Colchero, K. Hersey, and O. Gimenez. 2015. Disentangling 
# the effects of climate, density dependence, and harvest on an iconic large 
# herbivore’s population dynamics. Ecological Applications 

# Programmer: David N. Koons, 2014
# R version: 3.1.0

###############################################################################
# Annotated code for casting the top-ranked state-space model using R2jags    #
###############################################################################
library(R2jags)

# y: vector of population counts supplied as data below
# z: vector of standardized temperature covariate values supplied as data below
# extract: vector of known number of individuals harvested or translocated
#          supplied as data below. 

sink("gomtop.jags")
cat("
model {

### Priors and constraints ###

# Priors for the beta parameters. We used a strong prior on b0 (growth rate 
# from low density) based on extensive scientific literature. The prior on b1 
# Gompertz DD) was based on theoretical bounds for density regulation to 
# stable attractors (Dennis et al. 2006) with the mode of the probability 
# mass at 0. A vague prior distribution was used for b2, the covariate effect.
taub0 <- pow(0.039,-2)
b0 ~ dnorm(0.277,taub0)
taub1 <- pow(2,-2)
b1 ~ dnorm(0,taub1)I(-2,2)
taub2 <- pow(10,-2)
b2 ~ dnorm(0,taub2)
# Prior for the process variance - inverse gamma (0.001, 0.001); note that 
# tauPro is 1/variance.
tauPro ~ dgamma(0.001,0.001)
sigPro <- 1/pow(tauPro,0.5)
# Variance in observations is subsumed within the Poisson observation 
# model below.

# Initial population size was a known release, so we fixed it to the released
# abundance.
x[1] <- log(y[1]) 
N[1] <- exp(x[1])

### State-Space Model (the Likelihood) ###
# Observation process
for (t in 2:T) {
   y[t] ~ dpois(N[t])   
   }
# State process
for (t in 2:T){
   # solve for the logarithmic extraction term
   e[t-1] <- log(abs(1-extract[t-1]/N[t-1]))
   # Expected value of x[t]
   Ex[t] <- x[t-1] + e[t-1] + b0 + b1*(x[t-1]+e[t-1]) + b2*z[t]
   # Realized value of x[t] with process error
   x[t] ~ dnorm(Ex[t],tauPro)
   }

# Output population sizes on real scale & transform from log-link for the
# Poisson distributed observations.
for (t in 2:T) {
  N[t] <- exp(x[t])
  }
}
",fill = TRUE)
sink()

# Bundle the data for use in R2jags
jags.data <- list(y = Count, T = length(Count), z = tempcov, extract = removals)
  
# Initial values
inits <- function(){list(tauPro = 1, b0 = 0.3, b1 = -0.1, b2 = 0,
  x = c(NA,log(18),rep(NA,(length(bison$Year)-2))))}
  
# Parameters monitored
parameters <- c("b0","b1","b2","sigPro","N")

# MCMC settings
ni <- 150000
nt <- 50
nb <- 100000
nc <- 3

# Call JAGS from R using R2jags
set.seed(runif(1,1,999))  # sometimes this has to be reset
gomtop <- jags(jags.data,inits,parameters,"gomtop.jags",n.chains = nc,
  n.thin = nt,n.iter = ni,n.burnin = nb,working.directory = getwd())

print(gomtop, digits=3)
traceplot(gomtop)
save(gomtop, file="gomtop.Rdata")  

###############################################################################
# Annotated code for estimating posterior model probabilities among a set of  #
# hypothesized models.                                                        #
###############################################################################
# The iterative process below will create some models not of interest (e.g.,
# a model with the interaction term but not the additive effects), and one can
# discard them afterward.
# z: the design matrix of focal variables to consider in the set of models
p <- dim(z)[2]
models <- 2^p  

sink("gomboth.jags")
cat("
model {
    
    #### Priors and constraints ###
    
    # ...
    # implement code like that above, and then following the priors for growth 
    # rate from low density and the parameter for density dependence, insert 
    # code like the following. 
    
    # Kuo & Mallick VS (with random effect) for posterior variable inclusion 
    # probabilities across the alternative environmental covariates 
    for (j in 1:p){
    benv[j] ~ dnorm(0,tauEnv)
    ind[j] ~ dbern(0.5)
    }
    tauEnv <- pow(sigEnv,-2)
    sigEnv ~ dunif(0,2)     #Using a semi-informative prior for proper mixing among models 
    # Monitor posterior model probabilities for a limited number of models
    # according to Ntzoufras 2002
    for (j in 1:p) {
    index[j] <- pow(2,j-1)
    }
    mdl <- 1 + inprod(ind[],index[])
    for (m in 1:models){
    pmdl[m] <- equals(m,mdl)
    }
    
    # The remainder of the code follows that above 
    # ...
    }
    ",fill = TRUE)
sink()

# Bundle the data for use in R2jags
jags.data <- list(y = Count, T = length(Count), models = models, z = z, p = p, extract = removals)

# Initial values
inits <- function(){list(tauPro = 1, sigEnv = 1, b0 = 0.3, b1 = -0.1, benv = rep(0,p), 
                         ind = rep(1,p), x = c(NA,log(18),rep(NA,(length(bison$Year)-2))))}

# Parameters monitored
parameters <- c("N","b0","b1","benv","sigPro","sigEnv","ind","pmdl")

# MCMC settings
ni <- 1000000
nt <- 90
nb <- 100000
nc <- 3

# Call JAGS from R
set.seed(runif(1,1,999))  # sometimes this has to be reset
gomboth <- jags(jags.data,inits,parameters,"gomboth.jags",n.chains = nc,
                n.thin = nt,n.iter = ni,n.burnin = nb,working.directory = getwd())

print(gomboth,digits=3)