#===============================
#Specht_BakkenShorebirds_Model Code: 
# Multi-species SSVS occupancy code. Model code for full multi-species
#breeding pair occupancy model with stochastic search variable selection indicators.
# Final model code for assessment discarded detection parameters that were not favored by SSVS (gamma<0.5)
#===============================

#1. Prep data

dat<-read.csv("Specht_BakkenShorebirds_data.csv")


# Multispecies detection matrix
DM<-dat[,c("MAGO_BP.1","MAGO_BP.2","MAGO_BP.3","NOPI_BP.1","NOPI_BP.2","NOPI_BP.3","UPSA_BP.1","UPSA_BP.2","UPSA_BP.3","WILL_BP.1","WILL_BP.2","WILL_BP.3","WIPH_BP.1","WIPH_BP.2","WIPH_BP.3")]

#covariate matrices

#habitat suitability
suit<-cbind(dat$ST_S_MAGO, dat$ST_S_NOPI, dat$ST_S_UPSA, dat$ST_S_WILL, dat$ST_S_WIPH)

suit<-as.data.frame(suit)
#known previous occupancy by species
P.Occ<-as.data.frame(cbind(dat$early.MAGO, dat$early.NOPI, dat$early.UPSA, dat$early.WILL, dat$early.WIPH))
#survey type
survey.type<-c(0,0,1)

S=nrow(DM)
spp=5

#2. JAGS model
library(jagsUI)

sink("ms_FULL_ind.txt")
cat("
    model {
    
# Uniform priors (logit scale) for across species means
    psi.mean ~ dunif(-4, 4) # mean occupancy across species
    psi.SD ~ dunif(0.05, 2) # mean sd of occupancy across species
    psi.tau <- pow(psi.SD, -2)
    p.mean ~ dunif(-4, 4) # mean detection probability across species
    p.SD ~ dunif(0.05, 2) # mean sd of detection probability across species
    p.tau <- pow(p.SD, -2)

  # generate species specific means on logit scale
    for (j in 1:spp){ # for each species
    logit.psi.spp[j] ~ dnorm(psi.mean,psi.tau)
    logit.p.spp[j] ~ dnorm(p.mean,p.tau)

  #===== Species-specific Covariate Priors=======#

  #Occupancy Priors

    alpha0[j]<-logit.psi.spp[j] # Intercept
    alpha.Suit[j]~dunif(-20,20) # Habitat Suitability (Species specific varible)
    alpha.Traffic[j]~dunif(-20,20) # Traffic volume (site specific variable)
    alpha.Wells200[j]~dunif(-20,20) # Number of wells (site specific variable)
    alpha.Jdate[j]~dunif(-20,20) # Date to account for progression of breeding season behaviors
    alpha.PrevOcc[j] ~dunif(-20,20) # Previous known occupancy by that species at a site
    
  #Detection Priors without indicators
    beta0[j]<-logit.p.spp[j] # intercept
    beta.Survey[j]~dunif(-20,20) # indicator for transect survey type

  #Indicators on detection priors. The “gamma.___” terms indicate that a variable should be
      #included if >0.5, functioning as a similar metric to p-values based on alpha=0.5
    
    beta.Traffic[j]~dnorm(0,tau[ind.Traffic[j]]) # traffic volume
    gamma.Traffic[j]<-ind.Traffic[j]-1
    ind.Traffic[j]~ dcat(p_ind[ ])
   
    beta.Topo[j]~dnorm(0,tau[ind.Topo[j]]) #Topography of site
    gamma.Topo[j]<-ind.Topo[j]-1
    ind.Topo[j]~ dcat(p_ind[ ])
    
    beta.Wind[j]~dnorm(0,tau[ind.Wind[j]]) # ~Wind speed during survey
    gamma.Wind[j]<-ind.Wind[j]-1
    ind.Wind[j]~ dcat(p_ind[ ])
    
    beta.Temp[j]~dnorm(0,tau[ind.Temp[j]]) # Temperature at start of survey
    gamma.Temp[j]<-ind.Temp[j]-1
    ind.Temp[j]~ dcat(p_ind[ ])
    
    beta.Jdate[j]~dnorm(0,tau[ind.Jdate[j]]) # Date, to account for phenology, behavior
    gamma.Jdate[j]<-ind.Jdate[j]-1
    ind.Jdate[j]~ dcat(p_ind[ ])
    
    beta.Jdate2[j]~dnorm(0,tau[ind.Jdate2[j]]) # Date, quadratic term
    gamma.Jdate2[j]<-ind.Jdate2[j]-1
    ind.Jdate2[j]~ dcat(p_ind[ ])
    
    beta.Survey.Traffic[j]~dnorm(0,tau[ind.Survey.Traffic[j]]) #Survey type-traffic interaction
    gamma.Survey.Traffic[j]<-ind.Survey.Traffic[j]-1
    ind.Survey.Traffic[j]~ dcat(p_ind[ ])
    
    beta.ASM[j]~dnorm(0,tau[ind.ASM[j]]) # Observer
    gamma.ASM[j]<-ind.ASM[j]-1
    ind.ASM[j]~ dcat(p_ind[ ])
    
    beta.LRL[j]~dnorm(0,tau[ind.LRL[j]]) # Observer
    gamma.LRL[j]<-ind.LRL[j]-1
    ind.LRL[j]~ dcat(p_ind[ ])
    
    beta.MEW[j]~dnorm(0,tau[ind.MEW[j]]) # Observer
    gamma.MEW[j]<-ind.MEW[j]-1
    ind.MEW[j]~ dcat(p_ind[ ])
    
    beta.SMC[j]~dnorm(0,tau[ind.SMC[j]]) # Observer
    gamma.SMC[j]<-ind.SMC[j]-1
    ind.SMC[j]~ dcat(p_ind[ ])
    }
    
  # Indicator priors
    p_ind[1]<-1/2
    p_ind[2]<-1-p_ind[1]
    tau[1]<-tau_in
    tau[2]<-tau_in/1000
    tau_in<-pow(sd_beta,-2)
    sd_beta~dunif(0,100)
  
  # Likelihood =============================================
    for (i in 1:S) { # sites
    for (j in 1:spp){ # species
  
  # Occupancy Function (species specific): is site i occupied by spp j?
    z[i,j] ~ dbern(psi[i,j])
    logit(psi[i,j])<-alpha0[j] + alpha.Suit[j]*suit[i,j] + alpha.Traffic[j]*Traffic[i]
    + alpha.Wells200[j]*Wells200[i]
    + alpha.Jdate[j]*Jdate[i]+beta.Jdate2[j]*(Jdate[i]*Jdate[i])
    + alpha.PrevOcc[j]*PrevOcc[i,j]
  
  # Detection Function (species specific)
  
    for (r in 1:3) { # Surveys
  
    y[i,(j-1)*3+r] ~ dbern(z[i,j]*p1[i,(j-1)*3+r]) # prob of detecting spp
  
    logit(p1[i,(j-1)*3+r])<-beta0[j] + beta.Survey[j]*survey.type[r] +
    beta.ASM[j]*ASM[i] + beta.LRL[j]*LRL[i] + beta.MEW[j]*MEW[i] +
    beta.SMC[j]*SMC[i] +
    beta.Traffic[j]*Traffic[i] +
    beta.Topo[j]*Topo[i] + beta.Wind[j]*Wind[i] + beta.Temp[j]*Temp[i] +
    beta.Survey.Traffic[j]*survey.type[r]*Traffic[i]
  
        } # end r= Surveys
      } # end j= Species
    } # end S= Sites
    
  # Derived Quantities (back transformed occupancy and detection estimates)
    
    for(l in 1:spp){
    logit(psi.spp.pred[l])<- logit.psi.spp[l] #occupancy
    logit(p.spp.pred.pc[l])<-logit.p.spp[l] # detection- point count
    logit(p.spp.pred.tr[l])<-logit.p.spp[l]+beta.Survey[l] # detection- transect
    }

    } # end jags model
    ",fill = TRUE)

sink()


# MCMC settings 
ni <- 40000  # total iterations
nt <- 2     # thinning rate (save every nth simulation)
nb <- 2000  # discard first 2000 as burn-in (crappy estimates prior to convergence)
nc <- 3     # replicate chains (to assess convergence)
na <- 5000  # n.adapt

# Specify data, initial values, and parameters for storage
data <- list(y=DM, S=nrow(DM), spp=5, suit=suit, Traffic=dat$ST_Traffic, Wells200=dat$ST_W200, Jdate=dat$ST_Jdate, PrevOcc=P.Occ,
             survey.type=survey.type, Topo=dat$ST_Topo, Wind=dat$ST_Wind, Temp=dat$ST_Temp,
             ASM=dat$ASM, LRL=dat$LRL, MEW=dat$MEW, SMC=dat$SMC) # everything used inside jags model must be specified here


zst <- matrix(NA, nrow = S, ncol = spp)    # create matrix for initial values
for (i in 1:S){
  for (j in 1:spp){
    zst[i,j] <- ifelse(is.infinite(max(DM[i,((j-1)*3+1):((j-1)*3+3)], na.rm=T)),NA,
                       max(DM[i,((j-1)*3+1):((j-1)*3+3)], na.rm=T))# if seen, set initial value for z[i,j] to 1
  }
} 
inits <- function() list(z=zst, psi.mean=runif(1,-4,4), p.mean=runif(1,-4,4),
                         psi.SD=runif(1,0.05,2), p.SD=runif(1,0.05,2)) 
parms <- c("psi.mean", "psi.SD", "p.mean", "p.SD", "psi.spp.pred", "p.spp.pred.pc", "p.spp.pred.tr",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", "alpha.Jdate", "alpha.PrevOcc",
           "gamma.Traffic", "gamma.Topo","gamma.Time","gamma.Temp", "gamma.Wind", "gamma.Survey.Traffic",
           "gamma.ASM", "gamma.LRL", "gamma.MEW", "gamma.SMC", "gamma.Jdate", "gamma.Jdate2",
           "beta0", "beta.Survey","beta.Traffic", "beta.Topo","beta.Time", "beta.Temp", "beta.Wind", "beta.Jdate", "beta.Jdate2","beta.Survey.Traffic", 
           "beta.ASM", "beta.LRL", "beta.MEW","beta.SMC") # specify parameters to store

# Call jags from R 
ms_FULL_ind.out <- jags(data, inits, parms, "ms_FULL_ind.txt", n.chains = nc, 
                        n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)

save(ms_FULL_ind.out,file="MSOcc_ind.Rda")

# Examine output
traceplot(ms_FULL_ind.out)
print(ms_FULL_ind.out) 


############################################################
#  Full model with reduced detection structure

sink("ms_FULL_red_MAGO.txt")
cat("
    model {
    
    
    
    #Occupancy
    alpha0~dunif(-20,20)
    alpha.Suit~dunif(-20,20)
    alpha.Traffic~dunif(-20,20)
    alpha.Wells200~dunif(-20,20)
    alpha.Jdate~dunif(-20,20)
    alpha.PrevOcc~dunif(-20,20)
    
    #Detection
    beta0~dunif(-20,20)
    beta.Survey~dunif(-20,20)
    beta.Topo~dunif(-20,20)
    beta.Jdate~dunif(-20,20)
    beta.SMC~dunif(-20,20)
    
    
    
    # Likelihood
    for (i in 1:S) { # sites
    
    # Occupancy 
    z[i] ~ dbern(psi[i]) 
    logit(psi[i])<-alpha0 + alpha.Suit*suit[i] + alpha.Traffic*Traffic[i] + alpha.Wells200*Wells200[i] + 
    alpha.Jdate*Jdate[i] + alpha.PrevOcc* Prev.occ[i]
    
    
    # Detection
    
    for (r in 1:3) {  # Surveys
    y[i,r] ~ dbern(p.eff[i,r])
    p.eff[i,r]<-z[i]*p[i,r] # prob of detecting spp
    
    logit(p[i,r])<-beta0 + beta.Survey*survey.type[r] + beta.Jdate*Jdate[i]+
    beta.SMC*SMC[i] +
    beta.Topo*Topo[i] 
    
    } # end r= Surveys
    } # end S= Sites
    
    # Derived Quantities
    logit(psi.pred)<- alpha0 #back-transform
    logit(p.pred.pc)<-beta0
    logit(p.pred.tr)<-beta0+beta.Survey
    
    
    } # end jags model
    ",fill = TRUE)
sink()

# MCMC settings 
ni <- 20000  # total iterations
nt <- 2     # thinning rate (save every nth simulation)
nb <- 1000  # discard first 2000 as burn-in (crappy estimates prior to convergence)
nc <- 3     # replicate chains (to assess convergence)
na <- 4000  # 


# Specify data, initial values, and parameters for storage
data <- list(y=DM[,1:3], S=nrow(DM[,1:3]), suit=suit[,1], Traffic= dat$ST_Traffic, 
             Wells200=dat$ST_W200, Jdate=dat$ST_Jdate, 
             survey.type=survey.type, Topo=dat$ST_Topo, 
             SMC=dat$SMC ,Prev.occ=dat$early.MAGO) # everything used inside jags model must be specified here


zst <- ifelse(is.infinite(apply(DM[,1:3],1,max,na.rm=T)), NA, apply(DM[,1:3],1,max,na.rm=T))   # create matrix for initial values

inits <- function() list(z=zst, alpha0=runif(1,-3,3), beta0=runif(1,-3,3),alpha.Suit=runif(1,-3,3),
                         alpha.Traffic=runif(1,-3,3), alpha.Wells200=runif(1,-3,3),alpha.Jdate=runif(1,-3,3)) 
parms <- c("psi.pred", "p.pred.pc", "p.pred.tr",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", 
           "alpha.Jdate", "alpha.PrevOcc",
           "beta0", "beta.Survey", "beta.Topo", "beta.Jdate",
           "beta.SMC") # specify parameters to store

# Call jags from R 
ms_FULL_red_MAGO.out <- jags(data, inits, parms, "ms_FULL_red_MAGO.txt", n.chains = nc, 
                             n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)


# Examine output
traceplot(ms_FULL_red_MAGO.out)
print(ms_FULL_red_MAGO.out) 

############################################################
#  Full model with reduced detection structure

sink("ms_FULL_red_NOPI.txt")
cat("
    model {
    
    
    
    #Occupancy
    alpha0~dunif(-20,20)
    alpha.Suit~dunif(-20,20)
    alpha.Traffic~dunif(-20,20)
    alpha.Wells200~dunif(-20,20)
    alpha.Jdate~dunif(-20,20)
    alpha.PrevOcc~dunif(-20,20)
    
    #Detection
    beta0~dunif(-20,20)
    beta.Survey~dunif(-20,20)
    beta.Topo~dunif(-20,20)
    #beta.Jdate~dunif(-20,20)
    beta.Wind~dunif(-20,20)
    #beta.Temp~dunif(-20,20)
    #beta.Survey.Traffic[5]~dunif(-20,20)
    
    beta.ASM~dunif(-20,20)
    
    beta.LRL~dunif(-20,20)
    beta.SMC~dunif(-20,20)
    
    
    # Likelihood
    for (i in 1:S) { # sites
    
    # Occupancy 
    z[i] ~ dbern(psi[i]) 
    logit(psi[i])<-alpha0 + alpha.Suit*suit[i] + alpha.Traffic*Traffic[i] + alpha.Wells200*Wells200[i] 
    + alpha.Jdate*Jdate[i]+ alpha.PrevOcc*P.Occ[i]
    
    
    # Detection
    
    for (r in 1:3) {  # Surveys
    y[i,r] ~ dbern(p.eff[i,r])
    p.eff[i,r]<-z[i]*p[i,r] # prob of detecting spp
    
    logit(p[i,r])<-beta0 + beta.Survey*survey.type[r] +
    beta.ASM*ASM[i]+ beta.LRL*LRL[i] +beta.SMC*SMC[i] + 
    beta.Topo*Topo[i] + beta.Wind*Wind[i] #+ beta.Jdate*Jdate[i]
    
    } # end r= Surveys
    } # end S= Sites
    
    # Derived Quantities
    logit(psi.pred)<- alpha0 #back-transform
    logit(p.pred.pc)<-beta0
    logit(p.pred.tr)<-beta0+beta.Survey
    
    
    } # end jags model
    ",fill = TRUE)
sink()

# MCMC settings 
ni <- 50000  # total iterations
nt <- 2     # thinning rate (save every nth simulation)
nb <- 5000  # discard first 2000 as burn-in (crappy estimates prior to convergence)
nc <- 3     # replicate chains (to assess convergence)
na <- 10000  # 


# Specify data, initial values, and parameters for storage
data <- list(y=DM[,4:6], S=nrow(DM[,4:6]), suit=suit[,2], Traffic=dat$ST_Traffic, Wells200=dat$ST_W200, Jdate=dat$ST_Jdate, 
             survey.type=survey.type, Topo=dat$ST_Topo, Wind=dat$ST_Wind,P.Occ=dat$early.NOPI,
             LRL=dat$LRL, SMC=dat$SMC, ASM=dat$ASM) # everything used inside jags model must be specified here


zst <- ifelse(is.infinite(apply(DM[,4:6],1,max,na.rm=T)), NA, apply(DM[,4:6],1,max,na.rm=T))   # create matrix for initial values
inits <- function() list(z=zst, alpha0=runif(1,-3,3), beta0=runif(1,-3,3),alpha.Suit=runif(1,-3,3),
                         alpha.Traffic=runif(1,-3,3), alpha.Wells200=runif(1,-3,3),alpha.Jdate=runif(1,-3,3)) 

parms <- c("psi.pred", "p.pred.pc", "p.pred.tr",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", "alpha.Jdate", "alpha.PrevOcc",
           "beta0", "beta.Survey", "beta.Topo", "beta.Wind", #"beta.Jdate",
           "beta.ASM", "beta.LRL","beta.SMC") # specify parameters to store

# Call jags from R 
ms_FULL_red_NOPI.out <- jags(data, inits, parms, "ms_FULL_red_NOPI.txt", n.chains = nc, 
                             n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)


# Examine output
traceplot(ms_FULL_red_NOPI.out)
print(ms_FULL_red_NOPI.out) 


############################################################
#  Full model with reduced detection structure

sink("ms_FULL_red_UPSA.txt")
cat("
    model {
    
    
    
    #Occupancy
    alpha0~dunif(-20,20)
    alpha.Suit~dunif(-20,20)
    alpha.Traffic~dunif(-20,20)
    alpha.Wells200~dunif(-20,20)
    alpha.Jdate~dunif(-20,20)
    alpha.PrevOcc~dunif(-20,20)
    
    #Detection
    beta0~dunif(-20,20)
    beta.Survey~dunif(-20,20)
    #beta.Topo~dunif(-20,20)
    #beta.Traffic~dunif(-20,20)
    beta.Wind~dunif(-20,20)
    #beta.Temp~dunif(-20,20)
    #beta.Survey.Traffic[5]~dunif(-20,20)
    
    beta.ASM~dunif(-20,20)
    beta.SMC~dunif(-20,20)
    
    
    
    # Likelihood
    for (i in 1:S) { # sites
    
    # Occupancy 
    z[i] ~ dbern(psi[i]) 
    logit(psi[i])<-alpha0 + alpha.Suit*suit[i] + alpha.Traffic*Traffic[i] + alpha.Wells200*Wells200[i] + alpha.Jdate*Jdate[i]
    + alpha.PrevOcc* Prev.occ[i]
    
    # Detection
    
    for (r in 1:3) {  # Surveys
    y[i,r] ~ dbern(p.eff[i,r])
    p.eff[i,r]<-z[i]*p[i,r] # prob of detecting spp
    
    logit(p[i,r])<-beta0 + beta.Survey*survey.type[r] +
    beta.ASM*ASM[i] +beta.SMC*SMC[i] + 
    beta.Wind*Wind[i] 
    
    } # end r= Surveys
    } # end S= Sites
    
    # Derived Quantities
    logit(psi.pred)<- alpha0 #back-transform
    logit(p.pred.pc)<-beta0
    logit(p.pred.tr)<-beta0+beta.Survey
    
    
    } # end jags model
    ",fill = TRUE)
sink()

# MCMC settings 
ni <- 20000  # total iterations
nt <- 2     # thinning rate (save every nth simulation)
nb <- 1000  # discard first 2000 as burn-in (crappy estimates prior to convergence)
nc <- 3     # replicate chains (to assess convergence)
na <- 4000  # 


# Specify data, initial values, and parameters for storage
data <- list(y=DM[,7:9], S=nrow(DM[,7:9]), suit=suit[,3], Traffic=dat$ST_Traffic, Wells200=dat$ST_W200, Jdate=dat$ST_Jdate, 
             survey.type=survey.type, Wind=dat$ST_Wind,
             ASM=dat$ASM, SMC=dat$SMC, Prev.occ=dat$early.UPSA) # everything used inside jags model must be specified here


zst <- ifelse(is.infinite(apply(DM[,7:9],1,max,na.rm=T)), NA, apply(DM[,7:9],1,max,na.rm=T))   # create matrix for initial values

inits <- function() list(z=zst, alpha0=runif(1,-3,3), beta0=runif(1,-3,3),alpha.Suit=runif(1,-3,3),
                         alpha.Traffic=runif(1,-3,3), alpha.Wells200=runif(1,-3,3),alpha.Jdate=runif(1,-3,3)) 
parms <- c("psi.pred", "p.pred.pc", "p.pred.tr",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", "alpha.Jdate", "alpha.PrevOcc",
           "beta0", "beta.Survey", "beta.Wind", 
           "beta.ASM","beta.SMC") # specify parameters to store

# Call jags from R 
ms_FULL_red_UPSA.out <- jags(data, inits, parms, "ms_FULL_red_UPSA.txt", n.chains = nc, 
                             n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)


# Examine output
traceplot(ms_FULL_red_UPSA.out)
print(ms_FULL_red_UPSA.out) 



############################################################
#  Full model with reduced detection structure

sink("ms_FULL_red_WILL.txt")
cat("
    model {
    
    #Occupancy
    alpha0~dunif(-20,20)
    alpha.Suit~dunif(-20,20)
    alpha.Traffic~dunif(-20,20)
    alpha.Wells200~dunif(-20,20)
    alpha.Jdate~dunif(-20,20)
    #alpha.PrevOcc~dunif(-20,20)
    
    #Detection
    beta0~dunif(-20,20)
    beta.Survey~dunif(-20,20)
    beta.Topo~dunif(-20,20)
    #beta.Jdate~dunif(-20,20)
    #beta.Jdate2~dunif(-20,20)
    #beta.Traffic~dunif(-20,20)
    #beta.Wind~dunif(-20,20)
    #beta.Temp~dunif(-20,20)
    #beta.Survey.Traffic[5]~dunif(-20,20)
    
    #beta.ASM~dunif(-20,20)
    
    beta.LRL~dunif(-20,20)
    beta.SMC~dunif(-20,20)
    
    
    
    # Likelihood
    for (i in 1:S) { # sites
    
    # Occupancy 
    z[i] ~ dbern(psi[i]) 
    logit(psi[i])<-alpha0 + alpha.Suit*suit[i] + alpha.Traffic*Traffic[i] + alpha.Wells200*Wells200[i] 
    + alpha.Jdate*Jdate[i] #+ alpha.PrevOcc*PrevOcc[i]
    
    # Detection
    
    for (r in 1:3) {  # Surveys
    y[i,r] ~ dbern(p.eff[i,r])
    p.eff[i,r]<-z[i]*p[i,r] # prob of detecting spp
    
    logit(p[i,r])<-beta0 + beta.Survey*survey.type[r] +
    beta.LRL*LRL[i] +beta.SMC*SMC[i] + 
    beta.Topo*Topo[i] 
    #+ beta.Jdate*Jdate[i]+ beta.Jdate2*(Jdate[i]*Jdate[i])
    
    } # end r= Surveys
    } # end S= Sites
    
    # Derived Quantities
    logit(psi.pred)<- alpha0 #back-transform
    logit(p.pred.pc)<-beta0
    logit(p.pred.tr)<-beta0+beta.Survey
    cum.p<-p.pred.pc+((1-p.pred.pc)*p.pred.pc)+ ((1-(p.pred.pc+((1-p.pred.pc)*p.pred.pc)))*p.pred.tr)
    
    
    } # end jags model
    ",fill = TRUE)
sink()

# MCMC settings 
ni <- 25000  # total iterations
nt <- 2     # thinning rate (save every nth simulation)
nb <- 3000  # discard first 2000 as burn-in (crappy estimates prior to convergence)
nc <- 3     # replicate chains (to assess convergence)
na <- 4000  # 


# Specify data, initial values, and parameters for storage
data <- list(y=DM[,10:12], S=nrow(DM[,10:12]), suit=suit[,4], Traffic=dat$ST_Traffic, Wells200=dat$ST_W200, Jdate=dat$ST_Jdate, 
             survey.type=survey.type, Topo=dat$ST_Topo,#PrevOcc=dat$early.WILL,
             LRL=dat$LRL, SMC=dat$SMC) # everything used inside jags model must be specified here


zst <- ifelse(is.infinite(apply(DM[,10:12],1,max,na.rm=T)), NA, apply(DM[,10:12],1,max,na.rm=T))   # create matrix for initial values

inits <- function() list(z=zst, alpha0=runif(1,-3,3), beta0=runif(1,-3,3),alpha.Suit=runif(1,-3,3),
                         alpha.Traffic=runif(1,-3,3), alpha.Wells200=runif(1,-3,3),alpha.Jdate=runif(1,-3,3)) 

parms <- c("psi.pred", "p.pred.pc", "p.pred.tr","cum.p",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", "alpha.Jdate",#"alpha.PrevOcc",
           "beta0", "beta.Survey", "beta.Topo", #"beta.Jdate", "beta.Jdate2",
           "beta.LRL","beta.SMC") # specify parameters to store

# Call jags from R 
ms_FULL_red_WILL.out <- jags(data, inits, parms, "ms_FULL_red_WILL.txt", n.chains = nc, 
                             n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)


# Examine output
traceplot(ms_FULL_red_WILL.out)
print(ms_FULL_red_WILL.out) 

############################################################
#  Full model with reduced detection structure

sink("ms_FULL_red_WIPH.txt")
cat("
    model {
    
    #Occupancy
    alpha0~dunif(-20,20)
    alpha.Suit~dunif(-20,20)
    alpha.Traffic~dunif(-20,20)
    alpha.Wells200~dunif(-20,20)
    alpha.Jdate~dunif(-20,20)
    alpha.PrevOcc~dunif(-20,20)
    
    #Detection
    beta0~dunif(-20,20)
    beta.Survey~dunif(-20,20)
    beta.Topo~dunif(-20,20)
    beta.Traffic~dunif(-20,20)
    beta.Time~dunif(-20,20)
    #beta.Temp~dunif(-20,20)
    
    #beta.ASM~dunif(-20,20)
    
    beta.LRL~dunif(-20,20)
    #beta.SMC~dunif(-20,20)
    
    
    
    # Likelihood
    for (i in 1:S) { # sites
    
    # Occupancy 
    z[i] ~ dbern(psi[i]) 
    logit(psi[i])<-alpha0 + alpha.Suit*suit[i] + alpha.Traffic*Traffic[i] + alpha.Wells200*Wells200[i] + alpha.Jdate*Jdate[i]
    + alpha.PrevOcc* Prev.occ[i]
    
    # Detection
    
    for (r in 1:3) {  # Surveys
    y[i,r] ~ dbern(p.eff[i,r])
    p.eff[i,r]<-z[i]*p[i,r] # prob of detecting spp
    
    logit(p[i,r])<-beta0 + beta.Survey*survey.type[r] +
    beta.LRL*LRL[i] +
    beta.Topo*Topo[i] + beta.Traffic*Traffic[i] + beta.Time*Time[i]
    
    } # end r= Surveys
    } # end S= Sites
    
    # Derived Quantities
    logit(psi.pred)<- alpha0 #back-transform
    logit(p.pred.pc)<-beta0
    logit(p.pred.tr)<-beta0+beta.Survey
    cum.p<-p.pred.pc+((1-p.pred.pc)*p.pred.pc)+ ((1-(p.pred.pc+((1-p.pred.pc)*p.pred.pc)))*p.pred.tr)
    
    
    } # end jags model
    ",fill = TRUE)
sink()

# MCMC settings 
ni <- 20000  # total iterations
nt <- 2     # thinning rate (save every nth simulation)
nb <- 1000  # discard first 2000 as burn-in (crappy estimates prior to convergence)
nc <- 3     # replicate chains (to assess convergence)
na <- 4000  # 


# Specify data, initial values, and parameters for storage
data <- list(y=DM[,13:15], S=nrow(DM[,13:15]), suit=suit[,5], Traffic=dat$ST_Traffic, Wells200=dat$ST_W200, Jdate=dat$ST_Jdate, 
             survey.type=survey.type, Topo=dat$ST_Topo, Time=dat$ST_Time,
             LRL=dat$LRL, Prev.occ=dat$early.WIPH) # everything used inside jags model must be specified here


zst <- ifelse(is.infinite(apply(DM[,13:15],1,max,na.rm=T)), NA, apply(DM[,13:15],1,max,na.rm=T))   # create matrix for initial values

inits <- function() list(z=zst, alpha0=runif(1,-3,3), beta0=runif(1,-3,3),alpha.Suit=runif(1,-3,3),
                         alpha.Traffic=runif(1,-3,3), alpha.Wells200=runif(1,-3,3),alpha.Jdate=runif(1,-3,3)) 
parms <- c("psi.pred", "p.pred.pc", "p.pred.tr","cum.p",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", "alpha.Jdate", "alpha.PrevOcc",
           "beta0", "beta.Survey", "beta.Topo", "beta.Traffic", "beta.Time",
           "beta.LRL") # specify parameters to store

# Call jags from R 
ms_FULL_red_WIPH.out <- jags(data, inits, parms, "ms_FULL_red_WIPH.txt", n.chains = nc, 
                             n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)


# Examine output
traceplot(ms_FULL_red_WIPH.out)
print(ms_FULL_red_WIPH.out) 




#===================================================================================
#  BROOD MODELS
#===================================================================================

dat.b<-read.csv("Specht_BakkenBrood_data.csv")
suit.b<-cbind(dat.b$B.ST_S_MAGO, dat.b$B.ST_S_NOPI, dat.b$B.ST_S_UPSA, dat.b$B.ST_S_WILL, dat.b$B.ST_S_WIPH)

suit.b<-as.data.frame(suit.b)

survey.type<-c(0,0,1)

Prev.Occ.B<-cbind(dat.b$early.MAGO, dat.b$early.NOPI, dat.b$early.UPSA, dat.b$early.WILL, dat.b$early.WIPH)

DM.b<-dat.b[,c("MAGO_Brood.1","MAGO_Brood.2","MAGO_Brood.3","NOPI_Brood.1","NOPI_Brood.2","NOPI_Brood.3", 
               "UPSA_Brood.1","UPSA_Brood.2","UPSA_Brood.3","WILL_Brood.1","WILL_Brood.2","WILL_Brood.3", "WIPH_Brood.1","WIPH_Brood.2","WIPH_Brood.3")]
DM.Ab<-dat.b[,c("All_Brood.1","All_Brood.2","All_Brood.3")] # all broods

spp=5
S=nrow(DM.b)



#===============================================
# Multi-species brood occupancy model code. Model code for full multispecies
#brood occupancy model with shared detection and occupancy covariates retained
#AFTER stochastic search variable selection.
#===============================================


sink("ms_FULL_Broods.txt")
cat("
    model {
    
# Uniform priors (logit scale)

  # Mean occupancy and detection estimates across species
    psi.mean ~ dunif(-4, 4)
    psi.SD ~ dunif(0.05, 2)
    psi.tau <- pow(psi.SD, -2)
    p.mean ~ dunif(-4, 4)
    p.SD ~ dunif(0.05, 2)
    p.tau <- pow(p.SD, -2)

  # generate species specific means on logit scale
    for (j in 1:spp){
    logit.psi.spp[j] ~ dnorm(psi.mean,psi.tau)

  # ===Species-specific Covariate Priors ==============

    #Occupancy
    alpha0[j]<-logit.psi.spp[j] # Intercept

  #Detection
    beta.Survey[j]~dunif(-20,20) # indicator for transect survey type
      }
    
    logit.p~ dnorm(p.mean,p.tau)
    beta0<-logit.p

  # == Non species specific priors==
    
    alpha.Suit~dunif(-20,20) # suitability- species specific variable, shared coefficient
    alpha.Traffic~dunif(-20,20) # Traffic volume
    alpha.Wells200~dunif(-20,20) # number of wells
    alpha.Jdate~dunif(-20,20) # date
    alpha.PrevOcc~dunif(-20,20) # PrevOcc
    beta.Topo~dunif(-20,20) # topography
    beta.LRL~dunif(-20,20) # observe
    
  # Likelihood===================================
    for (i in 1:S) { # sites
    
    for (j in 1:spp){ # species
  
  # Occupancy
    z[i,j] ~ dbern(psi[i,j]) # is site occupied by spp j?
    logit(psi[i,j])<-alpha0[j] + alpha.Suit*suit[i,j] + alpha.Traffic*Traffic[i] +
    alpha.Wells200*Wells200[i] + alpha.Jdate*Jdate[i] + alpha.PrevOcc*PrevOcc[i,j]

  # Detection
    for (r in 1:3) { # Surveys

    y[i,(j-1)*3+r] ~ dbern(z[i,j]*p1[i,(j-1)*3+r]) # prob of detecting spp
    logit(p1[i,(j-1)*3+r])<-beta0 + beta.Survey[j]*survey.type[r] +
    beta.Topo*Topo[i] + beta.LRL*LRL[i]

        } # end r= Surveys
      } # end j= Species
    } # end S= Sites
    
  # Derived Quantities
  
  for(l in 1:spp){
    logit(psi.spp.pred[l])<- logit.psi.spp[l] #back-transformed occupancy
    logit(p.spp.pred.tr[l])<-logit.p+beta.Survey[l] #back-transformed transect detection
    }
    logit(p.pred.pc)<-logit.p #back-transformed point count detection

    } # end jags model
    ",fill = TRUE)

    sink()
    # MCMC settings
    ni <- 50000 # total iterations
    nt <- 2 # thinning rate (save every nth simulation)
    nb <- 5000 # discard first 2000 as burn-in (crappy estimates prior to convergence)
    nc <- 3 # replicate chains (to assess convergence)
    na <- 1000 # n.adapt

    # Specify data, initial values, and parameters for storage
    data <- list(y=DM.b, S=nrow(DM.b), spp=5, suit=suit.b, Traffic=dat.b$B.ST_Traffic, 
                 Wells200=dat.b$B.ST_W200, Jdate=dat.b$B.ST_Jdate, 
                 survey.type=survey.type, Topo=dat.b$B.ST_Topo, LRL=dat.b$LRL, PrevOcc=Prev.Occ.B) # everything used inside jags model must be specified here
    
    
    #Initial values
zst <- matrix(NA, nrow = nrow(DM.b), ncol = spp) # create matrix for initial values
for (i in 1:nrow(DM.b)){
  for (j in 1:spp){
    zst[i,j] <- ifelse(is.infinite(max(DM.b[i,((j-1)*3+1):((j-1)*3+3)], na.rm=T)),NA,
                       max(DM.b[i,((j-1)*3+1):((j-1)*3+3)], na.rm=T))# if seen, set initial value for z[i,j] to 1
  }
}
inits <- function() list(z=zst, psi.mean=runif(1,-4,4), p.mean=runif(1,-4,4),
                         psi.SD=runif(1,0.05,2), p.SD=runif(1,0.05,2)) 
parms <- c("psi.mean", "psi.SD", "p.mean", "p.SD", "psi.spp.pred", "p.spp.pred.pc", "p.spp.pred.tr",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", "alpha.Jdate", "alpha.PrevOcc",
           "gamma.Traffic", "gamma.Topo","gamma.Temp", "gamma.Wind", "gamma.Time","gamma.Jdate",
           "gamma.ASM", "gamma.LRL",  "gamma.SMC","gamma.PrevOcc",
           "beta0", "beta.Survey","beta.Traffic", "beta.Topo", "beta.Temp", "beta.Wind", "beta.Time","beta.Jdate",
           "beta.ASM", "beta.LRL", "beta.SMC") # specify parameters to store

# Call jags from R 
ms_FULL_Broods.out <- jags(data, inits, parms, "ms_FULL_Broods.txt", n.chains = nc, 
                                n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)

# Examine output
traceplot(ms_FULL_Broods.out)
print(ms_FULL_Broods.out) 



#==============================
#  Brood models considering species-specific differences in Traffic intensity
#================================

#==============================================================
sink("ms_FULL_Broods_Traffic.txt")
cat("
    model {
    
    # Uniform priors (logit scale)
    psi.mean ~ dunif(-4, 4)
    psi.SD ~ dunif(0.05, 2) 
    psi.tau <- pow(psi.SD, -2)
    
    p.mean ~ dunif(-4, 4)
    p.SD ~ dunif(0.05, 2)
    p.tau <- pow(p.SD, -2)
    
    # generate species specific means on logit scale and backtransform
    for (j in 1:5){
    logit.psi.spp[j] ~ dnorm(psi.mean,psi.tau)
    
    
    # Species-specific Covariate Priors
    
    #Occupancy
    alpha0[j]<-logit.psi.spp[j]
    alpha.Traffic[j]~dunif(-20,20)
    
    #Detection
    beta.Survey[j]~dunif(-20,20)
    
    }
    
    logit.p~ dnorm(p.mean,p.tau)
    
    beta0<-logit.p
    
    alpha.Suit~dunif(-20,20)
    alpha.Wells200~dunif(-20,20)
    alpha.Jdate~dunif(-20,20)
    alpha.PrevOcc~dunif(-20,20)
    
    beta.Topo~dunif(-20,20)
    beta.LRL~dunif(-20,20)
    
    # Likelihood
    for (i in 1:S) { # sites
    for (j in 1:5){ # 5 species
    
    # Occupancy 
    z[i,j] ~ dbern(psi[i,j]) # is site occupied by spp j?
    logit(psi[i,j])<-alpha0[j] + alpha.Suit*suit[i,j] + alpha.Traffic[j]*Traffic[i] + 
    alpha.Wells200*Wells200[i] + alpha.Jdate*Jdate[i] + alpha.PrevOcc*PrevOcc[i,j]
    
    
    # Detection
    
    for (r in 1:3) {  # Surveys
    y[i,(j-1)*3+r] ~ dbern(z[i,j]*p1[i,(j-1)*3+r]) # prob of detecting spp
    logit(p1[i,(j-1)*3+r])<-beta0 + beta.Survey[j]*survey.type[r] +
    beta.Topo*Topo[i] + beta.LRL*LRL[i] 
    
    } # end r= Surveys
    } # end j= Species
    } # end S= Sites
    
    # Derived Quantities
    for(l in 1:5){
    logit(psi.spp.pred[l])<- logit.psi.spp[l] #back-transform
    logit(p.spp.pred.tr[l])<-logit.p+beta.Survey[l]
    }
    logit(p.pred.pc)<-logit.p
    
    } # end jags model
    ",fill = TRUE)
sink()

# MCMC settings 
# MCMC settings 
ni <- 50000  # total iterations
nt <- 2     # thinning rate (save every nth simulation)
nb <- 5000  # discard first 2000 as burn-in (crappy estimates prior to convergence)
nc <- 3     # replicate chains (to assess convergence)
na <- 7000  # n.adapt

# Specify data, initial values, and parameters for storage
data <- list(y=DM.b, S=nrow(DM.b),  suit=suit.b, Traffic=dat.b$B.ST_Traffic, Wells200=dat.b$B.ST_W200, Jdate=dat.b$B.ST_Jdate, 
             survey.type=survey.type, Topo=dat.b$B.ST_Topo, PrevOcc=Prev.Occ.B, LRL=dat.b$LRL) # everything used inside jags model must be specified here


zst <- matrix(NA, nrow = nrow(DM.b), ncol = 5)    # create matrix for initial values
for (i in 1:nrow(DM.b)){
  for (j in 1:spp){
    zst[i,j] <- ifelse(is.infinite(max(DM.b[i,((j-1)*3+1):((j-1)*3+3)], na.rm=T)),NA,
                       max(DM.b[i,((j-1)*3+1):((j-1)*3+3)], na.rm=T))# if seen, set initial value for z[i,j] to 1
  }
}
inits <- function() list(z=zst, psi.mean=runif(1,-4,4), p.mean=runif(1,-4,4),
                         psi.SD=runif(1,0.05,2), p.SD=runif(1,0.05,2)) 
parms <- c("psi.mean", "psi.SD", "p.mean", "p.SD", "psi.spp.pred", "p.pred.pc", "p.spp.pred.tr",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", "alpha.Jdate", "alpha.PrevOcc",
           "beta0", "beta.Survey","beta.Topo", "beta.LRL") # specify parameters to store

# Call jags from R 
ms_FULL_Broods_Traffic.out <- jags(data, inits, parms, "ms_FULL_Broods_Traffic.txt", n.chains = nc, 
                                         n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)

# Examine output 
traceplot(ms_FULL_Broods_Traffic.out)
print(ms_FULL_Broods_Traffic.out) 

#=====================
# Brood model for looking at species-specific well density effects
#=====================

sink("ms_Broods_Wells.txt")

cat("
    model {
    
    # Uniform priors (logit scale)
    psi.mean ~ dunif(-4, 4)
    psi.SD ~ dunif(0.05, 2) 
    psi.tau <- pow(psi.SD, -2)
    
    p.mean ~ dunif(-4, 4)
    p.SD ~ dunif(0.05, 2)
    p.tau <- pow(p.SD, -2)
    
    # generate species specific means on logit scale and backtransform
    for (j in 1:spp){
    logit.psi.spp[j] ~ dnorm(psi.mean,psi.tau)
    
    
    # Species-specific Covariate Priors
    
    #Occupancy
    alpha0[j]<-logit.psi.spp[j]
    alpha.Wells200[j]~dunif(-20,20)
    
    #Detection
    beta.Survey[j]~dunif(-20,20)
    
    }
    
    logit.p~ dnorm(p.mean,p.tau)
    
    beta0<-logit.p
    
    alpha.Traffic~dunif(-20,20)
    
    alpha.Suit~dunif(-20,20)
    alpha.Jdate~dunif(-20,20)
    alpha.PrevOcc~dunif(-20,20)
    
    beta.Topo~dunif(-20,20)
    beta.LRL~dunif(-20,20)
    # Likelihood
    for (i in 1:S) { # sites
    for (j in 1:spp){
    
    # Occupancy 
    z[i,j] ~ dbern(psi[i,j]) # is site occupied by spp j?
    logit(psi[i,j])<-alpha0[j] + alpha.Suit*suit[i,j] + alpha.Traffic*Traffic[i] + 
    alpha.Wells200[j]*Wells200[i] + alpha.Jdate*Jdate[i] +alpha.PrevOcc*PrevOcc[i,j]
    
    
    # Detection
    
    for (r in 1:3) {  # Surveys
    y[i,(j-1)*3+r] ~ dbern(z[i,j]*p1[i,(j-1)*3+r]) # prob of detecting spp
    logit(p1[i,(j-1)*3+r])<-beta0 + beta.Survey[j]*survey.type[r] +
    beta.Topo*Topo[i] + beta.LRL*LRL[i] 
    
    } # end r= Surveys
    } # end j= Species
    } # end S= Sites
    
    # Derived Quantities
    for(l in 1:spp){
    logit(psi.spp.pred[l])<- logit.psi.spp[l] #back-transform
    logit(p.spp.pred.tr[l])<-logit.p+beta.Survey[l]
    }
    logit(p.pred.pc)<-logit.p
    
    } # end jags model
    ",fill = TRUE)
sink()

# MCMC settings 
ni <- 50000  # total iterations
nt <- 2     # thinning rate (save every nth simulation)
nb <- 5000  # discard first 2000 as burn-in (crappy estimates prior to convergence)
nc <- 3     # replicate chains (to assess convergence)
na <- 7000  # n.adapt

# Specify data, initial values, and parameters for storage
data <- list(y=DM.b, S=nrow(DM.b), spp=5, suit=suit.b, Traffic=dat.b$B.ST_Traffic, Wells200=dat.b$B.ST_W200, Jdate=dat.b$B.ST_Jdate, 
             survey.type=survey.type, Topo=dat.b$B.ST_Topo, PrevOcc=Prev.Occ.B, LRL=dat.b$LRL) # everything used inside jags model must be specified here


zst <- matrix(NA, nrow = nrow(DM.b), ncol = spp)    # create matrix for initial values
for (i in 1:nrow(DM.b)){
  for (j in 1:spp){
    zst[i,j] <- ifelse(is.infinite(max(DM.b[i,((j-1)*3+1):((j-1)*3+3)], na.rm=T)),NA,
                       max(DM.b[i,((j-1)*3+1):((j-1)*3+3)], na.rm=T))# if seen, set initial value for z[i,j] to 1
  }
}
inits <- function() list(z=zst, psi.mean=runif(1,-4,4), p.mean=runif(1,-4,4),
                         psi.SD=runif(1,0.05,2), p.SD=runif(1,0.05,2)) 
parms <- c("psi.mean", "psi.SD", "p.mean", "p.SD", "psi.spp.pred", "p.pred.pc", "p.spp.pred.tr",
           "alpha0", "alpha.Suit", "alpha.Traffic","alpha.Wells200", "alpha.Jdate", "alpha.PrevOcc",
           "beta0", "beta.Survey","beta.Topo", "beta.LRL") # specify parameters to store

# Call jags from R 
Broods_Wells.out <- jags(data, inits, parms, "ms_Broods_Wells.txt", n.chains = nc, 
                                       n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt=na)

# Examine output 
traceplot(Broods_Wells.out)
print(Broods_Wells.out)