##################### JAGS model code for Andean Boreal Migrants #############################################
### The code below is based on and closely follows the code in Sollmann et al. 2016. It was originally developed for line-transect
### distance sampling but is adapted here for point-count distance sampling. This version of the model includes linear relationships
### between abundance and the environmental covariates (climate, woody biomass).


#### Data to be provided to the model are (see attached file "R code to accompany Andean community model"):
##	spec = number of species in the community
##	OBSVAR = categorical observer covariate, here with 4 levels for the 4 observers
##	time = time of day continuous detectability covariate calculated as minutes from 0600
##	height = estimate of mean canopy height at each survey point
##	nG = number of distance intervals, here 8 based on 5m bins from 0 to 40m 
##	db = vector with break points of distance intevals, here (0,5,10,...40)
##	v = width of distance intervals, here constant at 5m
##	pi = proportion of area of each 5m band
##	VAR = environmental covariates (include land cover category, woody biomass at 1km radius, elevation, latitude, climate covariates)
##      nVAR = number of environmental covariates
##	y = matrix of species observations by site (sites as rows, species as columns)
##	nind = total number of observations across all species and sites
##	dclass = vector with distance class for each observation
##	species = vector with species index for each observation
##	site = vector with site index for each observation


model{
for (s in 1:spec){
	asig[s]~dnorm(mu_s, tau_s) # asig=species level detection intercept
	for (k in 1:nVAR){
	beta[s,k]~dnorm(mu_b[k],tau_b[k])
	}
	#beta[s,1]<-0
	alpha[s]~dnorm(mu_a,tau_a)
	csig[s]~dnorm(0,0.01) # prior for time of day effect on detectability
	dsig[s]~dnorm(0,0.01) # prior for canopy height effect on detectability
}

# community hyperpriors for the detection(s) and abundance(a) components of the model
mu_s~dnorm(0,0.01)
tau_s<-1/(sig_s*sig_s)
sig_s~dunif(0,20)
mu_a~dnorm(0,0.01)
sig_a<-1/sqrt(tau_a)
tau_a~dgamma(0.1,0.1)

for (k in 1:nVAR){
	mu_b[k]~dnorm(0,0.01)
	sig_b[k]<-1/sqrt(tau_b[k])
	tau_b[k]~dgamma(0.1,0.1)
}

# observer covariates - assumed to be a common effect across species
for (k in 2:4){
	bsig[k]~dnorm(0,0.01) # observer effect - categorical covariate with 4 levels
}

bsig[1]<-0  ##set beta for category 1 = 0

for (s in 1:spec){
	for (j in 1:nsites){
	sigma[s,j]<-exp(asig[s] + bsig[OBSVAR[j]] + csig[s]*time[j] + dsig[s]*height[j]) # log-linear predictor of detection parameter sigma
	f.0[s,j] <- 2 * dnorm(0,0, 1/sigma[s,j]^2)

### construct detection probabilities by distance class w/ half-normal model

  	for(k in 1:nG){
       		up[s,j,k]<-pnorm(db[k+1], 0, 1/sigma[s,j]^2) 
       		low[s,j,k]<-pnorm(db[k], 0, 1/sigma[s,j]^2) 
		p[s,j,k]<- 2 * ( up[s,j,k] - low[s,j,k])
   		f[s,j,k]<- p[s,j,k]/f.0[s,j]/v # detection prob. in distance category k                      
    		fc[s,j,k]<- f[s,j,k] * pi[k] # pi=percent area of k
    		fct[s,j,k]<-fc[s,j,k]/sum(fc[s,j,1:nG])  
  	}

 	pcap[s,j]<-sum(fc[s,j,1:nG]) # overall detection probability

     	    lambda[j,s]<- exp(alpha[s] + inprod(beta[s,],VAR[j,])) # log-linear predictor of abundance parameter lambda
     	    y[j,s]~ dbin(pcap[s,j],N[j,s])
     	    N[j,s]~dpois(lambda[j,s])

	    ###create replicate abundances for fit statistics, species by site
	    Nnew[j,s]~dpois(lambda[j,s])

	    ### calculate species and site specific residuals
	    FT1[j,s]<-pow(sqrt(N[j,s])-sqrt(lambda[j,s]),2)
	    FT1new[j,s]<-pow(sqrt(Nnew[j,s])-sqrt(lambda[j,s]),2)
	}
        ##sum residuals over sites
	T1p[s]<-sum(FT1[1:nsites,s])
	T1newp[s]<-sum(FT1new[1:nsites,s])
}

# Bayesian p-value
Bp.N<-sum(T1newp[1:spec])>sum(T1p[1:spec])

################################################################################################
### observation model

      for(i in 1:nind){
      	dclass[i] ~ dcat(fct[species[i],site[i],1:nG]) 

      	#### fit statistic for observation model
      	### draw new data set
     	dclassnew[i] ~ dcat(fct[species[i],site[i],1:nG]) 
      	Tobsp[i]<- pow(1- sqrt(fct[species[i],site[i],dclass[i]]),2)
      	Tobspnew[i]<- pow(1- sqrt(fct[species[i],site[i],dclassnew[i]]),2)
      }

Bp.Obs<-sum(Tobspnew[1:nind])>sum(Tobsp[1:nind])


################################################################################################
###derived parameters - total abundance and predicted abundance across land cover and elevation

for (i in 1:spec){
	Nspec[i]<-sum(N[1:nsites,i]) # total abundance across all sites, avg per site would be Nspec[i]/nsites 
	## code below calculates the predicted abundance of each species at each elevation and in each land cover class
	for (j in 1:56){
		MF.eff[j,i]<-exp(alpha[i] + beta[i,7]*m[j] + beta[i,8]*m2[j]) # mature forest
		SG.eff[j,i]<-exp(alpha[i] + beta[i,1] + beta[i,7]*m[j] + beta[i,8]*m2[j]) # secondary forest
		SC.eff[j,i]<-exp(alpha[i] + beta[i,2] + beta[i,7]*m[j] + beta[i,8]*m2[j]) # shade coffee
		SUN.eff[j,i]<-exp(alpha[i] + beta[i,3] + beta[i,7]*m[j] + beta[i,8]*m2[j]) # sun coffee
		FFRS.eff[j,i]<-exp(alpha[i] + beta[i,4] + beta[i,7]*m[j] + beta[i,8]*m2[j]) # forest fragments
		LFP.eff[j,i]<-exp(alpha[i] + beta[i,5] + beta[i,7]*m[j] + beta[i,8]*m2[j]) # live fences
		PAR.eff[j,i]<-exp(alpha[i] + beta[i,6] + beta[i,7]*m[j] + beta[i,8]*m2[j]) # paramo
		} 
	## following code calculates the average abundance in each land cover class across 1000-2000m for all classes except paramo (estimated>3275m)
	Mature.forest[i]<-mean(MF.eff[12:29,i]) 
	Second.growth[i]<-mean(SG.eff[12:29,i])
	Shade.coffee[i]<-mean(SC.eff[12:29,i]) 
	Sun.coffee[i]<-mean(SUN.eff[12:29,i])
 	Fragment[i]<-mean(FFRS.eff[12:29,i]) 
	Live.fence[i]<-mean(LFP.eff[12:29,i])
	Paramo[i]<-mean(PAR.eff[49:56,i])
	}
# Total community density in each land cover class for all 23 species
Mature<-sum(Mature.forest[1:23])
Sec.Growth<-sum(Second.growth[1:23])
Shade.Coff<-sum(Shade.coffee[1:23])
Sun.Coff<-sum(Sun.coffee[1:23])
Forest.Frag<-sum(Fragment[1:23])
Fence<-sum(Live.fence[1:23])
Para<-sum(Paramo[1:23])

# Estimated density across all elevations (275-3700m) using mature forest as the land cover class
for (j in 1:56){
	elev.tot[j]<-sum(MF.eff[j,1:23])
	}

# Calculate predicted linear relationships with woody biomass (tree.eff) and the three climate covariates for each species
# Standardized scale included 91 values but predictions were later restricted to the range observed for each variable
for (i in 1:spec){
	for (j in 1:91){
		tree.eff[j,i]<-exp(alpha[i] + beta[i,9]*c[j])
		precip.eff[j,i]<-exp(alpha[i] + beta[i,10]*c[j])
		seas.eff[j,i]<-exp(alpha[i] + beta[i,11]*c[j])
		cloud.eff[j,i]<-exp(alpha[i] + beta[i,12]*c[j])
		} 
	}



