theta.logistic.emgo.2016.jags

    model{
    # Priors and constraints
    q ~ dunif(0.01, 1) #harvest/population adjustment factor, #dunif(0.01, 1) worked
    #dbeta(1, 5) or dunif(0, 1) or dbeta(2, 5)--did not converge in 1mil.   
    c ~ dbeta(50, 150)    # c <- 0.25 for constant, distribution approx. from Dooley
    r.max ~  dlnorm(log(0.09), 0.3)T(0, 0.2) #Prior for r max, approximated from Dooley report & 
    #  pers. comm. on 5/9/16; allows some very high values to r.max > 0.3, try truncation;
    #  also tried: dunif(0.05, 0.25), dunif(0, 0.2), and dunif(0, 0.3)
    theta ~ dlnorm(log(1.76), 0.25)T(1, 3)  #theta-logistic parameter from Dooley report & pers. comm. on 5/9/16
    # allow much too high theta values > 100 in posterior, try truncation 
    #also used dunif(1,3)
    CC ~ dunif(50000, 400000) #<-- converged but gives very high CC and upper tail on N.tot
    #dunif(100000,250000) works but dlnorm(log(200000), 0.2)T(0, 400000) does not 
    sigma.proc ~ dunif(0, 0.3)                    # Prior for sd of state process
    tau.proc <- pow(sigma.proc, -2)
    #hierarchical model of missing harvest data: 1988, 2003, 2012, 2014, 2020, 2021, 2022
    m.har ~ dnorm(4081, 0.00001) #mean harvest before AMBCC season; parameterized as kill in 2016
    sigma.har ~ dunif(0, 3000) #process SD in harvest across years, shared between all years
    for(t in 1:3){ 
      tau.sur[t] <- pow(sigma.sur[t], -2)
      har[t] ~ dnorm(m.har, 1/sigma.har^2) #latent state process harvest
      har.sur[t] ~ dnorm(har[t], tau.sur[t]) #likelihood
      kill[t] <- har[t]/(1 - c) #translate harvest to kill
    }
    har[4] ~ dnorm(m.har, 1/sigma.har^2) #year 1988, latent state process harvest
    kill[4] <- har[4]/(1 - c)
    for(t in 5:18){ 
      tau.sur[t] <- pow(sigma.sur[t], -2)
      har[t] ~ dnorm(m.har, 1/sigma.har^2) #latent state process harvest
      har.sur[t] ~ dnorm(har[t], tau.sur[t]) #likelihood
      kill[t] <- har[t]/(1 - c) #translate harvest to kill
    }
    har[19] ~ dnorm(m.har, 1/sigma.har^2) #year 2003
    kill[19] <- har[19]/(1 - c)
    for(t in 20:27){ 
      tau.sur[t] <- pow(sigma.sur[t], -2)
      har[t] ~ dnorm(m.har, 1/sigma.har^2) #latent state process harvest
      har.sur[t] ~ dnorm(har[t], tau.sur[t]) #likelihood
      kill[t] <- har[t]/(1 - c) #translate harvest to kill
    }
    
    har[28] ~ dnorm(m.har, 1/sigma.har^2) #year 2012 missing,
    kill[28] <- har[28]/(1 - c)
    tau.sur[29] <- pow(sigma.sur[29], -2)
    har[29] ~ dnorm(m.har, 1/sigma.har^2) #latent state process harvest year 2013 in data,
    har.sur[29] ~ dnorm(har[29], tau.sur[29]) #likelihood
    kill[29] <- har[29]/(1 - c)
    har[30] ~ dnorm(m.har, 1/sigma.har^2) #year 2014 missing,
    kill[30] <- har[30]/(1 - c)
    
    for(t in 31:T){ #T is last data year before start of AMBCC season
      tau.sur[t] <- pow(sigma.sur[t], -2)
      har[t] ~ dnorm(m.har, 1/sigma.har^2) #latent state process harvest
      har.sur[t] ~ dnorm(har[t], tau.sur[t]) #likelihood
      kill[t] <- har[t]/(1 - c) #translate harvest to kill
    }
    mu.green ~ dlnorm(pmu.mean, pmu.tau) #prior for harvest under green policy

    # Likelihood
    # State process
    P[1] ~ dunif(0, 0.5)                      # Prior for initial population size
    P.true[1] ~ dlnorm(log(P[1]), tau.proc)
    for (t in 1:(T-1)){
    fk[t] <- (1 - exp(-0.0001*P.true[t]*CC))*kill[t]  #functional response of hunters
    P[t+1] <- max(P.true[t] + r.max*P.true[t]*(1 - P.true[t]^theta) - fk[t]/CC, 1e-5)
    P.true[t+1] ~ dlnorm(log(P[t+1]), tau.proc)
    }
    # Observation process
    for (t in 1:T) {
    logN.est[t] <- log(q) + log(P.true[t]) + log(CC)
    mu[t] <- logN.est[t] 
    tau.obs[t] <- pow(sigma.obs[t],-2)
    y[t] ~ dnorm(exp(mu[t]), tau.obs[t])
    }
    # Population sizes on real scale
    for (t in 1:T) {
    N.est[t] <- exp(logN.est[t])
    N.tot[t] <- N.est[t]/q
    }
    }