#Code for a hierarchical occupancy model for strips estimating detection probability

setwd('C:/Analyses/chancho')
library(R2WinBUGS)

mazama <- read.table('chancho.txt', header=T)  #leo arc

dist.mean = mean(mazama$dstance)
dist.sd = sd(mazama$dstance)
under.mean = mean(mazama$und_dens)
under.sd = sd(mazama$und_dens)
canopy.mean = mean(mazama$cnp_cvr)
canopy.sd = sd(mazama$cnp_cvr)


trans <- as.numeric(factor(mazama$trans))   #variable trans is not number

data <- list(
  n=length(mazama$Site), #I specify the length so it can reiterate
  j=mazama$count,#suma de cantidad de detecciones/camara
  y=mazama$total,
  distance=((mazama$dstance-dist.mean)/dist.sd),
  canopy=((mazama$cnp_cvr-canopy.mean)/canopy.sd),
  understory=((mazama$und_dens - under.mean)/under.sd),
  trans=trans,
  ntrans=max(trans)
  )               

inits <- function(){       #cargo inits para WinBugs
  list(
  bMean=runif(1),
  b1=rnorm(1),
  aMean=runif(1),
  a1=rnorm(1),
  a2=rnorm(1),
  z=as.integer(mazama$total>0),
  transect.effect=rnorm(max(trans)),
  sd.transect=runif(1,0,5)
  )
  }

  params <- list(
    'transect.effect',
  'b0',
  'b1',
  'a0',
  'a1',
  'a2',
  'bMean',
  'occ.fs'
    )

# Define model
cat("
      model{
      
#Prior distributions for fixed effects  (B&A)

bMean ~ dunif(0,1)
b0 <- log(bMean) - log(1-bMean)       #log transform parameter
b1 ~ dnorm(0,.001)                      #set priors for b
aMean ~ dunif(0,1)
a0 <- log(aMean) - log(1-aMean)
        #set prior for a
a1 ~ dnorm(0,.001)
a2 ~ dnorm(0,.001)

#Prior distributions for random effect variance parameters
sd.transect ~ dunif(0,10)
tau.transect <- pow(sd.transect,-2)


#Prior distributions for random effects
for(r in 1:ntrans){
transect.effect[r] ~ dnorm(b0, tau.transect)}


for (i in 1:n) { #start initial loop over the n sites
	# True state model for the partially observed true state
		z[i] ~ dbern(psi[i])        # True occupancy z at site i
     	psi[i] <- 1/(1+exp(-fu[i]))
       fu[i]<- b1*distance[i]  + transect.effect[trans[i]]

		
	#observation model 
		mu.y[i] <- z[i] * p[i]
     p[i]<- 1/(1+exp(-cu[i]))
        cu[i] <- a0 +  a1*understory[i] + a2*canopy[i] 

    
   		y[i] ~ dbin(mu.y[i],j[i])     
		
	}#for loop

   
    occ.fs <- (sum(z[]))/n # Number of occupied sites 
  
  } #model

  ",fill=TRUE, file="trainingmodel.txt")

                                                        
out = bugs(data, inits, params,
model.file='trainingmodel.txt',
debug=T, n.chains=3, n.iter=400, n.burnin=70, n.thin=3, DIC=F, 
working.directory=getwd())