###########################################################################
###JAGS model code for Anderson, Mott, Hartman and Whiteman. 2017 Copeia###
###########################################################################
#
#Author: T.L. Anderson 
#
#Contact me (anderstl@gmail.com) if you have any questions
#
#Note: The JAGS code below largely follows examples and code in Kery (2010) and Kery and Schaub (2012)

# Abundance Models in Period 1 --------------------------------------------

#Ambystoma maculatum
sink('abun.macP1.txt')
cat("
model{
    
    #Covariate Priors
      b0~dnorm(0,0.001) #intercept
      b1.can~dnorm(0,0.001) #canopy cover
      b2.ar~dnorm(0,0.001) #pond area
      b3.talp~dnorm(0,0.001) #talpoideum density
      b4.date~dnorm(0,0.001) #sampling date
      b5.for~dnorm(0,0.001) #forested habitat surrounding pond
      b6.depth~dnorm(0,0.001) #sample depth
      b7.clust~dnorm(0,0.001) #pond clustering
    
    #Random effect priors
      tau.pond ~ dunif(0,100)
      tau.year ~ dunif(0,100)
      for(i in 1:ponds){
        eps.pond[i] ~ dnorm(0,tau.pond) #Pond random effect
      }      
      for(i in 1:years){
        eps.year[i] ~ dnorm(0,tau.year) #year random effect
      }  
      eps[i]~dnorm(0,tau) #individual random effect
      tau~dunif(0,100)
    
    # Likelihood for Abundance
    for (i in 1:R) {
      n.mac[i] ~ dpois(lambda.mac[i])
      log(lambda.mac[i]) <- vol[i]+b0+b1.can*canopy[i]+b2.ar*area[i]+b3.talp*talp[i]+b4.date*Jdate[i]+b5.for*forest[i]+b6.depth*depth[i]+b7.clust*clust[i]+eps.year[year[i]]+eps.pond[pond[i]]+eps[i]


    # Model fit assessments
    Presi[i]<-(n.mac[i]-lambda.mac[i])/sqrt(lambda.mac[i]) #pearsons residuals
    Mac.new[i]~dpois(lambda.mac[i])
    Presi.new[i]<-(Mac.new[i]-lambda.mac[i])/sqrt(lambda.mac[i])
    D[i]<-pow(Presi[i],2)
    D.new[i]<-pow(Presi.new[i],2)
    } #close i
    
    # add up discrepancy measures
    fit<-sum(D[])
    fit.new<-sum(D.new[]) 
    test<-step(fit.new-fit)
    bpvalue<-mean(test) #Bayesian p-value
}", fill=TRUE)
sink()

#Ambystoma talpoideum
sink('abun.talpP1.txt')
cat("
model {
    
    #Covariate priors
      b0~dnorm(0,0.001)
      b1.can~dnorm(0,0.001)
      b2.ar~dnorm(0,0.001)
      b3.mac~dnorm(0,0.001)
      b4.date~dnorm(0,0.001)
      b5.for ~ dnorm(0,0.001)
      b6.depth~dnorm(0,0.001)
      b7.clust~dnorm(0,0.001)

    #Random effect priors
      tau.pond ~ dunif(0,10)    
      tau.year ~ dunif(0,10)
      for(i in 1:ponds){
        eps.pond[i] ~ dnorm(0,tau.pond) #pond random effect
      }      
      for(i in 1:years){
        eps.year[i] ~ dnorm(0,tau.year) #year random effect
      }
      eps[i]~dnorm(0,tau)
      tau~dunif(0,100)

    # Likelihood for Abundance
    for (i in 1:R) {
      n.talp[i] ~ dpois(lambda.talp[i])
      log(lambda.talp[i]) <- vol[i]+b0+b1.can*canopy[i]+b2.ar*area[i]+b3.mac*mac[i]+b4.date*Jdate[i]+b5.for*forest[i]+b6.depth*depth[i]+b7.clust*clust[i]+eps.pond[pond[i]]+eps.year[year[i]]+eps[i]
      
    
    # Model fit assessments
      Presi[i]<-(n.talp[i]-lambda.talp[i])/sqrt(lambda.talp[i]) #pearsons residuals
      Talp.new[i]~dpois(lambda.talp[i])
      Presi.new[i]<-(Talp.new[i]-lambda.talp[i])/sqrt(lambda.talp[i])
      D[i]<-pow(Presi[i],2)
      D.new[i]<-pow(Presi.new[i],2)
    } #close i
    
    #add up discrepancy measures
    fit<-sum(D[])
    fit.new<-sum(D.new[])   
    test<-step(fit.new-fit)
    bpvalue<-mean(test) #bayesian p-value
}", fill=TRUE)
sink()


# Size Models in Period 1 -------------------------------------------------

#Ambystoma maculatum
sink('size.macP1.txt')
cat("
    model {
    
    #Covariate Priors
      b0~dnorm(0,0.001) #intercept
      b1.can~dnorm(0,0.001) #canopy cover
      b2.ar~dnorm(0,0.001) #pond area
      b3.talp~dnorm(0,0.001) #talpoideum density
      b4.date~dnorm(0,0.001) #sampling date
      b5.mac~dnorm(0,0.001) #maculatum density
      b6.for~dnorm(0,0.001) #forested habitat around pond
      b7.depth~dnorm(0,0.001) #sampling depth
      b8.clust~dnorm(0,0.001) #pond clustering
      tau.mac~dunif(0,100)

    #Random effect priors
      tau.pond ~ dgamma(0.001,0.001)    
      tau.year ~ dgamma(0.001,0.001)
      for(i in 1:ponds){
        eps.pond[i] ~ dnorm(0,tau.pond) #pond random effect
      }      
      for(i in 1:years){
        eps.year[i] ~ dnorm(0,tau.year) #year random effect
      } 

    # Likelihood for Size
    for (i in 1:R) {
      size.mac[i] ~ dnorm(mu.mac[i],tau.mac)
      mu.mac[i] <- b0+b1.can*canopy[i]+b2.ar*area[i]+b3.talp*talp[i]+b4.date*Jdate[i]+b5.mac*mac[i]+b6.for*forest[i]+b7.depth*depth[i]+b8.clust*clust[i]+eps.pond[pond[i]]+eps.year[year[i]]
    
    # Model fit assessments 
    residual[i]<-size.mac[i]-mu.mac[i]
    predicted[i]<-mu.mac[i]
    sq[i]<-pow(residual[i],2)
    y.new[i]~dnorm(mu.mac[i],tau.mac)
    sq.new[i]<-pow(y.new[i]-predicted[i],2)
    }

    #add up discrepancy measures
    fit<-sum(sq[])
    fit.new<-sum(sq.new[])
    test<-step(fit.new-fit)
    bpvalue<-mean(test) #bayesian p-value
    }", fill=TRUE)
sink()

#Ambystoma talpoideum
sink('size.talpP1.txt')
cat("
model {
    
    #Covarite priors
      b0~dnorm(0,0.001) #intercept
      b1.can~dnorm(0,0.001) #canopy cover
      b2.ar~dnorm(0,0.001) #pond area
      b3.talp~dnorm(0,0.001) #talpoideum density
      b4.date~dnorm(0,0.001) #sampling date
      b5.mac~dnorm(0,0.001) #maculatum density
      b6.for~dnorm(0,0.001) #forested habitat 
      b7.depth~dnorm(0,0.001) #sample depth
      b8.clust~dnorm(0,0.001) #pond clustering
      tau.talp~dunif(0,100)
    
    #Random effect priors
      tau.pond ~ dgamma(0.001,0.001)    
      tau.year ~ dgamma(0.001,0.001)
      for(i in 1:ponds){
        eps.pond[i] ~ dnorm(0,tau.pond) #pond random effect
      }       
      for(i in 1:years){
        eps.year[i] ~ dnorm(0,tau.year) #year random effect
      }  
    
    # Likelihood for Size
    for (i in 1:R) {
      size.talp[i] ~ dnorm(mu.talp[i],tau.talp)
      mu.talp[i] <- b0+b1.can*canopy[i]+b2.ar*area[i]+b3.talp*talp[i]+b4.date*Jdate[i]+b5.mac*mac[i]+b6.for*forest[i]+b7.depth*depth[i]+b8.clust*clust[i]+eps.pond[pond[i]]+eps.year[year[i]]
    
    #Model fit assessments
      residual[i]<-size.talp[i]-mu.talp[i]
      predicted[i]<-mu.talp[i]
      sq[i]<-pow(residual[i],2)
      y.new[i]~dnorm(mu.talp[i],tau.talp)
      sq.new[i]<-pow(y.new[i]-predicted[i],2)
    } #close i

    #Add up discrepancy measures
    fit<-sum(sq[])
    fit.new<-sum(sq.new[])
    test<-step(fit.new-fit)
    bpvalue<-mean(test) #bayesian p-value
}", fill=TRUE)
sink()


# Abundance Models in Period 2 --------------------------------------------
##Note the abundacne models in Period 2 are formatted for data in long-format

sink('abun.macP2.txt')
cat("model {
    
    # Detection covariate priors
      a0 ~ dunif(-10,10)
      a1.meth ~ dunif(-10,10)
    
    # Abundance covariate priors
      b0~dnorm(0,0.001) #intercept
      b1.can~dnorm(0,0.001) #canopy cover
      b2.ar~dnorm(0,0.001) #pond area
      b3.talp~dnorm(0,0.001) #talpoideum density
      b4.date~dnorm(0,0.001) #sampling date
      b5.for~dnorm(0,0.001) #forested habitat 
      b6.clust~dnorm(0,0.001) #pond clustering
    
    #Randon effect priors
      tau~dunif(0,100)
      tau.p~dunif(0,100)
      tau.pond ~ dunif(0,100)    
      tau.year ~ dunif(0,100)
      for(i in 1:n.ponds){
        eps.pond[i] ~ dnorm(0,tau.pond) #pond random effect
      }      
      for(i in 1:years){
        eps.year[i] ~ dnorm(0,tau.year) #year random effect
      }
    
    # Abundance likelihood
      for (i in 1:n.ponds) {
        n.mac[i] ~ dpois(lambda.mac[i])
        log(lambda.mac[i]) <- log(eff[i])+b0+b1.can*canopy[i]+b2.area*area[i]+b3.talp*talp[i]+b4.date*date[i]+b5.forest*forest[i]+b6.clust*clust[i]+eps.pond[pond[i]] +eps.year[year[i]]
      }

    # Detection likelihood
      for(i in 1:n.samples){
        y.mac[i] ~ dbin(p.mac[i],n.mac[sample[i]])
        p.mac[i] <- exp(lp.lim[i])/(1 + exp(lp.lim[i]))   
        lp.lim[i] <- min(999, max(-999, lp[i]))         # Help stabilize the logit
        lp[i]~dnorm(FIT[i], tau.p)                      # Random noise in detection modeled here
        FIT[i] <-a0+a1.meth*method[sample[i]] 
    
    # Assess model fit using Chi-squared discrepancy
      # Compute fit statistic for observed data
        eval[i] <- p.mac[i] * n.mac[sample[i]]
        Residuals[i] <- y.mac[i] - eval[i]
        E[i] <- pow((Residuals[i]),2) / (eval[i] + 0.5)
    
      # Generate replicate data and compute fit stats for them
        y.new[i] ~ dbin(p.mac[i], n.mac[sample[i]])
        Residuals.new[i] <- y.new[i] - eval[i]
        E.new[i] <- pow((Residuals.new[i]),2) / (eval[i] + 0.5) 
        } #close I        
    
    # Add up individual chi2 values for overall fit statistic 
      fit.actual <- sum(E[]) # Fit statistic for actual data set 
      fit.sim <- sum(E.new[]) # Fit statistic for a fitting model 
      c.hat <- fit.actual / fit.sim # c-hat estimate 
      bpv <- step(fit.actual-fit.sim) 
      bpval<-mean(bpv) # Bayesian p-value
    
}",fill=TRUE)
sink()

#Ambystoma talpoideum
sink('abun.talpP2.txt')
cat("model {
    
    # Detection covariate priors
      a0 ~ dunif(-10,10)
      a1.meth ~ dunif(-10,10)
      tau.p~dunif(0,100)

    # Abudnance covariate priors
      b0~dnorm(0,0.001) #intercept
      b1.can~dnorm(0,0.001) #canopy cover
      b2.ar~dnorm(0,0.001) #pond area
      b3.mac~dnorm(0,0.001) #maculatum density
      b4.date~dnorm(0,0.001) #sampling date
      b5.for~dnorm(0,0.001) #forested habitat 
      b6.clust~dnorm(0,0.001) #pond clustering
    
    #Random effect priors
      tau~dunif(0,100)
      tau.pond ~ dunif(0,100)    
      tau.year ~ dunif(0,100)
      for(i in 1:n.ponds){
        eps.pond[i] ~ dnorm(0,tau.pond) #pond random effect
      }      
      for(i in 1:years){
        eps.year[i] ~ dnorm(0,tau.year) #year random effect
      }
    
    # Abundance likelihood
      for (i in 1:n.ponds) {
        n.mac[i] ~ dpois(lambda.talp[i])
        log(lambda.mac[i]) <- log(eff[i])+b0+b1.can*canopy[i]+b2.area*area[i]+b3.mac*mac[i]+b4.date*date[i]+b5.forest*forest[i]+b6.clust*clust[i]+eps.pond[pond[i]] +eps.year[year[i]]
      }
    
    # Detection likelihood
      for(i in 1:n.samples){
        y.talp[i] ~ dbin(p.talp[i],n.talp[sample[i]])
        p.talp[i] <- exp(lp.lim[i])/(1 + exp(lp.lim[i]))   
        lp.lim[i] <- min(999, max(-999, lp[i]))         # Help stabilize the logit
        lp[i]~dnorm(FIT[i], tau.p)                      # Random noise in detection modeled here
        FIT[i] <-A0+A1*method[sample[i]] 
    
    # Assess model fit using Chi-squared discrepancy
      # Compute fit statistic for observed data
        eval[i] <- p.talp[i] * n.talp[sample[i]]
        Residuals[i] <- y.talp[i] - eval[i]
        E[i] <- pow((Residuals[i]),2) / (eval[i] + 0.5)
    
      # Generate replicate data and compute fit stats for them
        y.new[i] ~ dbin(p.talp[i], n.talp[sample[i]])
        Residuals.new[i] <- y.new[i] - eval[i]
        E.new[i] <- pow((Residuals.new[i]),2) / (eval[i] + 0.5) 
    } #close I        
    
    # Add up individual chi2 values for overall fit statistic 
      fit.actual <- sum(E[]) # Fit statistic for actual data set 
      fit.sim <- sum(E.new[]) # Fit statistic for a fitting model 
      c.hat <- fit.actual / fit.sim # c-hat estimate 
      bpv <- step(fit.actual-fit.sim) 
      bpval<-mean(bpv) # Bayesian p-value
}",fill=TRUE)
sink()


# Size Models for Period 2 ------------------------------------------------

#Ambystoma maculatum
sink('size.macP2.txt')
cat("
model {
    
    #Size covariate priors
      b0~dnorm(0,0.001) #intercept
      b1.can~dnorm(0,0.001) #canopy cover
      b2.ar~dnorm(0,0.001) #pond area
      b3.talp~dnorm(0,0.001) #talpoideum density
      b4.date~dnorm(0,0.001) #sampling date
      b5.mac~dnorm(0,0.001) #maculatum density
      b6.for~dnorm(0,0.001) #forested habitat around pond
      b7.clust~dnorm(0,0.001) #pond clustering
      tau.mac~dunif(0,100)

    #Random effect priors
      tau.pond ~ dgamma(0.001,0.001)    
      tau.year ~ dgamma(0.001,0.001)
      for(i in 1:ponds){
        eps.pond[i] ~ dnorm(0,tau.pond)
      }      
      for(i in 1:years){
        eps.year[i] ~ dnorm(0,tau.year)
      }  

    # Likelihood
      for (i in 1:R) {
        size.mac[i] ~ dnorm(mu.mac[i],tau.mac)
        mu.mac[i] <- b0+b1.can*canopy[i]+b2.ar*area[i]+b3.talp*talp[i]+b4.date*Jdate[i]+b5.mac*mac[i]+b6.for*forest[i]+b7.clust*clust[i]+eps.pond[pond[i]]+eps.year[year[i]]
    
      #Model fit assessment
        residual[i]<-size.mac[i]-mu.mac[i]
        predicted[i]<-mu.mac[i]
        sq[i]<-pow(residual[i],2)
        y.new[i]~dnorm(mu.mac[i],tau.mac)
        sq.new[i]<-pow(y.new[i]-predicted[i],2)
      } #close i

    #Add up discrepancy measures
    fit<-sum(sq[])
    fit.new<-sum(sq.new[])
    test<-step(fit.new-fit)
    bpvalue<-mean(test) #Bayesian P-value
}", fill=TRUE)
sink()

#Ambystoma talpoideum
sink('size.talpP2.txt')
cat("
    model {
    
    #Covariate Priors
      b0~dnorm(0,0.001) #intercept
      b1.can~dnorm(0,0.001) #canopy cover
      b2.ar~dnorm(0,0.001) #pond area
      b3.talp~dnorm(0,0.001) #talpoideum density
      b4.date~dnorm(0,0.001) #sampling date
      b5.mac~dnorm(0,0.001) #maculatum density
      b6.for~dnorm(0,0.001) #forested habitat around pond
      b7.clust~dnorm(0,0.001) #pond clustering
      tau.talp~dunif(0,100)
    
    #Random effect priors
      tau.pond ~ dgamma(0.001,0.001)    
      tau.year ~ dgamma(0.001,0.001)
      for(i in 1:ponds){
        eps.pond[i] ~ dnorm(0,tau.pond)
      }      
      for(i in 1:years){
        eps.year[i] ~ dnorm(0,tau.year)
      }  
    
    # Likelihood
      for (i in 1:R) {
        size.talp[i] ~ dnorm(mu.talp[i],tau.talp)
        mu.talp[i] <- b0+b1.can*canopy[i]+b2.ar*area[i]+b3.talp*talp[i]+b4.date*Jdate[i]+b5.mac*mac[i]+b6.for*forest[i]+b7.clust*clust[i]+eps.pond[pond[i]]+eps.year[year[i]]
    
      #Model fit assessment
        residual[i]<-size.talp[i]-mu.talp[i]
        predicted[i]<-mu.talp[i]
        sq[i]<-pow(residual[i],2)
        y.new[i]~dnorm(mu.talp[i],tau.talp)
        sq.new[i]<-pow(y.new[i]-predicted[i],2)
    } #close i
    
    #Add up discrepancy measures
      fit<-sum(sq[])
      fit.new<-sum(sq.new[])
      test<-step(fit.new-fit)
      bpvalue<-mean(test) #Bayesian P-value
}", fill=TRUE)
sink()