index-jags.Rmd

---
title: "Distribution of Mediterranean monk seals in Greece using occupancy models"
author: "Olivier Gimenez"
date: "September 2021, updated May 2022"
output:
  html_document:
    code_folding: show
    df_print: paged
    highlight: tango
    number_sections: yes
    theme: united
    toc: yes
    toc_depth: 1
    keep_md: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, 
                      dpi = 300, 
                      message = FALSE, 
                      warning = FALSE)
```

# Introduction

This is an analysis of the RINT presence-only data gathered by the [MOm NGO](http://www.mom.gr/homepage.asp?ITMID=101&LANG=EN) through a citizen-science program. The objective is to map the distribution of monk seals by fitting occupancy models to these data. 

Load all the packages we will need. 
```{r}
library(tidyverse)
library(lubridate)
theme_set(theme_light())
library(sf)
library(viridis)
library(gridExtra)
library(transformr)
library(here)
```

# Exploratory data analyses

Let's load and explore the data to get a feeling of the info we have (and we do not have). 
```{r}
load(here::here("data","monseal.RData"))
```

First, the number of sightings over time.
```{r}
monkseal %>% 
  ggplot() +
  aes(x = year) + 
  geom_bar()
```

The sightings are done along the year (all years are pooled together). August gets maximum sightings.
```{r}
monkseal %>% 
  ggplot() +
  aes(x = factor(month)) + 
  geom_bar() + 
  xlab('month')
```

Regarding observers, sightings are mostly done by local people, then to a lesser extent tourists, spear gun fishers, professional fishermen and sailmen, port police and a few others. 
```{r}
monkseal %>% 
  filter(!is.na(observer)) %>%
  count(observer) %>%
  ggplot() +
  aes(x = n, y = fct_reorder(observer, n)) %>%
  geom_col() +
  labs(x = NULL, y = NULL)
```

Last, regarding where the seals were when they were spotted, we see that the sightings are mostly at sea (h), from human settlement (i) and on beach (f). 
```{r}
monkseal %>% 
  filter(!is.na(where)) %>%
  count(where) %>%
  ggplot() +
  aes(x = n, y = fct_reorder(where, n)) %>%
  geom_col() +
  labs(x = NULL, y = NULL) 
```

# Data mapping

Before diving deep into the statistical analyses, let's put the sightings on a map. 

First, let's get a map of Greece.
```{r}
greece <- st_read(here::here("shapefiles","greece.shp"))
coastlines <- st_read(here::here("shapefiles","coastlines.shp"))
ggplot() + 
  geom_sf(data = greece) + 
  geom_sf(data = coastlines, color = 'red') + 
  labs(title = 'Map of Greece and coastlines', x = NULL, y = NULL)
```

And our grid (including a 20-km buffer around coastlines).
```{r}
grid <- st_read(here::here("shapefiles","grid.shp"))
ggplot() + 
  geom_sf(data = greece) + 
  geom_sf(data = grid, lwd = 0.1) + 
  labs(title = 'Gridded map of Greece', x = NULL, y = NULL)
```

Note that a few sites are small, more precisely `r round(mean(as.numeric(grid$areakm2)<10))` percent of the `r nrow(grid)` (i.e. `r sum(as.numeric(grid$areakm2)<10)` sites) are less than $10km^2$.

Get all sightings and map them:
```{r}
counts <- monkseal %>%
  group_by(utm) %>%
  count()
names(counts) <- c('UTM_NAME','val')
grid_wcounts <- right_join(grid, counts)
classes <- classInt::classIntervals(grid_wcounts$val, n = 5, style = "jenks")
bwr <- grDevices::colorRampPalette(rev(c('red', 'orange', 'yellow', 'cyan', 'blue')))
grid_wcounts <- grid_wcounts %>%
  mutate(valdiscrete = cut(val, classes$brks, include.lowest = T))
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = grid, colour = "grey50", fill = "white", lwd = 0.1) + 
  geom_sf(data = grid_wcounts, lwd = 0.1, aes(fill = valdiscrete)) + 
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
  scale_fill_manual(
    name = 'counts',
    values = bwr(5),
    guide = guide_legend(reverse = TRUE)) + 
  labs(title = '', x = NULL, y = NULL)
ggsave(here::here("figures","map1.png"), dpi = 600)
```

# Occupancy analysis

## Data formating

Prepare the data.
```{r}
grid <- grid %>% 
  filter(as.numeric(areakm2)>10) # filter out small sites
mm <- c('Mar','Apr','May','Jun','Jul','Aug','Sep','Oct') # extract month with data
indyear <- c(2000:2007, 2013:2020) # define two periods
ids <- unique(grid$UTM_NAME)
mom_detections <- list()
for (k in 1:length(indyear)){
  mom_detections[[k]] <- matrix(0,nrow=length(ids),ncol=length(mm))
  data <- filter(monkseal,year==indyear[k])
  ind <- 1
  for (i in mm){
    temp <- filter(data, month==i)
    utm_temp <- unique(temp$utm)
    for (j in utm_temp){
      # 1 if a sighting was made in square j in month i
      mom_detections[[k]][ids==j,ind] <- 1
    }
    ind <- ind + 1
  }
}
mom_det <- do.call(cbind, mom_detections) # bind data from all year in columns
# convert the occupancy dataset in a 3D array:
y <- list()
ind <- 0
for (i in 1:length(indyear)){
	mask <- (ind + i):(ind + i + length(mm) - 1)
	y[[i]] <- mom_det[,mask]
	ind <- ind + length(mm) - 1
}

# convert list into array (https://stackoverflow.com/questions/37433509/convert-list-to-a-matrix-or-array)
y <- array(unlist(y), dim = c(nrow(y[[1]]), ncol(y[[1]]), length(y)))
dim(y)
```

Now we reformat the data as follows: first, for each year, we determine whether the species was ever detected or not; then for each period, we count the number of times the species was detected:
```{r}
new_y <- NULL
for (i in 1:dim(y)[3]){ # loop over years
  new_y <- cbind(new_y,apply(y[,1:8,i],1,sum))
}
new_y <- (new_y > 0) 
y1 <- apply(new_y[,1:8],1,sum)
y2 <- apply(new_y[,9:16],1,sum)
y <- cbind(y1,y2)
dim(y)
range(y)
```

## Inference

We now specify our model in BUGS language:
```{r}
sink(here::here("codes","dynoccct.jags"))
cat("
model{

# Specify priors
mupsi ~ dnorm(-5,1)
for(i in 1:nsite){
cloglog(psi[i]) <- mupsi + offset[i]
}

sdgam ~ dunif(0,3)
taugam <- pow(sdgam,-1)
sdeps ~ dunif(0,3)
taueps <- pow(sdeps,-1)
sdp ~ dunif(0,10)
taup <- pow(sdp,-2)

for(i in 1:nyear){
    mup.prob[i] ~ dunif(0,1)
    logit(mup[i]) <- mup.prob[i]
}

for(i in 1:nsite){
    lp[i] ~ dnorm(0,taup)
        for(t in 1:nyear){
            logit(p[i,t]) <- mup[t] + lp[i]
      }
}

mueps ~ dnorm(-5,1)
mugam ~ dnorm(-5,1)

# Ecological submodel: Define state conditional on parameters
for(i in 1:nsite){
   z[i,1] ~ dbern(psi[i])
   leps[i] ~ dnorm(0,1)
   lgam[i] ~ dnorm(0,1)
   cloglog(gamma[i]) <- mugam + taugam * lgam[i] + offset[i]
   cloglog(eps[i]) <- mueps + taueps * leps[i] + offset[i]
   muZ[i,2] <- z[i,1] * (1-eps[i]) + (1-z[i,1]) * gamma[i]
   z[i,2] ~ dbern(muZ[i,2])
}

# Observation model
for (i in 1:nsite){
		for (k in 1:nyear){
			y[i,k] ~ dbin(muy[i,k],8)
			muy[i,k] <- z[i,k]*p[i,k]
		} #k
} #i

}
", fill = TRUE)
sink()
```

Specify data, initial values, parameters to be monitored and various MCMC details:
```{r}
# data
offset <- log(as.numeric(grid$areakm2))
win.data <- list(y = y, nsite = dim(y)[1], nyear = dim(y)[2], offset = offset)

# initial values
zst <- cbind(as.numeric(y1 > 0), as.numeric(y2 > 0)) # observed occurrence as inits for z 
inits <- function() {list(z = zst, 
                          sdp = runif(1,1,9), 
                          sdgam = runif(1,0.1,2.5), 
                          sdeps = runif(1,0.1,2.5))}

# parameters monitored
params <- c("psi","mupsi","sdgam","sdeps","sdp","mup","lp","mugam","lgam","mueps","leps","z")

# MCMC settings
ni <- 15000
nb <- 5000
nc <- 2
ni <- 250
nb <- 50
nc <- 2
```

Run JAGS from R:
```{r eval=FALSE}
library(R2jags)
ptm <- proc.time()
out <- jags(data = win.data, 
            inits = inits, 
            parameters.to.save = params, 
            model.file = here::here("codes","dynoccct.jags"), 
            n.chains = nc, 
            n.iter = ni, 
            n.burnin = nb,
            n.thin = 1)
x <- proc.time() -  ptm
save(out, x, file = here::here("models","monksealscloglog.RData"))
```

The code above takes some time to run. I run it once, saved the results and use them from here:
```{r}
load(here::here("models","monksealscloglog.RData"))
```

Check convergence:
```{r}
library(MCMCvis)
MCMCtrace(object = out,
          pdf = FALSE,
          ind = TRUE,
          Rhat = TRUE, # add Rhat
          n.eff = TRUE, # add eff sample size
          params = c("sdgam","sdeps","sdp",
                     "mup","mugam","mueps"))
```

OK, random effects on extinction/col are not needed, also we add CAR on occupancy. 

We now specify our model in BUGS language:
```{r}
sink(here::here("codes","dynocc.jags"))
cat("
model{

# Specify priors for occupancy
mupsi ~ dnorm(-5,1)
tau ~ dunif(0,10)#dgamma(2,2) 
alpha ~ dunif(0,1)#1	    	

# Occupancy
for(i in 1:nsite){
cloglog(psi[i]) <- mupsi + phi[i] + offset[i]
mu[i] <- 0
}
# CAR
inv.c[1:nsite,1:nsite] <- tau * (D[1:nsite,1:nsite] - alpha * W[1:nsite,1:nsite])
phi[1:nsite] ~ dmnorm(mu[1:nsite],inv.c[1:nsite,1:nsite])

# Specify priors for detection
sdp ~ dunif(0,10)
taup <- pow(sdp,-2)
for(i in 1:nyear){
    mup.prob[i] ~ dunif(0,1)
    logit(mup[i]) <- mup.prob[i]
}
# Detection
for(i in 1:nsite){
    lp[i] ~ dnorm(0,taup)
        for(t in 1:nyear){
            logit(p[i,t]) <- mup[t] + lp[i]
      }
}

# prior
mueps ~ dnorm(-5,1)
mugam ~ dnorm(-5,1)

# Dynamic state model
for(i in 1:nsite){
   z[i,1] ~ dbern(psi[i])
   leps[i] ~ dnorm(0,1)
   lgam[i] ~ dnorm(0,1)
   cloglog(gamma[i]) <- mugam + offset[i]
   cloglog(eps[i]) <- mueps + offset[i]
   muZ[i,2] <- z[i,1] * (1-eps[i]) + (1-z[i,1]) * gamma[i]
   z[i,2] ~ dbern(muZ[i,2])
}

# Observation model
for (i in 1:nsite){
		for (k in 1:nyear){
			y[i,k] ~ dbin(muy[i,k],8)
			muy[i,k] <- z[i,k]*p[i,k]
		} #k
} #i

}
", fill = TRUE)
sink()
```

Specify data, initial values, parameters to be monitored and various MCMC details:
```{r}
# spatial ingredients
neighbours <- spdep::poly2nb(pl = grid)
#plot(grid,border="grey")
#plot(neighbours,spdep::coordina(grid),pch=16,add=TRUE)
# proximity matrix
weights <- spdep::nb2mat(neighbours,style="B")

# alpha = 1/max(eigen(weights)$values)
# tau2 = 1
# P = (diag(rowSums(weights)) - alpha * weights)/tau2
## Is P valid?
# library(MASS)
# Sigma = chol2inv(P)
# mvrnorm(n=1,mu=rep(0,dim(Sigma)[1]),Sigma=Sigma) #Useful for simulations later on

# data
offset <- log(as.numeric(grid$areakm2))
win.data <- list(y = y, 
                 nsite = dim(y)[1], 
                 nyear = dim(y)[2], 
                 offset = offset,
                 D = diag(rowSums(weights)),
                 W = weights)

# initial values
zst <- cbind(as.numeric(y1 > 0), as.numeric(y2 > 0)) # observed occurrence as inits for z 
inits <- function() {list(z = zst, 
                          sdp = runif(1,1,9))}

# parameters monitored
params <- c("psi","mupsi","sdp","mup","lp","mugam","lgam","mueps","leps","z","alpha","tau")

# MCMC settings
ni <- 1500
nb <- 500
nc <- 2
```

Run JAGS from R:
```{r eval=FALSE}
library(R2jags)
ptm <- proc.time()
out <- jags(data = win.data, 
            inits = inits, 
            parameters.to.save = params, 
            model.file = here::here("codes","dynocc.jags"), 
            n.chains = nc, 
            n.iter = ni, 
            n.burnin = nb,
            n.thin = 1)
x <- proc.time() -  ptm
save(out, x, file = here::here("models","monksealscloglogspatial.RData"))
```

The code above takes some time to run. I run it once, saved the results and use them from here:
```{r}
load(here::here("models","monksealscloglogspatial.RData"))
```

Check convergence:
```{r}
library(MCMCvis)
MCMCtrace(object = out,
          pdf = FALSE,
          ind = TRUE,
          Rhat = TRUE, # add Rhat
          n.eff = TRUE, # add eff sample size
          params = c("sdp","mup","mugam","mueps"))
```


Print results:
```{r, eval = FALSE}
out$BUGSoutput$sims.matrix %>%
  as_tibble() %>%
  pivot_longer(cols = everything(),  values_to = "value", names_to = "parameter") %>%
  filter(str_detect(parameter, "mu") | str_detect(parameter, "sd")) %>%
  group_by(parameter) %>%
  summarize(median = median(value),
            lci = quantile(value, probs = 2.5/100),
            uci = quantile(value, probs = 97.5/100))
```

## Post-process results

Check out the number of occupied UTM squares by the species over time. First work out the naive occupancy, that is the number of sites that were observed occupied over the two periods:
```{r}
naiveocc <- rep(NA,2) 
for (i in 1:2){
  naiveocc[i] <- sum(y[,i]>0)
}
naiveocc
```

Now the estimated occupancy:
```{r}
zz <- out$BUGSoutput$sims.list$z
dim(zz) # 2000 x 2850 x 2 (nbMCMC x nSquares x nyears)
nbpixocc <- apply(zz,c(1,3),sum) # 2000 x 2
meannbpixocc <- apply(nbpixocc,2,median)
sdnbpixocc <- apply(nbpixocc,2,sd)
qulnbpixocc <- apply(nbpixocc,2,quantile,probs=2.5/100)
quunbpixocc <- apply(nbpixocc,2,quantile,probs=97.5/100)
```

Put everything together:
```{r}
occupancy <- data.frame(period = c(1,2),
  naiveocc = naiveocc, 
  medianocc = meannbpixocc,
  qulnbpixocc = qulnbpixocc,
  quunbpixocc = quunbpixocc)
occupancy2 <- data.frame(period = rep(c(1,2),2),
  occ = c(occupancy$naiveocc, occupancy$medianocc), 
  quantlower = c(rep(NA,2),occupancy$qulnbpixocc),
  quantupper = c(rep(NA,2), occupancy$quunbpixocc),
  Occupancy = c(rep('naive',2),rep('estimated',2)))
occupancy2
```

Now plot:
```{r}
occupancy2 %>% 
  ggplot() + 
  aes(x = period, y = occ, group = Occupancy, linetype = Occupancy) + 
  geom_ribbon(data = occupancy2,
              aes(ymin = quantlower, ymax = quantupper),
              alpha = 0.1, 
              show.legend = FALSE) + 
  geom_line() +
  scale_linetype_manual(values = c("solid", "dashed")) +  
  geom_point() + 
  labs(x = 'Year', y = 'Number of occupied sites')
```

Now display local extinction, colonization and species detection probabilities estimates with credible intervals:
```{r}
#install.packages("remotes")
#remotes::install_github("mikemeredith/MMmisc")
epsmean <- MMmisc::invcloglog(out$BUGSoutput$mean$mueps + mean(offset))
epsql <- MMmisc::invcloglog(out$BUGSoutput$mean$mueps + mean(offset) - 2*out$BUGSoutput$sd$mueps)
epsqu <- MMmisc::invcloglog(out$BUGSoutput$mean$mueps + mean(offset) + 2*out$BUGSoutput$sd$mueps)

gammean <- MMmisc::invcloglog(out$BUGSoutput$mean$mugam + mean(offset))
gamql <- MMmisc::invcloglog(out$BUGSoutput$mean$mugam + mean(offset) - 2*out$BUGSoutput$sd$mugam)
gamqu <- MMmisc::invcloglog(out$BUGSoutput$mean$mugam + mean(offset) + 2*out$BUGSoutput$sd$mugam)

pmean <- plogis(out$BUGSoutput$mean$mup)
pqu <- plogis(out$BUGSoutput$mean$mup + 2*out$BUGSoutput$sd$mup)
pql <- plogis(out$BUGSoutput$mean$mup - 2*out$BUGSoutput$sd$mup)

eps <- data.frame(
  param = epsmean, 
  qlo = epsql,
  qup = epsqu)

gam <- data.frame(
  param = gammean, 
  qlo = gamql,
  qup = gamqu)

det <- data.frame(period = c(1,2),
  param = pmean,
  qlo = pql,
  qup = pqu)

eps
gam
det
```

Let's map the estimated distribution of monk seals. First, compute realized occupancy per site and per period, and put altogether in a big table:
```{r}
zz <- out$BUGSoutput$sims.list$z # 2000 x 2850 x 2 (nbMCMC x nSquares x nyears)
meanz <- apply(zz,c(2,3),mean)
dim(meanz)
toplot <- data.frame(
  meanz = c(meanz[,c(1,2)]), 
  UTM_NAME = rep(unique(grid$UTM_NAME),2), 
  yearr = c(rep('2000-2007',nrow(meanz)), rep('2013-2020',nrow(meanz)))) 
ring_occ <- left_join(grid,toplot)
```

```{r}
nrow(meanz)
sum(meanz[,2] == meanz[,1])
sum(meanz[,2] > meanz[,1])
sum(meanz[,2] < meanz[,1])
```

Then, plot two maps, one for each period:
```{r}
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = grid, colour = "grey50", fill = "white", lwd = 0.1) + 
  geom_sf(data = ring_occ, lwd = 0.1, aes(fill = as.numeric(meanz))) + 
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
  scale_fill_viridis(
    name = 'Pr(occupancy)', 
    direction = -1,
    alpha = 0.7) + 
  labs(title = '') +
  xlab("") + ylab("") + 
  facet_wrap(~ yearr, ncol = 2) + 
  theme(legend.position = 'bottom', 
        legend.title = element_text(size=8), 
        legend.text = element_text(size=8),
        legend.key.size = unit(0.5, "cm"))
ggsave(here::here("figures","map2.png"), dpi = 600)
```

We'd also like to have two maps of the estimated occupied sites for the two periods. To do so, we say that all squares with prob of occupancy greater than $25\%$ are occupied. This gives the following maps: 
```{r}
zz <- out$BUGSoutput$sims.list$z # 2000 x 2850 x 2 (nbMCMC x nSquares x nyears)
sumz <- apply(zz,c(2,3),sum)
meanz <- apply(zz,c(2,3),sum)/dim(zz)[1]*100
z_occupied <- (meanz > 25) # 
toplot <- data.frame(
  z_occupied = c(z_occupied[,c(1,2)]), 
  UTM_NAME = rep(unique(grid$UTM_NAME),2), 
  yearr = c(rep('2000-2007',nrow(z_occupied)), rep('2013-2020',nrow(z_occupied)))) 
ring_occ <- left_join(grid,toplot) %>%
  mutate(z_occupied = if_else(z_occupied == TRUE, 'used', 'unused by the species'))
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = ring_occ, lwd = 0.1, aes(fill = z_occupied)) + 
    geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
    scale_fill_manual(values = c('white','blue'),
                    name = "Site is") + 
  xlab("") + ylab("") +
  facet_wrap(~ yearr, ncol = 2) + 
  theme(legend.position = "none")
ggsave(here::here("figures","map3.png"), dpi = 600)
```

```{r echo = FALSE}
left_join(grid,toplot) %>%
  filter(z_occupied == TRUE) %>%
  group_by(yearr) %>%
  summarize(total_area = sum(areakm2)) 
```


We'd also like to map the difference in occupancy probabilities.
```{r}
delta <- as.matrix(meanz[,2]/100 - meanz[,1]/100)
ring_occ2 <- ring_occ %>%
  filter(yearr == '2013-2020')
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = grid, colour = "grey50", fill = "white", lwd = 0.1) + 
  geom_sf(data = ring_occ2, lwd = 0.1, aes(fill = delta)) + 
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(15,"RdYlGn"),
                       name = '') + 
  labs(title = '') +
  xlab("") + ylab("") + 
  theme(legend.position = 'bottom', 
        legend.title = element_text(size=8), 
        legend.text = element_text(size=8),
        legend.key.size = unit(0.5, "cm"))
ggsave(here::here("figures","map4.png"), dpi = 600)
```

# Intersect with protected areas

```{r}
pa <- st_read(here::here("shapefiles","N2000_v31_marine_WGS84.shp"))
pa <- st_transform(pa,2100)
sum(rapply(st_geometry(pa), nrow))
pa <- st_simplify(pa,dTolerance=50)
sum(rapply(st_geometry(pa), nrow))
pa %>% 
  mutate(area = st_area(.),
         areakm2 = units::set_units(area, km^2)) %>%
  summarize(total = sum(areakm2))
```

Let's have a look to the map of marine protected areas:
```{r}
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "grey50") + 
  geom_sf(data = pa, colour = "red", fill = "transparent", lwd = 0.1) 
ggsave(here::here("figures","mpa.png"), dpi = 600)
```

Look at where the high probability of presence on monk seals pups coincides with marine protected areas, for the period 2013-2020.
```{r}
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = ring_occ2, lwd = 0.1, aes(fill = z_occupied)) +
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
    scale_fill_manual(values = c('white','blue'),
                    name = "Site is") + 
  geom_sf(data = pa, colour = "red", fill = "transparent", lwd = 0.1) +
  xlab("") + ylab("") +
  theme(legend.position = "none")
ggsave(here::here("figures","mpamonkseals.png"), dpi = 600)
```

Out of the total number of cells with high probability of occurrence of monk seals, this number intersects with marine protected areas.
```{r}
mask <- ring_occ %>%
  filter(z_occupied == 'used') %>%
  filter(yearr == '2013-2020') %>%
  st_intersects(pa) %>%
  lengths > 0
sum(mask)
ring_occ2$map <- ring_occ2$z_occupied
mask2 <- (ring_occ2$map == 'used')
ring_occ2$map[mask2][mask] <- 'used, overlap w/ mpa'
ring_occ2$map[mask2][!mask] <- 'used, no overlap w/ mpa'
```

On a map, this gives.
```{r}
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = ring_occ2, lwd = 0.1, aes(fill = map)) +
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
    scale_fill_manual(values = c('white','blue','red'),
                    name = "Site is") + 
  geom_sf(data = pa, colour = "red", fill = "transparent", lwd = 0.1) +
  xlab("") + ylab("") +
  theme(legend.position = "bottom")
ggsave(here::here("figures","mpamonksealsdetails.png"), dpi = 600)
```

# Focus on pups

We do the same analyses with pups.

```{r}
monkseal <- monkseal %>%
  filter(stage=='PUP0' | stage=='PUP1' | stage=='PUP2' | stage=='PUP3' | stage=='PUP4')
```

The number of sightings.
```{r}
monkseal %>% 
  ggplot() +
  aes(x = year) + 
  geom_bar()
```

The sightings are done along the year (all years are pooled together). October gets maximum sightings.
```{r}
monkseal %>% 
  ggplot() +
  aes(x = factor(month)) + 
  geom_bar() + 
  labs(x = NULL)
```

Regarding the observers, we see that the sightings are mostly done by local people, then to a lesser extent tourists, and a few others. 
```{r}
monkseal %>% 
  filter(!is.na(observer)) %>%
  count(observer) %>%
  ggplot() +
  aes(x = n, y = fct_reorder(observer, n)) %>%
  geom_col() +
  labs(y = NULL, x = NULL)
```

Last, regarding where the seals were when they were spotted, we see that the sightings are mostly on beach (f). 
```{r}
monkseal %>% 
  filter(!is.na(where)) %>%
  count(where) %>%
  ggplot() +
  aes(x = n, y = fct_reorder(where, n)) %>%
  geom_col() +
  labs(y = NULL, x = NULL)
```

Map counts.
```{r}
counts <- monkseal %>%
  group_by(utm) %>%
  count()
names(counts) <- c('UTM_NAME','val')
grid_wcounts <- right_join(grid, counts)
classes <- classInt::classIntervals(grid_wcounts$val, n = 3, style = "jenks")
bwr <- grDevices::colorRampPalette(rev(c('red', 'orange', 'yellow', 'cyan', 'blue')))
grid_wcounts <- grid_wcounts %>%
  mutate(valdiscrete = cut(val, classes$brks, include.lowest = T))
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = grid, colour = "grey50", fill = "white", lwd = 0.1) + 
  geom_sf(data = grid_wcounts, lwd = 0.1, aes(fill = valdiscrete)) + 
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
  scale_fill_manual(
    name = 'counts',
    values = bwr(5),
    guide = guide_legend(reverse = TRUE)) + 
  labs(title = '', x = NULL, y = NULL)
ggsave(here::here("figures","map1pups.png"),dpi=600)
```

Format data for occupancy analyses.
```{r}
mm <- c('Sep','Oct','Nov','Dec','Jan') # extract data
indyear <- c(2000:2007,2013:2020)
ids <- unique(grid$UTM_NAME)
mom_detections <- list()
for (k in 1:length(indyear)){
  mom_detections[[k]] <- matrix(0,nrow=length(ids),ncol=length(mm))
  data <- filter(monkseal,year==indyear[k])
  ind <- 1
  for (i in mm){
    temp <- filter(data, month==i)
    utm_temp <- unique(temp$utm)
    for (j in utm_temp){
      # 1 if a sighting was made in square j in month i
      mom_detections[[k]][ids==j,ind] <- 1
    }
    ind <- ind + 1
  }
}
mom_det <- do.call(cbind, mom_detections) # bind data from all year in columns
# convert the occupancy dataset in a 3D array:
y <- list()
ind <- 0
for (i in 1:length(indyear)){
	mask <- (ind + i):(ind + i + length(mm) - 1)
	y[[i]] <- mom_det[,mask]
	ind <- ind + length(mm) - 1
}

# convert list into array (https://stackoverflow.com/questions/37433509/convert-list-to-a-matrix-or-array)
y <- array(unlist(y), dim = c(nrow(y[[1]]), ncol(y[[1]]), length(y)))
dim(y)

new_y <- NULL
for (i in 1:dim(y)[3]){ # loop over years
  new_y <- cbind(new_y,apply(y[,1:5,i],1,sum))
}
new_y <- (new_y > 0) 
y1 <- apply(new_y[,1:8],1,sum)
y2 <- apply(new_y[,9:16],1,sum)
y <- cbind(y1,y2)
dim(y)
range(y)
```

We specify our model in BUGS language.
```{r}
sink(here::here("codes","dynoccctpups.jags"))
cat("
model{

# Specify priors
mupsi ~ dnorm(-5,1)
for(i in 1:nsite){
cloglog(psi[i]) <- mupsi + offset[i]
}

sdgam ~ dunif(0,10)
taugam <- pow(sdgam,-2)
sdeps ~ dunif(0,10)
taueps <- pow(sdeps,-2)
sdp ~ dunif(0,10)
taup <- pow(sdp,-2)

for(i in 1:nyear){
    mup.prob[i] ~ dunif(0,1)
    logit(mup[i]) <- mup.prob[i]
}

for(i in 1:nsite){
    lp[i] ~ dnorm(0,taup)
        for(t in 1:nyear){
            logit(p[i,t]) <- mup[t] + lp[i]
      }
}

mueps ~ dnorm(-5,1)
mugam ~ dnorm(-5,1)

# Ecological submodel: Define state conditional on parameters
for(i in 1:nsite){
   z[i,1] ~ dbern(psi[i])
   leps[i] ~ dnorm(0,taueps)
   lgam[i] ~ dnorm(0,taugam)
   cloglog(gamma[i]) <- mugam + lgam[i] + offset[i]
   cloglog(eps[i]) <- mueps + leps[i] + offset[i]
   muZ[i,2] <- z[i,1] * (1-eps[i]) + (1-z[i,1]) * gamma[i]
   z[i,2] ~ dbern(muZ[i,2])
}

# Observation model
for (i in 1:nsite){
		for (k in 1:nyear){
			y[i,k] ~ dbin(muy[i,k],5)
			muy[i,k] <- z[i,k]*p[i,k]
		} #k
} #i

}
", fill = TRUE)
sink()
```

Specify data, initial values, parameters to be monitored and various MCMC details:
```{r}
# data
offset = log(as.numeric(grid$areakm2))
win.data <- list(y = y, nsite = dim(y)[1], nyear = dim(y)[2], offset = offset)

# initial values
zst <- cbind(as.numeric(y1 > 0), as.numeric(y2 > 0)) # observed occurrence as inits for z 
inits <- function() {list(z = zst, sdp = runif(1,1,9), sdgam = runif(1,1,9), sdeps = runif(1,1,9))}

# parameters monitored
params <- c("psi","sdgam","sdeps","sdp","mup","lp","mugam","lgam","mueps","leps","z")

# MCMC settings
ni <- 25000
nb <- 10000
nc <- 2
```

Run JAGS from R:
```{r eval=FALSE}
library(R2jags)
ptm <- proc.time()
out <- jags(data = win.data, inits = inits, parameters.to.save = params, model.file = here::here("codes","dynoccctpups.jags"), n.chains = nc, n.iter = ni, n.burnin = nb)
x <- proc.time() -  ptm
save(out, x, file = here::here("models","monksealscloglogpups.RData"))
```

The code above takes some time to run. I run it once, saved the results and use them from here:
```{r}
load(here::here("models","monksealscloglogpups.RData"))
```

Print results:
```{r}
out$BUGSoutput$sims.matrix %>%
  as_tibble() %>%
  pivot_longer(cols = everything(),  values_to = "value", names_to = "parameter") %>%
  filter(str_detect(parameter, "mu") | str_detect(parameter, "sd")) %>%
  group_by(parameter) %>%
  summarize(median = median(value),
            lci = quantile(value, probs = 2.5/100),
            uci = quantile(value, probs = 97.5/100))
```

Check out the number of occupied UTM squares by the species over time. First the naive occupancy, that is the number of sites that were observed occupied over the two periods:
```{r}
naiveocc <- rep(NA,2) 
for (i in 1:2){
  naiveocc[i] <- sum(y[,i]>0)
}
naiveocc
```

Now the estimated occupancy:
```{r}
zz <- out$BUGSoutput$sims.list$z
dim(zz) # 2000 x 2850 x 2 (nbMCMC x nSquares x nyears)
nbpixocc <- apply(zz,c(1,3),sum) # 2000 x 2
meannbpixocc <- apply(nbpixocc,2,median)
sdnbpixocc <- apply(nbpixocc,2,sd)
qulnbpixocc <- apply(nbpixocc,2,quantile,probs=2.5/100)
quunbpixocc <- apply(nbpixocc,2,quantile,probs=97.5/100)
```

Put everything together:
```{r}
occupancy <- data.frame(period = c(1,2),
  naiveocc = naiveocc, 
  medianocc = meannbpixocc,
  qulnbpixocc = qulnbpixocc,
  quunbpixocc = quunbpixocc)
occupancy2 <- data.frame(period = rep(c(1,2),2),
  occ = c(occupancy$naiveocc, occupancy$medianocc), 
  quantlower = c(rep(NA,2),occupancy$qulnbpixocc),
  quantupper = c(rep(NA,2), occupancy$quunbpixocc),
  Occupancy = c(rep('naive',2),rep('estimated',2)))
occupancy2
```

Now plot:
```{r}
occupancy2 %>% 
  ggplot() + 
  aes(x = period, y =  occ, group = Occupancy, linetype = Occupancy) + 
    geom_ribbon(data=occupancy2,aes(ymin=quantlower,ymax=quantupper),alpha=0.1, show.legend = FALSE) + 
 geom_line() +
  scale_linetype_manual(values=c("solid", "dashed")) +  
  geom_point() + 
  xlab('Year') + 
  ylab('Number of occupied sites')
```

Now display local extinction, colonization and species detection probabilities estimates with credible intervals:
```{r}
#install.packages("remotes")
#remotes::install_github("mikemeredith/MMmisc")
epsmean <- MMmisc::invcloglog(out$BUGSoutput$mean$mueps + mean(offset))
epsql <- MMmisc::invcloglog(out$BUGSoutput$mean$mueps + mean(offset) - 2*out$BUGSoutput$sd$mueps)
epsqu <- MMmisc::invcloglog(out$BUGSoutput$mean$mueps + mean(offset) + 2*out$BUGSoutput$sd$mueps)

gammean <- MMmisc::invcloglog(out$BUGSoutput$mean$mugam + mean(offset))
gamql <- MMmisc::invcloglog(out$BUGSoutput$mean$mugam + mean(offset) - 2*out$BUGSoutput$sd$mugam)
gamqu <- MMmisc::invcloglog(out$BUGSoutput$mean$mugam + mean(offset) + 2*out$BUGSoutput$sd$mugam)

pmean <- plogis(out$BUGSoutput$mean$mup)
pqu <- plogis(out$BUGSoutput$mean$mup + 2*out$BUGSoutput$sd$mup)
pql <- plogis(out$BUGSoutput$mean$mup - 2*out$BUGSoutput$sd$mup)

eps <- data.frame(
  param = epsmean, 
  qlo = epsql,
  qup = epsqu)

gam <- data.frame(
  param = gammean, 
  qlo = gamql,
  qup = gamqu)

det <- data.frame(period = c(1,2),
  param = pmean,
  qlo = pql,
  qup = pqu)

eps
gam
det
```

Let's map the estimated distribution of monk seals. First, compute realized occupancy per site and per period, and put altogether in a big table:
```{r}
zz <- out$BUGSoutput$sims.list$z # 2000 x 2850 x 2 (nbMCMC x nSquares x nyears)
meanz <- apply(zz,c(2,3),mean)
dim(meanz)
toplot <- data.frame(
  meanz = c(meanz[,c(1,2)]), 
  UTM_NAME = rep(unique(grid$UTM_NAME),2), 
  yearr = c(rep('2000-2007',nrow(meanz)), rep('2013-2020',nrow(meanz)))) 
ring_occ <- left_join(grid,toplot)
```

```{r}
nrow(meanz)
sum(meanz[,2] == meanz[,1])
sum(meanz[,2] > meanz[,1])
sum(meanz[,2] < meanz[,1])
```

Then, plot two maps, one for each multi-year period:
```{r}
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = grid, colour = "grey50", fill = "white", lwd = 0.1) + 
  geom_sf(data = ring_occ, lwd = 0.1, aes(fill = as.numeric(meanz))) + 
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
  scale_fill_viridis(
    name = 'Pr(occupancy)', 
    direction = -1,
    alpha = 0.7) + 
  labs(title = '', x = NULL, y = NULL) +
  facet_wrap(~ yearr, ncol = 2) + 
  theme(legend.position = 'bottom', 
        legend.title = element_text(size=8), 
        legend.text = element_text(size=8),
        legend.key.size = unit(0.5, "cm"))
ggsave(here::here("figures","map2pups.png"),dpi=600)
```

We'd also like to have two maps of the estimated occupied sites for the two periods. To do so, we say that all squares with prob of occupancy greater than $25\%$ are occupied. This gives the following maps: 
```{r}
zz <- out$BUGSoutput$sims.list$z # 2000 x 2850 x 2 (nbMCMC x nSquares x nyears)
sumz <- apply(zz,c(2,3),sum)
meanz <- apply(zz,c(2,3),sum)/dim(zz)[1]*100
z_occupied <- (meanz > 25) # 
toplot <- data.frame(
  z_occupied = c(z_occupied[,c(1,2)]), 
  UTM_NAME = rep(unique(grid$UTM_NAME),2), 
  yearr = c(rep('2000-2007',nrow(z_occupied)), rep('2013-2020',nrow(z_occupied)))) 
ring_occ <- left_join(grid,toplot) %>%
  mutate(z_occupied = if_else(z_occupied == TRUE, 'used', 'unused by the species'))

ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = ring_occ, lwd = 0.1, aes(fill = z_occupied)) + 
    geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
    scale_fill_manual(values = c('white','blue'),
                    name = "Site is") + 
  labs(x = NULL, y = NULL) +
  facet_wrap(~ yearr, ncol = 2) + 
  theme(legend.position = "none")
ggsave(here::here("figures","map3pups.png"),dpi=600)
```

```{r echo = FALSE}
left_join(grid,toplot) %>%
  filter(z_occupied == TRUE) %>%
  group_by(yearr) %>%
  summarize(total_area = sum(areakm2)) 
```


We'd like to consider the difference in occupancy probabilities.
```{r}
delta <- as.matrix(meanz[,2]/100 - meanz[,1]/100)
ring_occ2 <- ring_occ %>%
  filter(yearr == '2013-2020')
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = grid, colour = "grey50", fill = "white", lwd = 0.1) + 
  geom_sf(data = ring_occ2, lwd = 0.1, aes(fill = delta)) + 
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(15,"RdYlGn"),
                       name = '') + 
  labs(title = '', x = NULL, y = NULL) + 
  theme(legend.position = 'bottom', 
        legend.title = element_text(size=8), 
        legend.text = element_text(size=8),
        legend.key.size = unit(0.5, "cm"))
ggsave(here::here("figures","map4pups.png"),dpi=600)
```

Look at where the high probability of presence on monk seals pups coincides with marine protected areas, for the period 2013-2020.
```{r}
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = ring_occ2, lwd = 0.1, aes(fill = z_occupied)) +
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
    scale_fill_manual(values = c('white','blue'),
                    name = "Site is") + 
  geom_sf(data = pa, colour = "red", fill = "transparent", lwd = 0.1) +
  xlab("") + ylab("") +
  theme(legend.position = "none")
ggsave(here::here("figures","mpamonksealspups.png"),dpi=600)
```

Out of the total number of cells with high probability of occurrence of monk seals pups, this number intersects with marine protected areas.
```{r}
mask <- ring_occ %>%
  filter(z_occupied == 'used') %>%
  filter(yearr == '2013-2020') %>%
  st_intersects(pa) %>%
  lengths > 0
sum(mask)
ring_occ2$map <- ring_occ2$z_occupied
mask2 <- (ring_occ2$map == 'used')
ring_occ2$map[mask2][mask] <- 'used, overlap w/ mpa'
ring_occ2$map[mask2][!mask] <- 'used, no overlap w/ mpa'
```

On a map, this gives.
```{r}
ggplot() + 
  geom_sf(data = greece, colour = "grey50", fill = "white") + 
  geom_sf(data = ring_occ2, lwd = 0.1, aes(fill = map)) +
  geom_sf(data = greece, colour = "grey50", fill = "white",lwd = 0.2) + 
    scale_fill_manual(values = c('white','blue','red'),
                    name = "Site is") + 
  geom_sf(data = pa, colour = "red", fill = "transparent", lwd = 0.1) +
  xlab("") + ylab("") +
  theme(legend.position = "bottom")
ggsave(here::here("figures","mpamonksealsdetailspups.png"),dpi=600)
```