rm(list=ls())

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

## --------------------------------------------------
## 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
my.data <- ## FILL IN HERE
Z.init <- ## FILL IN HERE
## --------------------------------------------------

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

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

  ## PRIORS
  psi ~ dunif(0,1)
  p ~ dunif(0,1)
  
  ## LIKELIHOOD
  for(site in 1:nsite) {
    ## Occurrence
    Z[site] ~ dbern(psi)
    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', '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]

## plot posteriors

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

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

hist(bugs$sims.array[,,'num.occ'],
     prob=TRUE, main='', las=1,
     xlab='Number of occupied sites',
     xlim=c(0,44), col='red')