##################################################################################
############# Occupancy and relative abundances of introduced ungulates  ##################
############# on New Zealand’s public conservation land, 2012–2018           ##################
############# Moloney et al, R script                                                                     ###################
##################################################################################

#######################
#### Packages to load ####
#######################

library(tidyverse);library(brms);library(viridis);library(prettymapr);library(RColorBrewer);
theme_set(theme_bw())

###########################
#### User defined functions ####
###########################

HurdleFit <- function(data, model){
  modelU <- fitted(model, newdata=data, summary=F)
  modelH <- 1-fitted(model, newdata=data, dpar='hu', summary=F)
  modelC <- modelU / modelH
  modelC0 <- t(apply(modelC,2,CI95))
  modelH0 <- t(apply(modelH,2,CI95))
  modelU0 <- t(apply(modelU,2,CI95))
  data <- data.frame(data, EstU=modelU0[,1], LBU=modelU0[,2], UBU=modelU0[,3],
                     EstC=modelC0[,1], LBC=modelC0[,2], UBC=modelC0[,3],
                     EstH=modelH0[,1], LBH=modelH0[,2], UBH=modelH0[,3])
  return(data)
}


HurdleFit2 <- function(data, model){
  modelU <- fitted(model, newdata=data, summary=F)
  modelH <- 1-fitted(model, newdata=data, dpar='hu', summary=F)
  modelC <- modelU / modelH
  MmodelU <- rowMeans(modelU)
  MmodelH <- rowMeans(modelH)
  MmodelC <- rowMeans(modelC)
  data <- c(CI95(MmodelU), CI95(MmodelH), CI95(MmodelC))
  return(data)
}

###################
#### Data loading ####
###################

NZUngData <- read.csv('NZUngData.csv')
UbT <- read.csv('UngulateByTransect.csv')

######################
#### Data definitions ####
######################

# Plot: the unique plot ID
# Season: the season the survey was conducted
# FPC: faecal pellet count
# Effort: the total length of transects at this plot
# FPI: the raw Faecal Pellet Index (FPI) before transect length is accounted for.
# FPIc: the FPI corrected for transect length, (FPI/Effort*600)
# PlotProp: proportion of quadrats that had pellets
# lEffort: log of Effort
# EastingNZTM: Easting for the plot from the New Zealand Transverse Mercator 2000
# NorthingNZTM: Northing for the plot from the New Zealand Transverse Mercator 2000
# Ung_A: 1 if Sustained aerial control at this plot, 0 otherwise
# Year: the time of year that the survey season started (eg 2011-2012 season is 2011.5)
# Forest: "Yes" if plot was designated forest, "No" otherwise
# FeralGoat: 1 if the plot is within the known Feral Goat area, 0 otherwise
# Tahr: 1 if the plot is within the known Himilayan Tahr area, 0 otherwise
# SikaDeer: 1 if the plot is within the known Sika Deer area, 0 otherwise
# RedDeer: 1 if the plot is within the known Red Deer area, 0 otherwise
# Island: which island (North, South or Stewart) that the plot is on
# WARONorth: 1 if the plot is on the North Island and has WARO (wild animal recovery operation),
#            0 otherwise
# WAROSouth: 1 if the plot is on the South Island and has WARO (wild animal recovery operation),
#            0 otherwise
# ForestNorth: 1 if the plot is on the North Island and is designated forest, 0 otherwise
# ForestSouth: 1 if the plot is on the South Island and is designated forest, 0 otherwise
# SFCont: 1 if the plot is on the South Island and is designated forest and has sustained aerial
#         ungulate control at this point, 0 otherwise
# AB, DE, MN & PI: individual transects at the plot

######################
#### Ungulate model ####
######################

UngulateModel <-brm(bf(FPIc ~ t2(EastingNZTM ,NorthingNZTM, by=factor(Island)) + Island*zYear +
                         WARONorth + WAROSouth + Ung_A,
                       hu ~ t2(EastingNZTM, NorthingNZTM, by=factor(Island)) + Island*zYear +
                         WARONorth + WAROSouth + ForestNorth + ForestSouth + Ung_A+SFCont),
                    data=NZUngData, family='hurdle_lognormal', cores=8,
                    control=list(adapt_delta=0.99, max_treedepth = 20))
bayes_R2(UngulateModel)
pp_check(UngulateModel, nsamples=100)

UngulatePred <- HurdleFit(NZUngData, UngulateModel)

######################
#### Red Deer model ####
######################
RedDeerData <- NZUngData %>% filter(RedDeer==1 & Tahr==0 & FeralGoat==0  & SikaDeer==0) %>%
  mutate(Island=factor(Island, levels=c('North', 'South')))


RedDeerModel <-brm(bf(FPIc ~ t2(EastingNZTM ,NorthingNZTM, by=factor(Island)) + Island*zYear +
                        WARONorth + WAROSouth + Ung_A,
                      hu ~ t2(EastingNZTM, NorthingNZTM, by=factor(Island)) + Island*zYear +
                        WARONorth + WAROSouth + ForestNorth + ForestSouth + Ung_A + SFCont),
                   data=RedDeerData, family='hurdle_lognormal', cores=8,
                   control=list(adapt_delta=0.99, max_treedepth = 20))
bayes_R2(RedDeerModel)
pp_check(RedDeerModel, nsamples=100)

RedDeerPred <- HurdleFit(RedDeerData, RedDeerModel)

########################
#### Site plot (Figure 1) ####
########################

par(mar=c(0.5,1.2,1.5,1)+0.2)
plot(NZRegmap[-OutsideIs,], col='lightgrey')
points(UbT$EastingNZTM, UbT$NorthingNZTM, pch=15, col=UbT$Col2, cex=0.6)
legend('topleft', legend=c('Woody site', 'Non-woody site'), pch=15,
       col=c('darkgreen', 'darkblue'), bty='n', cex=2)
addscalebar(plotepsg=2193, pos='bottomright', label.cex=2)
addnortharrow('topright', scale=2)

###############################
#### Bayesian occupancy model ####
###############################

UngDataJAGS <- list(N=nrow(UbT), Obs=as.matrix(UbT[,3:6]), Year=UbT$Year,
                    Island=unclass(factor(UbT$Island,
                                          levels=c('North', 'South', 'Stewart Island'))),
                    Forest=ifelse(UbT$Forest=='Yes', 1, 0))
str(UngDataJAGS)

UngOcc.jags=function(){
  #Priors and parameter mapping
  for (r in 1:3){      # where r is region, 1 = North Is, 2 = South Is, 3 = Stewart Is
    alpha[r] ~ dunif(-10, 10) # where alpha is the intercept of the occupancy rate
    beta[r] ~ dunif(-10, 10) # where beta1 is the effect of forest on occupancy
  } #r
  beta[4] ~ dunif(-10, 10) # where beta4 is the rate of change of occupancy over time
                           # (given no site is surveyed twice)
  p ~ dbeta(1,1) # where p is the probablity that a pellets are detected on a transect at an
                 # occupied site
  #Ecological sub-model: Define the state conditions on parameters
  for(i in 1:N){
    z[i] ~ dbern(mu[i]) #true state for site i
    logit(mu[i]) <- alpha[Island[i]] + beta[Island[i]]*Forest[i] + beta[4]*Year[i]
  } #s
  #Obersvational model
  for(j in 1:N){
    for(k in 1:4){
      Obs[j,k]~dbern(z[j]*p)
    }#k
  }#j
  #Derived parameters
  p4 <- 1-pow(1-p, 4) # prob. that ungulate is detected at at least one of the four transects
  logit(psi[1,1]) <- alpha[1] # North island, nonforested
  logit(psi[1,2]) <- alpha[1] + beta[1] # North island, forest
  logit(psi[2,1]) <- alpha[2] # South island, nonforested
  logit(psi[2,2]) <- alpha[2] + beta[2] # South island, forest
  logit(psi[3,1]) <- alpha[3] # Stewart island, nonforested
  logit(psi[3,2]) <- alpha[3] + beta[3] # Stewart island, forest
  ForDif[1] <- beta[2]-beta[1]
  ForDif[2] <- beta[3]-beta[1]
  ForDif[3] <- beta[3]-beta[2]
}

init=function(){
  list(z=rep(1,nrow(ungRep)))
}

parms=c('psi','alpha','beta','p', 'p4', 'ForDif')
# MCMC settings
ni <- 1200
nt <- 1
nb <- 200
nc <- 3
sn <- 1
UngMC=jags(data=UngDataJAGS, inits=init,parameters.to.save=parms, model.file=UngOcc.jags,
           n.chains=nc, n.thin=nt*sn ,n.iter=ni*sn, n.burnin = nb*sn)

#Non-forest total
round(quantile((UngMC$BUGSoutput$sims.list$psi[, 1, 1]*56 +
                  UngMC$BUGSoutput$sims.list$psi[, 2, 1]*501 +
                  UngMC$BUGSoutput$sims.list$psi[, 3, 1]*2)/(56 + 501 + 2),
               probs=c(0.025, 0.5, 0.975)), 3)
#Forest total
round(quantile((UngMC$BUGSoutput$sims.list$psi[, 1, 2]*260 +
                  UngMC$BUGSoutput$sims.list$psi[, 2 ,2]*504+
                  UngMC$BUGSoutput$sims.list$psi[, 3, 2]*23)/(260 + 504 + 23),
               probs=c(0.025, 0.5, 0.975)), 3)
#Overall
round(quantile((UngMC$BUGSoutput$sims.list$psi[, 1, 1]*56 +
                  UngMC$BUGSoutput$sims.list$psi[, 2, 1]*501 +
                  UngMC$BUGSoutput$sims.list$psi[, 3, 1]*2 +
                  UngMC$BUGSoutput$sims.list$psi[, 1, 2]*260 +
                  UngMC$BUGSoutput$sims.list$psi[, 2, 2]*504 +
                  UngMC$BUGSoutput$sims.list$psi[, 3, 2]*23)/(56 + 501 + 2 + 260 + 504 + 23),
               probs=c(0.025, 0.5, 0.975)), 3)
#North island
round(quantile((UngMC$BUGSoutput$sims.list$psi[, 1, 1]*56 +
                  UngMC$BUGSoutput$sims.list$psi[, 1, 2]*260)/(56 + 260),
               probs=c(0.025, 0.5, 0.975)), 3)
#South island
round(quantile((UngMC$BUGSoutput$sims.list$psi[, 2, 1]*501 +
                  UngMC$BUGSoutput$sims.list$psi[, 2, 2]*504)/(501 + 504),
               probs=c(0.025, 0.5, 0.975)), 3)
#Stewart Island
round(quantile((UngMC$BUGSoutput$sims.list$psi[, 3, 1]*2 +
                  UngMC$BUGSoutput$sims.list$psi[, 3, 2]*23)/(2 + 23),
               probs=c(0.025, 0.5, 0.975)), 3)



######################################
#### Ungulate Occupancy plots (Figure 3) ####
######################################

GreenPal <- colorRampPalette(c('lightgreen', 'darkgreen'))(7)
BluePal <- colorRampPalette(c('lightblue', 'darkblue'))(7)
ptext <- c('0.0-0.1', '0.1-0.2', '0.2-0.4', '0.4-0.6', '0.6-0.8', '0.8-0.9', '0.9-1.0')
pcuts <- c(0, 0.1, 0.2, 0.4, 0.6, 0.8, 0.9, 1)

UngulatePred$pCol[which(UngulatePred$Forest=='Yes')] <-
  GreenPal[unclass(cut(UngulatePred$EstH[which(UngulatePred$Forest=='Yes')], pcuts,
                       include.lowest=T, right=F))]
UngulatePred$pCol[which(UngulatePred$Forest=='No')] <-
  BluePal[unclass(cut(UngulatePred$EstH[which(UngulatePred$Forest=='No')], pcuts,
                      include.lowest=T, right=F))]

par(mar=c(0.5,1.2,1.5,1)+0.2)
plot(NZRegmap[-OutsideIs,], col='lightgrey')
points(UngulatePred$EastingNZTM, UngulatePred$NorthingNZTM, pch=15, col=UngulatePred$pCol)
legend(x=1030000, y=6200000, rev(ptext), fill=rev(GreenPal), bty='n', title='Occupancy   ',
       cex=2)
legend(x=990000, y=6200000, legend=rep('', length(ptext)), fill=rev(BluePal), bty='n',
       title='  ', cex=2)
addscalebar(plotepsg=2193, pos='bottomright', label.cex=2)
addnortharrow('topright', scale=2)






#############################################
#### Ungulate Faecal Pellet Index plots (Figure 4) ####
#############################################

GreenPal <- colorRampPalette(c('lightgreen', 'darkgreen'))(7)
BluePal <- colorRampPalette(c('lightblue', 'darkblue'))(7)
utext <- c('0-25', '25-40', '40-50', '50-60', '60-75', '75-150', '150-300')
ucuts <- c(0, 25, 40, 50, 60, 75, 150, 310)

UngulatePred$uCol[which(UngulatePred$Forest=='Yes')] <-
  GreenPal[unclass(cut(UngulatePred$EstU[which(UngulatePred$Forest=='Yes')], ucuts,
                       include.lowest=T, right=F))]
UngulatePred$uCol[which(UngulatePred$Forest=='No')] <-
  BluePal[unclass(cut(UngulatePred$EstU[which(UngulatePred$Forest=='No')], ucuts,
                      include.lowest=T, right=F))]

par(mar=c(0.5,1.2,1.5,1)+0.2)
plot(NZRegmap[-OutsideIs,], col='lightgrey')
points(UngulatePred$EastingNZTM, UngulatePred$NorthingNZTM, pch=15, col=UngulatePred$uCol)
legend(x=990000, y=6200000, rev(utext), fill=rev(GreenPal), bty='n',
       title='Faecal Pellet Index', cex=2)
legend(x=990000, y=6200000, legend=rep('', length(utext)), fill=rev(BluePal), bty='n',
       title='  ', cex=2)
addscalebar(plotepsg=2193, pos='bottomright', label.cex=2)
addnortharrow('topright', scale=2)






#########################################################################
#### Temporal changes in the ungulate occupancy and Faecal Pellet Index (Figure 5) ####
#########################################################################

UngDataGroups <- NZUngData %>% filter(Island !='Stewart Island') %>%
  mutate(Group=fct_collapse(Season, '2011-13' = c('2011-12', '2012-13'), '2013-14'='2013-14',
                            '2014-15'='2014-15', '2015-16'='2015-16', '2016-17'='2016-17',
                            '2017-18'='2017-18'), GI = paste(Group, Island))
UngDataGroups$Ung_A=factor(UngDataGroups$Ung_A)

FirstLast <- data.frame(GI=unique(UngDataGroups$GI), MedianU=NA, LBU=NA, UBU=NA, MedianH=NA,
                        LBH=NA, UBH=NA, MedianC=NA, LBC=NA, UBC=NA) %>%
  separate(GI, c('Season', 'Island'), sep=' ', remove=FALSE)
for(i in 1:nrow(FirstLast)){
  FirstLast[i,-1:-3] <- HurdleFit2(UngDataGroups %>% filter(GI==FirstLast$GI[i]), UngulateModel)
}

FirstLastLong <- FirstLast %>% rename(Median_U='MedianU', LB_U='LBU', UB_U='UBU',
                                      Median_H='MedianH', LB_H='LBH', UB_H='UBH',
                                      Median_C='MedianC', LB_C='LBC', UB_C='UBC') %>%
  pivot_longer(Median_U:UB_C, names_to=c('.value', 'Statistic'), names_pattern='(.+)_(.+)') %>%
  mutate(Island=if_else(Island=='North', 'North Island', 'South Island'),
         Statistic=factor(Statistic, levels=c('H', 'C', 'U')))


blank_data <- data.frame(Statistic=c('H', 'H', 'C', 'C'), x=1, y=c(0, 1, 0, 300),
                         LB=c(0, 1, 0, 300), UB=c(0, 1, 0, 300), Island='North Island')
ggplot(FirstLastLong %>% filter(Statistic != 'U'), aes(x=Season, y=Median, ymin=LB, ymax=UB,
                                                       color=Island, shape=Island)) +
  geom_pointrange(position=position_dodge(0.4)) +
  geom_blank(data = blank_data, aes(x = x, y = y)) +
  facet_wrap(~Statistic, scales='free_y', nrow=2, strip.position='left',
             labeller=as_labeller(c(C='Conditional Faecal Pellet Index', H='Occupancy'))) +
  ylab(NULL) + theme(strip.background = element_blank(), strip.placement = "outside",
                     legend.title=element_blank()) +
  geom_text(data=labels, aes(x = -Inf, y = Inf, label = label, group = Statistic),
            size = 5, hjust = -0.5, vjust = 1.4, inherit.aes = FALSE)

######################################
#### Red Deer Occupancy plots (Figure 6) ####
######################################

GreenPal <- colorRampPalette(c('lightgreen', 'darkgreen'))(7)
BluePal <- colorRampPalette(c('lightblue', 'darkblue'))(7)
pRtext <- c('0.0-0.2', '0.2-0.4', '0.4-0.6', '0.6-0.7', '0.7-0.8', '0.8-0.9', '0.9-1.0')
pRcuts <- c(0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1)

RedDeerPred$pCol[which(RedDeerPred$Forest=='Yes')] <-
  GreenPal[unclass(cut(RedDeerPred$EstH[which(RedDeerPred$Forest=='Yes')], pRcuts,
                       include.lowest=T, right=F))]
RedDeerPred$pCol[which(RedDeerPred$Forest=='No')] <-
  BluePal[unclass(cut(RedDeerPred$EstH[which(RedDeerPred$Forest=='No')], pRcuts,
                      include.lowest=T, right=F))]

par(mar=c(0.5,1.2,1.5,1)+0.2)
plot(NZRegmap[-OutsideIs,], col='lightgrey')
points(RedDeerPred$EastingNZTM, RedDeerPred$NorthingNZTM, pch=15, col=RedDeerPred$pCol)
legend(x=1030000, y=6200000, rev(pRtext), fill=rev(GreenPal), bty='n', title='Occupancy   ',
       cex=2)
legend(x=990000, y=6200000, legend=rep('', length(pRtext)), fill=rev(BluePal), bty='n',
       title='  ', cex=2)
addscalebar(plotepsg=2193, pos='bottomright', label.cex=2)
addnortharrow('topright', scale=2)

#############################################
#### Red Deer Faecal Pellet Index plots (Figure 7) ####
#############################################

GreenPal <- colorRampPalette(c('lightgreen', 'darkgreen'))(7)
BluePal <- colorRampPalette(c('lightblue', 'darkblue'))(7)
utext <- c('0-25', '25-40', '40-50', '50-60', '60-75', '75-150', '150-300')
ucuts <- c(0, 25, 40, 50, 60, 75, 150, 310)

RedDeerPred$uCol[which(RedDeerPred$Forest=='Yes')] <-
  GreenPal[unclass(cut(RedDeerPred$EstU[which(RedDeerPred$Forest=='Yes')], ucuts,
                       include.lowest=T, right=F))]
RedDeerPred$uCol[which(RedDeerPred$Forest=='No')] <-
  BluePal[unclass(cut(RedDeerPred$EstU[which(RedDeerPred$Forest=='No')], ucuts,
                      include.lowest=T, right=F))]

par(mar=c(0.5,1.2,1.5,1)+0.2)
plot(NZRegmap[-OutsideIs,], col='lightgrey')
points(RedDeerPred$EastingNZTM, RedDeerPred$NorthingNZTM, pch=15, col=RedDeerPred$uCol)
legend(x=990000, y=6200000, rev(utext), fill=rev(GreenPal), bty='n',
       title='Faecal Pellet Index', cex=2)
legend(x=990000, y=6200000, legend=rep('', length(utext)), fill=rev(BluePal), bty='n',
       title='  ', cex=2)
addscalebar(plotepsg=2193, pos='bottomright', label.cex=2)
addnortharrow('topright', scale=2)



#######################################################
#### Observed ungulate FPI plots (Suppliemtnary material S7) ####
#######################################################

GreenPal <- colorRampPalette(c('lightgreen', 'darkgreen'))(7)
BluePal <- colorRampPalette(c('lightblue', 'darkblue'))(7)
otext <- c('0', '1-10', '10-25', '25-50', '50-100', '100-200', '>200')
ocuts <- c(0, 1, 10, 25, 50, 100, 200, 1200)

RawClass <- data.frame(EastingNZTM=NZUngData$EastingNZTM, NorthingNZTM=NZUngData$NorthingNZTM,
                       FPIc=NZUngData$FPIc, Forest=NZUngData$Forest, Col=NA)

RawClass$Col[which(RawClass$Forest=='Yes')] <-
  GreenPal[unclass(cut(RawClass$FPIc[which(RawClass$Forest=='Yes')], ucuts, include.lowest=T,
                       right=F))]
RawClass$Col[which(RawClass$Forest=='No')] <-
  BluePal[unclass(cut(RawClass$FPIc[which(RawClass$Forest=='No')], ucuts, include.lowest=T,
                      right=F))]

plot(NZRegmap[-OutsideIs,], col='lightgrey')
points(RawClass$EastingNZTM, RawClass$NorthingNZTM, pch=15, col=RawClass$Col)
legend(x=990000, y=6200000, rev(otext), fill=rev(GreenPal), bty='n',
       title='Faecal Pellet Index', cex=2)
legend(x=990000, y=6200000, legend=rep('', length(otext)), fill=rev(BluePal), bty='n',
       title='  ', cex=2)
addscalebar(plotepsg=2193, pos='bottomright', label.cex=2)
addnortharrow('topright', scale=2)



##################################################################################
#### Observed Faecal Pellet Index by woody and non-woody sites (Suppliemtnary material S8) ####
##################################################################################

ggplot(NZUngData %>%
         mutate(Wood=factor(ifelse(Forest=='Yes', 'Woody', 'Non-woody'),
                            levels=c('Woody', 'Non-woody'))), aes(FPIc)) +
  geom_histogram(binwidth=10, aes(y=..density..)) + facet_wrap(Wood~., nrow=2) +
  scale_y_continuous(name='Density', limits=c(0, 0.06), expand=c(0,0)) +
  scale_x_continuous(name='Mean observed FPI', expand=c(0,0), limits=c(-10, 1200)) +
  theme(plot.margin=unit(c(0.5,1,0,0.5), 'cm'))
