inst/jags/cor_Occ_static.R

In crushing05/BayesCorrOcc: Fit Bayesian Correlated Occupancy Model

sink(file="inst/jags/cor_Occ_static.jags")
cat("
    model {

    #### Prior distributions

    ## Priors for linear indicator variables -- prob = 1 if quadratic term in model, 0.5 otherwise
    for(ii in 1:nPred){
      g[ii] ~ dbern(g[ii + nPred] + 0.5 * (1 - g[ii + nPred]))
    }


    ## Priors for quadratic indicator variables - prob = 0.5
    for(ii in (nPred + 1):(nPred * 2)){
      g[ii] ~ dbern(0.5)
    }


    ## Priors for betas - Normal(0, tau.beta) if g == 1; Normal(mu, 1/se^2) if g == 0
    for(ii in 1:(nPred * 2)){
      betaT[ii] ~ dnorm(bpriorm[ii], tprior[ii])T(-10, 10)
      beta[ii] <- g[ii] * betaT[ii]

      bpriorm[ii] <- (1 - g[ii]) * mu[ii]
      tprior[ii] <- g[ii] * tau.beta + (1 - g[ii]) * pow(se[ii], -2)
    }

    tau.beta <- pow(sigma.beta, -2)
    sigma.beta ~ dunif(0, 10)


    ## Detection priors
    alpha0 ~ dnorm(0, 0.1)T(-10, 10)
    alpha1 ~ dnorm(0, 0.1)T(-10, 10)
    alpha2 ~ dnorm(0, 0.1)T(-10, 10)

    for(ii in 1:nObs){
      eps[ii] ~ dnorm(0, tau.obs)
    }
    eps[(nObs + 1)] <- 0

    tau.obs <- pow(sigma.obs, -2)
    sigma.obs ~ dunif(0, 10)


    ## Spatial correlation priors
    xpsi[1] ~ dunif(0, 1)
    xpsi[2] ~ dunif(0, 1)


    #### GAM priors from mgcv::jagam()

    ## Parametric effect priors
    for (ii in 1:1) { b[ii] ~ dnorm(0,0.033) }


    ## prior for s(sim_dat$xy$x,sim_dat$xy$y)...
    K1 <- S1[1:29,1:29] * lambda[1]  + S1[1:29,30:58] * lambda[2]
    b[2:30] ~ dmnorm(zero[2:30],K1)


    ## smoothing parameter priors
    for (ii in 1:2) {
      lambda[ii] ~ dgamma(.05,.005)
      rho[ii] <- log(lambda[ii])
    }


    ## Logit predicted occupancy
    eta <- X %*% b + X1 %*% beta


    ### Likelihood
    for (ii in 1:nRoutes) {
      ## Detection probability
      p[ii, 1] <- 0
      logit(p[ii, 2]) <- alpha0 + alpha1 * Xp[ii] + alpha2 * nov[ii]+ eps[obs[ii]]


      ## Route-level occupancy probability
      logit(psi[ii]) <- eta[ii]


      ## Route-level occupancy
      z[ii] ~ dbern(psi[ii])


      ##  y: local availability at stop 1 -- 0 = locally unavailable,  1 = locally available
      y[ii, 1] ~ dbern(xpsi[1]/(xpsi[1] + (1 - xpsi[2])))
      h[ii, 1] ~ dbern(p[ii, z[ii] * y[ii, 1] + 1])


      ## Availability at stops 2-nStops
      for (jj in 2:nStops) {
        y[ii, jj] ~ dbern(xpsi[(z[ii] * y[ii, jj - 1] + 1)])
        h[ii, jj] ~ dbern(p[ii, z[ii] * y[ii, jj] + 1])
      } # jj
    } # ii

    }
    ", fill=TRUE)
sink()