rm(list=ls())

library(coda)
library(rjags)
library(R2jags)

## note - the solution contains the full analysis for 'wrma'

## --------------------------------------------------
## load the data for all species
load(url('https://www.sfu.ca/~lmgonigl/materials-qm/data/occ-data.RData'),
     verbose=TRUE)

## --------------------------------------------------
## package data for JAGS (Boat-billed Flycatcher)
sp <- 'wrma'
my.data <- list(X=X[,,sp],
                nsite=nsite,
                forest=forest,
                nvisit=nvisit)
## remember that the inits must also be subsetted to correct species
Z.init <- Z.init[,sp]
## --------------------------------------------------

## take a look at the data
head(my.data$X, 10)
my.data$nsite
my.data$forest
my.data$nvisit

## inits
my.inits <- function() {
  list(Z=Z.init)
}

## model
my.model <- function() {

  ## PRIORS
  psi.0 ~ dnorm(0,0.001)
  p.0 ~ dnorm(0,0.001)

  psi.forest ~ dnorm(0,0.001)

  for(site in 1:nsite) {
    logit(psi[site]) <- psi.0 + psi.forest*forest[site]
  }
  logit(p) <- p.0

  
  ## LIKELIHOOD
  for(site in 1:nsite) {
    ## Occurrence
    Z[site] ~ dbern(psi[site])
    p.eff[site] <- Z[site] * p
    for(visit in 1:nvisit) {
      ## Detection 
      X[site,visit] ~ dbern(p.eff[site])
    }
  }

  ## Derived quantities
  num.occ <- sum(Z[]) ## Number of occupied sites
}

## parameters to track
my.params <- c('psi', 'psi.forest', 'p', 'num.occ')

res <- jags(data=my.data,
            inits=my.inits,
            parameters.to.save=my.params,
            model.file=my.model,
            n.iter=110000,
            n.burnin=10000,
            n.thin=10,
            n.chains=3,
            working.directory=NULL)

## to make life easier
bugs <- res$BUGSoutput
summ <- res$BUGSoutput$summary

## inspect output
columns <- c('mean','sd','2.5%','97.5%', 'Rhat')
summ[,columns]

hist(bugs$sims.array[,,'psi[1]'],
     prob=TRUE, main='', las=1,
     xlab=expression('Occupancy, '~psi), col='red')

hist(bugs$sims.array[,,'psi.forest'],
     prob=TRUE, main='', las=1,
     xlab=expression('Occupancy, '~psi), col='red')

hist(bugs$sims.array[,,'p'],
     prob=TRUE, main='', las=1,
     xlab=expression('Detectability, '~italic(p)),
     col='red')