analyses.Rmd

---
title: "Predicting reintroduction sites using habitat selection in an endangered mesopredator"
author: "Kathryn Macpherson, Belinda Wilson, Shoshana Rapley"
date: "January 12, 2024"
output: 
  html_document:
    toc: true
    number_sections: true
    top_depth: 3
    toc_float:
      collapsed: true
    theme: cerulean
    highlight: pygments
editor_options:
  chunk_output_type: console
knit: (function(inputFile, encoding) { rmarkdown::render(inputFile, encoding = encoding, output_file = file.path(dirname(inputFile), 'tutorial.html')) })
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, eval=FALSE)
```

# **Background**

Reintroductions are a key strategy for reversing biodiversity loss (Sheean *et al*. 2012) and restoring ecological function (Lipsey & Child 2007). Their success can often be attributed to habitat suitability at the recipient site (Bennett *et al*. 2012), since it drives resource (i.e., food, shelter, protection) and mate availability (Su *et al*. 2021). It is important to confirm the habitat suitability of a recipient site, often through monitoring the movement and habitat selection of founders (individuals released for a reintroduction, Bennett *et al*. 2012).

We reintroduced eastern quolls (*Dasyurus viverrinus*) in a series of trials over three consecutive years (2016–2018) into a conservation-fenced sanctuary as part of the [Mulligans Flat-Goorooyarroo Woodland Experiment](https://www.coexistenceconservationlab.org/mulligans-flat-goorooyarroo-woodland-experiment). Throughout these trials, we radiotracked founders to reveal >200 diurnal dens.

Here we explore the habitat selection of reintroduced eastern quolls to predict suitable habitat for future translocations, using species distribution models (SDMs).

# **Setup**

First, we manually installed the [pacman Package Management Tool](https://cran.r-project.org/web/packages/pacman/index.html), which enables us to install and load subsequent packages in a condensed and efficient way.

We also manually installed `Rtools` (after having [downloaded the most up-to-date version from CRAN](https://cran.r-project.org/bin/windows/Rtools/)), and the `slga` (soils) package from GitHub.

These manual installations only need to be done once, upon first running the project.

```{r, eval=FALSE}
# Manually install pacman package
install.packages(pacman)

# Manually install Rtools package
installr::install.Rtools(check=TRUE, check_r_update=TRUE, GUI=TRUE)

# Manually install slga package
devtools::install_github("obrl-soil/slga")
```

```{r}
# Install and load required packages
pacman::p_load(beepr, data.table, devtools, dplyr, elevatr, forcats, ggmap, ggplot2, ggpubr, ggsn, ggspatial, freqtables, installr, janitor, lme4, lubridate, MuMIn, ncdf4, plotrix, plotROC, raster, rasterVis, readxl, reshape2, rgdal, rgl, rJava, rstudioapi, scales, SDMtune, sf, slga, st, stringr, terra, tidyr, tidyterra, tidyverse, viridis, zeallot)
```

We also set the working directory to where this R markdown is saved using the `rstudioapi` package.

```{r}
# Set working directory to where your Rmarkdown is saved
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
```

# **Data preparation**

Here we read in our eastern quoll den assessment data, and compare each of our *intrinsic* den covariates with what is available throughout MFWS. These are:

  1.  **Den type** (open, hole, grass, log, pipe, drain, or culvert)
  2.  **Den cover** (ground, rocky outcrop, log, stump, fallen tree roots, living tree, dead tree, termite mound, pipe, drain, culvert, or building)
  3.  **Den substrate** (natural, experiment, or artificial)
  4.  **Number of entrances** (1, 2, or ≥3)
  5.  **Plot type** (open, or kangaroo exclusion zone)
  6.  **Habitat** (open grassland, Eucalypt forest, or Eucalypt woodland)
  7.  **Activity** (active, inactive, or could not locate)

## Import

```{r}
# Assign raw data filename to an object
raw_data <- "data.xlsx"

# Read in data
data <- read_excel("input/raw data.xlsx", sheet="dens") %>%
  clean_names() %>%
  # Convert numeric dates to date format
  mutate(date = as.Date(date, origin="1899-12-30"), 
         entrances = as.factor(entrances_holes)) %>%
  dplyr::select(den_id, sector, "den_type"=type, 
                "den_cover"=under, "den_substrate"='in', 
                entrances_holes, plot_type, habitat, activity, 
                easting, northing, latitude, longitude) %>%
  subset(sector!="Outside")
```

Calculate how many dens were first found in each study year.

```{r}
years <- read_excel("input/raw data.xlsx", sheet="visits") %>%
  clean_names() %>%
  mutate(year=year(as.Date(date_tracked))) %>%
  group_by(eqd_code) %>%
  summarise(min(year)) %>%
  clean_names() %>%
  group_by(min_year) %>%
  count()
```

As expected, the number of dens first identified decreased throughout the years (119  in 2016, 79 in 2017, and 33 in 2018). 

### Den map

Here we generate a map of the den sites coloured in order of first use: oldest dens are purple, middle dens are green, and newest dens are yellow (ranging February 2016–November 2018).

You will need to [generate your own API key](https://developers.google.com/maps/documentation/javascript/get-api-key) to fetch the basemap from Google.

```{r,echo=FALSE, eval=FALSE}
# Read in MFWS fence shapefile for plotting
mfws <- readOGR(dsn="input/mfws_fence.shp", verbose=FALSE) %>%
  spTransform(CRS("+proj=utm +zone=55 +datum=WGS84")) %>%
  fortify(verbose=FALSE) %>% #transforms from from sp to df
  mutate(lat=as.numeric(lat + 10000000), 
         long=as.numeric(long))

# Use Google API to fetch a base map
#ggmap::register_google(key="[enter API key here]")
sat <- get_map(location=c(lat=-35.16580, lon=149.16454),
               zoom=14, source="google", maptype="satellite", crop=TRUE)

# Display the map
ggmap(sat)

# Read in MFWS fence shapefile for plotting
mf_outline <- readOGR(dsn="shapefiles/mfws_fence.shp", verbose=FALSE) %>%
  spTransform("EPSG:4326") %>%
  fortify(verbose=FALSE) #transforms from from sp to df

# Read in MFWS fence shapefile for plotting
mf_tracks <- readOGR(dsn="shapefiles/MFWS_tracks.shp", verbose=FALSE) %>%
  spTransform("EPSG:4326") %>%
  fortify(verbose=FALSE) #transforms from from sp to df

# Create a map
den_map_sat <- ggmap(sat) +
  geom_path(mf_outline, mapping=aes(x=long, y=lat, group=group),
            col="white", size=1) +
  geom_path(mf_tracks, mapping=aes(x=long, y=lat, group=group),
            col="white", size=0.5) +
  geom_point(data, mapping=aes(x=longitude, y=latitude, 
                               col=den_id), size=2, alpha=0.8) +
  ggsn::scalebar(y.min=-35.153, x.min=149.145, 
                 y.max=-35.182, x.max=149.190,
                 dist=500, height=3, dist_unit="m",
                 box.color="white", st.color="white",
                 st.size=3, anchor=c(y=-35.175, x=149.182),
                 st.dist=0.04, transform=TRUE) +
  theme(panel.grid.major=element_blank(),
        panel.background=element_rect(fill="white"),
        plot.margin=unit(c(0, 0, 0, 0), "cm"),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.x=element_blank(),
        axis.ticks.y=element_blank(),
        legend.position = "none", 
        legend.key=element_blank()) +
  scale_colour_viridis(discrete=T, option="D", begin=0, end=1) +
  xlab("") + ylab("")

# Display the plot
print(den_map_sat)
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(den_map_sat, filename="output/den satellite map.jpeg", 
       width=200, height=125, units="mm")
```

# **Intrinsic den covariates**

## Den type

  a) Plots

```{r}
# Transfer from wide to long format, and calculate percent for each level
long <- data %>%
  group_by(den_type) %>%
  summarise(percent = 100/nrow(data)*length(den_type), 
            n = length(den_type)) %>%
  na.omit() %>% mutate(den_type = factor(den_type, 
                                         levels=c("Natural", 
                                                  "Experiment", 
                                                  "Artificial")))

# Create histogram plot
den_type_hist <- ggplot(long, aes(x=den_type, y=percent, fill=den_type)) +
  geom_col() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Den type") + ylab("Percentage (%)") +
  # Viridis scale options are A–E
  scale_fill_viridis(discrete=T, option="D", begin=0.9, end=0.2) +
  scale_y_continuous(limits=c(0, 100), breaks=seq(0, 100, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))


# Histogram for publication
den_type_hist_pub <- ggplot(long, aes(x=den_type, y=percent, fill=den_type)) +
  geom_col() +
  annotate("text", x=3.5, y=90, label = parse(text = "bold('a')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Den type") + ylab("Percentage (%)") +
  # Viridis scale options are A–E
  scale_fill_viridis(discrete=T, option="D", begin=0.1, end=0.9) +
  scale_y_continuous(limits=c(0, 90), breaks=seq(0, 90, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Display the plot
print(den_type_hist_pub)

# Subset to those records with valid den substrates
plot_data <- subset(data, den_type!="NA")
View(data)
# Create map plot
den_type_map <- ggplot() + 
  geom_path(mfws, mapping=aes(x=long, y=lat, group=group), 
            col="grey15") +
  geom_point(plot_data, mapping=aes(x=easting, y=northing, 
                                    col=den_type), alpha=0.8) +
  coord_sf(xlim=c(695750, 699250), ylim=c(6104250, 6107750)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.70, 0.25),
        legend.title = element_text(size=9), 
        legend.text = element_text(size=9),
        legend.key.size = unit(0.5, 'cm'),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(-0.5, -1, -1, -1), "cm"),
        legend.key = element_blank()) +
  scale_colour_manual(values=viridis(length(unique(data$den_type)), 
                                     begin=0.9, end=0.2), name="Den type",
                      breaks=c('Natural', 'Experiment', 'Artificial')) +
  xlab("") + ylab("")

# Combine histogram and map
den_type_plots <- ggarrange(den_type_hist, den_type_map, 
                            ncol=2, nrow=1, widths=c(3,3))

# Display the plot
print(den_type_plots)
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(den_type_plots, 
       filename="output/den type histogram and map.jpeg", 
       width=200, height=75, units="mm")

# Export the plot to a jpeg file
ggsave(den_type_hist_pub, 
       filename="output/den type histogram publication.jpeg", 
       width=150, height=75, units="mm")
```

  b) Models

```{r}
# Chi-square test
den_type_mod <- data %>%
  subset(den_type!="N/A") %>%
  freq_table(den_type) %>%
  freq_test() %>%
  select(c(1:3, 5, 7, 12:14))
```

## Den cover

  a) Plots

```{r}
# Transfer from wide to long format, and calculate percent for each level
long <- data %>%
  group_by(den_cover) %>%
  summarise(percent = 100/nrow(data)*length(den_cover), 
            n = length(den_cover)) %>% 
  na.omit() %>%
  mutate(den_cover=factor(den_cover, 
                          levels=c("Ground", "Log", 
                                   "Fallen tree roots", "Living tree", 
                                   "Pipe, drain, or culvert", "Stump",
                                   "Dead tree", "Rocky outcrop", "Building",
                                   "Termite mound")))

# Create histogram plot
den_cover_hist <- ggplot(long, aes(x=den_cover, y=percent, fill=den_cover)) +
  geom_col() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        legend.position = "none",
        axis.text.x = element_text(angle=-45, vjust=0.5, hjust=0.2),
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey"),
  plot.margin = unit(c(0.5, 0.75, 0, 0.2), "cm")) +
  xlab("Den cover") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=1, end=0) +
  scale_y_continuous(limits=c(0, 50), breaks=seq(0, 100, 20)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))


# Histogram for publication
den_cover_hist_pub <- ggplot(long, aes(x=den_cover, y=percent, fill=den_cover)) +
  geom_col() +
  annotate("text", x=10.5, y=50, label = parse(text = "bold('c')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(),
        plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
        legend.position = "none",
        axis.text.x = element_text(angle=-45, vjust=0.5, hjust=0.2),
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Den cover") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=0, end=1) +
  scale_y_continuous(limits=c(0, 50), breaks=seq(0, 50, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Display the plot
print(den_cover_hist_pub)

# Subset to those records with valid den covers
plot_data <- subset(data, den_cover!="NA")

# Create map plot
den_cover_map <- ggplot() + 
  geom_path(mfws, mapping=aes(x=long, y=lat, group=group), col="grey15") +
  geom_point(plot_data, mapping=aes(x=easting, y=northing, 
                                    col=den_cover), alpha=0.8) +
  coord_sf(xlim=c(695750, 699250), ylim=c(6104250, 6107750)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.85, 0.25),
        legend.title = element_text(size=9), 
        legend.text = element_text(size=9),
        legend.key.size = unit(0.5, 'cm'),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(-0.5, -1, -1, -1), "cm"),
        legend.key = element_blank()) +
  scale_colour_manual(values=viridis(length(unique(data$den_cover)), 
                                     begin=1, end=0), name="Den cover") + 
  guides(colour=guide_legend(ncol=2)) +
  xlab("") + ylab("")

# Combine histogram and map into single plot
den_cover_plots <- ggarrange(den_cover_hist, den_cover_map, 
                            ncol=2, nrow=1, widths=c(3,3))

# Display the plot
print(den_cover_plots)
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(den_cover_plots, 
       filename="output/den cover histogram and map.jpeg", 
       width=275, height=125, units="mm")

# Export the plot to a jpeg file
ggsave(den_cover_hist_pub, 
       filename="output/den cover histogram publication.jpeg", 
       width=150, height=75, units="mm")
```

  b) Models

```{r}
# Chi-square test
den_cover_mod <- data %>%
  subset(den_cover!="N/A") %>%
  freq_table(den_cover) %>%
  freq_test() %>%
  select(c(1:3, 5, 7, 12:14))
```

## Den substrate

  a) Plots

```{r}
# Transfer from wide to long format, and calculate percent for each level
long <- data %>%
  group_by(den_substrate) %>%
  summarise(percent = 100/nrow(data)*length(den_substrate), 
            n = length(den_substrate)) %>%
  na.omit() %>%
  mutate(den_substrate=factor(den_substrate, 
                          levels=c("Rabbit warren", "Log", 
                                   "Grass", "Pipe, drain, or culvert", 
                                   "Open")))

# Use ifelse to define 'hole' as 'rabbit warren' instead
data$den_substrate <- ifelse(data$den_substrate == "Hole", "Rabbit warren", as.character(data$den_substrate))

# Create histogram plot
den_substrate_hist <- ggplot(long, aes(x=den_substrate, y=percent, 
                                       fill=den_substrate)) +
  geom_col() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Den substrate") + ylab("Percentage (%)") +
  # Options include magma, inferno, plasma, viridis, turbo
  scale_fill_viridis(discrete=T, option="D", begin=1, end=0) +
  scale_y_continuous(limits=c(0, 90), breaks=seq(0, 100, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Histogram for publication
den_substrate_hist_pub <- ggplot(long, aes(x=den_substrate, y=percent, 
                                       fill=den_substrate)) +
  geom_col() +
  annotate("text", x=5.5, y=90, label = parse(text = "bold('b')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(),
        plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Den substrate") + ylab("Percentage (%)") +
  # Options include magma, inferno, plasma, viridis, turbo
  scale_fill_viridis(discrete=T, option="D", begin=0, end=1) +
  scale_y_continuous(limits=c(0, 90), breaks=seq(0, 90, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Display the plot
print(den_substrate_hist_pub)

# Subset to those records with valid den substrates
plot_data <- subset(data, den_substrate!="NA")

# Create map plot
den_substrate_map <- ggplot() + 
  geom_path(mfws, mapping=aes(x=long, y=lat, group=group), 
            col="grey15") +
  geom_point(plot_data, mapping=aes(x=easting, y=northing, 
                                    col=den_substrate), alpha=0.8) +
  coord_sf(xlim=c(695750, 699250), ylim=c(6104250, 6107750)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.70, 0.25),
        legend.title = element_text(size=9), 
        legend.text = element_text(size=9),
        legend.key.size = unit(0.5, 'cm'),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(-0.5, -1, -1, -1), "cm"),
        legend.key = element_blank()) +
  scale_colour_manual(values=viridis(length(unique(data$den_substrate)), 
                                     begin=1, end=0), name="Den substrate") +
  xlab("") + ylab("")

# Combine histogram and map into single plot
den_substrate_plots <- ggarrange(den_substrate_hist, 
                                 den_substrate_map, 
                                 ncol=2, nrow=1, widths=c(3,2))

# Display the plot
print(den_substrate_plots)
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(den_substrate_plots, 
       filename="output/den substrate histogram and map.jpeg", 
       width=250, height=125, units="mm")

# Export the plot to a jpeg file
ggsave(den_substrate_hist_pub, 
       filename="output/den substrate histogram publication.jpeg", 
       width=150, height=75, units="mm")
```

  b) Models

```{r}
# Chi-square test
den_substrate_mod <- data %>%
  subset(den_substrate!="N/A") %>%
  freq_table(den_substrate) %>%
  freq_test() %>%
  select(c(1:3, 5, 7, 12:14))
```

## Number of entrances

```{r}
# Use ifelse to define entrance hole values
data_2 <- data %>%
  mutate(entrances = ifelse(entrances_holes==1, 1,
                     ifelse(entrances_holes==2, 2,
                     ifelse(entrances_holes=="≥3", "≥3", NA)))) %>%
  drop_na(entrances) %>%
  mutate(entrances = factor(entrances, levels=c("1", "2", "≥3")))
```

  a) Plots

```{r}
# Transfer from wide to long format, and calculate percent for each level
long <- data_2 %>%
  group_by(entrances) %>%
  summarise(percent = 100/nrow(data)*length(entrances), 
            n = length(entrances)) %>%
  na.omit()

# Create histogram plot
entrances_hist <- ggplot(long, aes(x=entrances, y=percent, 
                                       fill=entrances)) +
  geom_col() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Number of entrances") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=0.9, end=0.2) +
  scale_y_continuous(limits=c(0, 100), breaks=seq(0, 100, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Histogram for publication
entrances_hist_pub <- ggplot(long, aes(x=entrances, y=percent, 
                                       fill=entrances)) +
  geom_col() +
  annotate("text", x=3.5, y=80, label = parse(text = "bold('d')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(),
        plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Number of entrances") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=0.1, end=0.9) +
  scale_y_continuous(limits=c(0, 80), breaks=seq(0, 80, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Display the plot
print(entrances_hist_pub)

# Create map plot
entrances_map <- ggplot() + 
  geom_path(mfws, mapping=aes(x=long, y=lat, group=group), 
            col="grey15") +
  geom_point(data_2, mapping=aes(x=easting, y=northing, 
                                    col=entrances)) +
  coord_sf(xlim=c(695750, 699250), ylim=c(6104250, 6107750)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.75, 0.25), 
        legend.title = element_text(size=9), 
        legend.text = element_text(size=9),
        legend.key.size = unit(0.5, 'cm'),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(-0.5, -1, -1, -1), "cm"),
        legend.key = element_blank()) +
  scale_colour_manual(values=viridis(length(unique(data_2$entrances)), 
                                     begin=0.9, end=0.2), 
                                     name="Number of entrances") +
                                     xlab("") + ylab("")

# Combine histogram and map into single plot
entrances_plots <- ggarrange(entrances_hist, 
                                 entrances_map, 
                                 ncol=2, nrow=1, widths=c(3,2))

# Display the plot
print(entrances_plots)
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(entrances_plots, 
       filename="output/entrances histogram and map.jpeg", 
       width=250, height=125, units="mm")

# Export the plot to a jpeg file
ggsave(entrances_hist_pub, 
       filename="output/entrances histogram publication.jpeg", 
       width=150, height=75, units="mm")
```

  b) Models

```{r}
# Chi-square test
entrances_mod <- data %>%
  subset(entrances_holes!="N/A") %>%
  freq_table(entrances_holes) %>%
  freq_test() %>%
  select(c(1:3, 5, 7, 12:14))
```

## Plot type

### Plots

```{r}
# Transfer from wide to long format, and calculate percent for each level
long <- data %>%
  group_by(plot_type) %>%
  summarise(percent = 100/nrow(data)*length(plot_type), 
            n = length(plot_type)) %>%
  na.omit() %>% mutate(plot_type = factor(plot_type, 
                                          levels=c("Open", 
                                                   "Kangaroo exclosure", 
                                                   "Bettong exclosure")))

# Create histogram plot
plot_type_hist <- ggplot(long, aes(x=plot_type, y=percent, 
                                       fill=plot_type)) +
  geom_col() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Plot type") + ylab("Percentage (%)") +
  # Options include magma, inferno, plasma, viridis, turbo
  scale_fill_viridis(discrete=T, option="D", begin=0.9, end=0.1) +
  scale_y_continuous(limits=c(0, 100), breaks=seq(0, 100, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Histogram for publication
plot_type_hist_pub <- ggplot(long, aes(x=grazing_pressure, y=percent, 
                                       fill=grazing_pressure)) +
  geom_col() + 
  annotate("text", x=3.5, y=90, label = parse(text = "bold('b')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(),
        plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Grazing pressure") + ylab("Percentage (%)") +
  # Options include magma, inferno, plasma, viridis, turbo
  scale_fill_viridis(discrete=T, option="D", begin=0.1, end=0.9) +
  scale_y_continuous(limits=c(0, 100), breaks=seq(0, 100, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Display the plot
print(plot_type_hist_pub)

# Subset to those records with valid plot types
plot_data <- subset(data, plot_type!="NA")

# Create map plot
plot_type_map <- ggplot() + 
  geom_path(mfws, mapping=aes(x=long, y=lat, group=group), 
            col="grey15") +
  geom_point(plot_data, mapping=aes(x=easting, y=northing, 
                                    col=plot_type)) +
  coord_sf(xlim=c(695750, 699250), ylim=c(6104250, 6107750)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.75, 0.25),
        legend.title = element_text(size=9), 
        legend.text = element_text(size=9),
        legend.key.size = unit(0.5, 'cm'),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(-0.5, -1, -1, -1), "cm"),
        legend.key = element_blank()) +
  scale_colour_manual(values=viridis(length(unique(data$plot_type)), 
                                     begin=0.9, end=0.1), name="Plot type", 
                      breaks=c('Open', 'Kangaroo exclosure', 
                               'Bettong exclosure')) + xlab("") + ylab("")

# Combine histogram and map into single plot
plot_type_plots <- ggarrange(plot_type_hist, 
                                 plot_type_map, 
                                 ncol=2, nrow=1, widths=c(3,2))

# Display the plot
print(plot_type_plots)
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(plot_type_plots, 
       filename="output/plot type histogram and map.jpeg", 
       width=250, height=125, units="mm")

# Export the plot to a jpeg file
ggsave(plot_type_hist_pub, 
       filename="output/plot type histogram publication.jpeg", 
       width=150, height=75, units="mm")
```

### Models

```{r}
# Chi-square test
plot_type_mod <- data %>%
  subset(plot_type!="N/A") %>%
  freq_table(plot_type) %>%
  freq_test() %>%
  select(c(1:3, 5, 7, 12:14))
```

### Adjust for availability within MFWS
  
```{r}
# Read in MFWS fence shapefile
mfws <- st_read("input/mfws_fence.shp")

# Reproject and prepare shapefile
mfws <- mfws %>%
  st_set_crs(st_crs(mfws)) %>%
  st_transform("EPSG:32755") %>% 
  clean_names() %>%
  mutate(plot="open") %>%
  # Since the roo layer only includes three column, we select them here
  select(id, geometry, plot)

# Read in MFWS kangaroo exclusion zone shapefile using sf
kan <- st_read("input/MFGO_kangaroo_exclosures_polygons.shp")

# Reproject and prepare shapefile
kan <- kan %>%
  st_set_crs(st_crs(kan)) %>%
  st_transform("EPSG:32755") %>% 
  clean_names() %>%
  mutate(plot="kangaroo", 
         id=row_number()) %>%
  # Remove exclosures in the adjacent reserve
  filter(id <6) %>%
  # Mutate id so it has prefix
  mutate(id=paste0("K", row_number()))

# Read in MFWS kangaroo exclusion zone shapefile using sf
bet <- st_read("input/MFWS_bettong_exclosures_polygons.shp") 

# Reproject and prepare shapefile
bet <- bet %>%
  st_set_crs(st_crs(bet)) %>%
  st_transform("EPSG:32755") %>% 
  clean_names() %>%
  # Since the roo layer only includes three column, we select them here
  select(id, geometry) %>%
  # Mutate id so it has prefix
  mutate(id=paste0("B", row_number())) %>%
  # Assign bettong plots inside kangaroo plots as another category
  mutate(id = ifelse(id == "B1", "KB1", 
              ifelse(id == "B6", "KB6", 
              ifelse(id == "B8", "KB8", 
              ifelse(id == "B9", "KB9",
              ifelse(id == "B10", "KB10", 
              ifelse(id == "B12", "KB12", id))))))) %>%
  # Assign zone type to these plots
  mutate(plot = ifelse(str_starts(id, "^B"), "bettong", "both"))

# Crop MFWS polygon by the shape of kan and bet layers
mfws_clip <- st_difference(mfws, st_union(kan, bet)) %>%
  select(id=id.1, plot=plot.1, geometry)

# Plot the plot types
ggplot() + 
  geom_sf(data=mfws_clip, aes(fill=plot)) + 
  geom_sf_text(data=mfws_clip, aes(label=plot)) +

  geom_sf(data=kan, aes(fill=id)) + 
  geom_sf_text(data=kan, aes(label=id)) +
  geom_sf(data=bet, aes(fill=id)) + 
  geom_sf_text(data=bet, aes(label=id)) +
  #coord_sf(expand=FALSE, xlim=c(extent(mfws)[1]+100, 
  #                             extent(mfws)[2]+100), 
  #                       ylim=c(extent(mfws)[3]+100, 
  #                              extent(mfws)[4]+100)) +
  theme_minimal() + 
  theme(legend.position = "none", 
        axis.title.x = element_blank(), 
        axis.title.y = element_blank())

# Combine kangaroo and bettong layers and calculate area
plots <- rbind(kan, bet) %>%
  # Calculate area of each geometry (or polygon)
  mutate(area_km2=as.numeric(st_area(geometry))/1e6)

plots_sum <- plots %>%
  as.data.frame() %>% 
  group_by(plot) %>%
  summarise(area_km2=sum(area_km2)) %>%
  select(plot, area_km2)

# Calculate the open plot area
plots_sum <- plots_sum %>%
  add_row(plot="open", area_km2=mfws_area_minus_aquatic-
            sum(plots$area_km2))

# Read in den coordinates data and prepare for plotting
den <- read_excel("input/raw data.xlsx", sheet="dens") %>% 
  clean_names() %>% 
  mutate(easting=as.numeric(easting), 
         northing=as.numeric(northing))

# Convert to spatial points dataframe
den_sp <- SpatialPointsDataFrame( 
  data.frame(den$easting, den$northing), den,
  proj4string=CRS("EPSG:32755")) %>% 
  st_as_sf()

# Join den coordinates and plots layer
den_zone <- st_join(den_sp, plots, join=st_intersects) %>%
  st_drop_geometry() %>%
  as.data.frame() %>%
  mutate(plot = ifelse(is.na(plot), "open", plot)) %>%
  # Remove dens outside MFWS
  filter(sector != "Outside", 
         !is.na(sector))

# Extract exclosure zone for each den coordinate and calculate frequencies
den_zone_adj <- den_zone %>%
  group_by(plot) %>% 
  count() %>% 
  # Since only one den occurs in both, but is included in the kangaroo plot
  # count, we exclude it
  filter(plot != "both") %>%
  # Adjusted by frequency of dens
  mutate(den_freq=n) %>%
  select(-n) %>%
  mutate(den_perc=(den_freq/sum(den_zone_adj$den_freq))*100) %>%
  # Remove dens outside MFWS
  filter(!is.na(plot))

# Assign MFWS minus aquatic area value
mfws_area_minus_aquatic <- sum(mfws_group$area_km2)-mfws_group$area_km2[1]

# Assign correct number of dens (223)
dens_n <- length(unique(den_zone$den))

# Adjust frequencies and percentages to account for exclosure zone availability
den_zone_plot <- den_zone_adj %>%
  # Join with df to get area_km2 column
  left_join(plots, by="plot") %>%
  group_by(plot) %>% 
  summarize(plot = max(plot), den_freq = max(den_freq), 
             den_perc = max(den_perc), area_km2 = sum(area_km2)) %>% 
  mutate(area_km2 = ifelse(plot == "open", 
                       # Calculate open area by subtracting other plot types
                       mfws_area_minus_aquatic-
                         (sum(area_km2, na.rm=TRUE)), area_km2)) %>%
  # Calculate frequency of dens adjusted by available area 
  mutate(den_freq_adj=(den_freq/area_km2)*
                         mfws_area_minus_aquatic/
                            nrow(den_zone_adj)) %>%
  # Redistributing so the number of dens is 223
  mutate(den_freq_adj=den_freq_adj/sum(den_freq_adj)*dens_n) %>%
  # Calculate percentage of dens adjusted by available area
  mutate(den_perc_adj=(den_freq_adj/sum(den_zone_plot$den_freq_adj))*100, 
         plot=factor(plot, levels = c("open", "kangaroo", "bettong"), 
                           labels = c("Open", "Kangaroo", "Bettong")))

# Display chi-square test for unequal distribution of dens
# after having been adjusted for available area 
den_zone %>%
  # Remove any dens that were in zone NA (i.e., outside MFWS)
  filter(!is.na(plot), 
        plot != "both") %>%
  freq_table(plot) %>%
  mutate(n=den_zone_plot$den_freq_adj, 
         percent=den_zone_plot$den_perc_adj) %>%
  freq_test() %>%
  # Select the statistics we're interested in
  select(plot=cat, n, percent, t_crit, 
         chi2_pearson, df, p=p_chi2_pearson)
```

 Bar plot
 
```{r}
# Use ifelse to define plot type as grazing pressure instead
den_zone_plot <- den_zone_plot %>%
  mutate(grazing_pressure = ifelse(plot == "Open", "High",
                          ifelse(plot == "Kangaroo", "Medium",
                                 ifelse(plot == "Bettong", "Low", NA))))

# Ensure columns are ordered correctly for plotting
den_zone_plot$grazing_pressure <- factor(den_zone_plot$grazing_pressure,
                                         levels = c("High", "Medium", "Low"))
```
 
 
```{r}
# Plot raw den percentages per habitat type
zone_raw_barplot <- ggplot(data=den_zone_plot, 
  aes(x= grazing_pressure, y=den_perc, fill=grazing_pressure)) +
  geom_col() +
  annotate("text", x=3.5, y=90, label = parse(text = "bold('c')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"),
        axis.line.y = element_line(colour="grey")) +
  xlab("") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=0.1, end=0.9) +
  scale_y_continuous(limits=c(0, 90), breaks=seq(0, 90, 10)) 

# Display the plot
print(zone_raw_barplot)

# Plot adjusted den percentages per habitat type
zone_adj_barplot <- ggplot(data=den_zone_plot, 
  aes(x= grazing_pressure, y=den_perc_adj, fill=grazing_pressure)) +
  geom_col() +
  annotate("text", x=3.5, y=60, label = parse(text = "bold('d')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"),
        axis.line.y = element_line(colour="grey")) +
  xlab("") + ylab("Percentage adjusted for availability (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=0.1, end=0.9) +
  scale_y_continuous(limits=c(0, 60), breaks=seq(0, 60, 10)) 

# Display the plot
print(zone_adj_barplot)

# Combine the plots
zone_barplots <- ggarrange(zone_raw_barplot, zone_adj_barplot, ncol= 2, nrow= 1, 
                                   widths= c(1,1), heights = c(1,1)) %>%
  annotate_figure(bottom=text_grob("Grazing zone", hjust=0.5, vjust=0))


# Display the combined plot
print(zone_barplots)
```

```{r, echo=FALSE, eval=FALSE}
ggsave(zone_barplots, 
       filename="output/grazing pressure barplots.jpeg", 
       width=225, height=100, units="mm")
```

## Habitat

### Plots

```{r}
# Transfer from wide to long format, and calculate percent for each level
long <- plot_data %>%
  group_by(habitat) %>%
  summarise(percent = 100/nrow(plot_data)*length(habitat), 
            n = length(habitat))

# Create histogram plot
habitat_hist <- ggplot(long, aes(x=habitat, y=percent, fill=habitat)) +
  geom_col() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Habitat type") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=1, end=0) +
  scale_y_continuous(limits=c(0, 100), breaks=seq(0, 100, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Histogram for publication
habitat_hist_pub <- ggplot(long, aes(x=habitat, y=percent, fill=habitat)) +
  geom_col() +
  annotate("text", x=4.5, y=60, label = parse(text = "bold('a')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(),
        plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
        legend.position = "none",
        axis.line.x = element_line(colour="grey"), 
        axis.line.y = element_line(colour="grey")) +
  xlab("Habitat type") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=1, end=0) +
  scale_y_continuous(limits=c(0, 60), breaks=seq(0, 60, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Display the plot
print(habitat_hist_pub)

# Subset to those records with valid habitat types
plot_data <- subset(data, habitat!="NA")

# Create map plot
habitat_type_map <- ggplot() + 
  geom_path(mfws, mapping=aes(x=long, y=lat, group=group), col="grey15") +
  geom_point(plot_data, mapping=aes(x=easting, y=northing, 
                                    col=habitat), alpha=0.8) +
  coord_sf(xlim=c(695750, 699250), ylim=c(6104250, 6107750)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.70, 0.25),
        legend.title = element_text(size=9), 
        legend.text = element_text(size=9),
        legend.key.size = unit(0.5, 'cm'),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(-0.5, -1, -1, -1), "cm"),
        legend.key = element_blank()) +
  scale_colour_manual(values=viridis(length(unique(data$habitat)), 
                                     begin=1, end=0), name="Habitat type") +
  xlab("") + ylab("")

# Combine histogram and map into single plot
habitat_type_plots <- ggarrange(habitat_hist, habitat_type_map, 
                                ncol=2, nrow=1, widths=c(3,3))

# Display the plot
print(habitat_type_plots)
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(habitat_type_plots, 
       filename="output/habitat type histogram and map.jpeg", 
       width=250, height=100, units="mm")

# Export the plot to a jpeg file
ggsave(habitat_hist_pub, 
       filename="output/habitat type histogram publication.jpeg", 
       width=150, height=75, units="mm")
```

### Models

```{r}
# Chi-square test
habitat_type_mod <- data %>%
  subset(habitat!="N/A") %>%
  freq_table(habitat) %>%
  freq_test() %>%
  select(c(1:3, 5, 7, 12:14))
```

### Adjust for availability within MFWS

Map of NVIS habitat types (as per Wilson *et al*. *in review*).

```{r}
# Read in MFWS NVIS groups using sf
nvis <- st_read("input/mfws_nvis_vegetation_groups.shp") %>%
  st_transform("EPSG:4326") %>% 
  clean_names() %>% 
  mutate(mfws_name=factor(mfws_name, 
                          levels=c("Grassland", "Eucalypt forest", 
                                   "Eucalypt woodland", "Regrowth", 
                                   "Aquatic")))

# Map plot
nvis_mfws <- ggplot() + 
  geom_sf(st_zm(nvis), mapping=aes(fill=mfws_name), 
          col="grey30", lwd=0.5, alpha=0.8) + 
  annotation_scale(style="ticks", tick_height=0.5,
                   width_hint=0.3, 
                   pad_x=unit(1.55, "cm"),
                   pad_y=unit(0, "cm")) +
  #annotate(geom="text", y=-35.157, x=149.186, label="b)", fontface="bold", size=6)
  theme(axis.line.x = element_blank(),
        axis.line.y = element_blank(),
        legend.position = c(0.75, 0.15),
        legend.title = element_text(size=9), 
        legend.text = element_text(size=9),
        legend.key.size = unit(0.5, 'cm'),
        plot.margin = unit(c(0, 0, 1, 0), "cm"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.length = unit(0, "pt"),
        plot.title = element_text(hjust=0.5),
        text = element_text(colour="black"),
        axis.text.y = element_blank(),
        axis.text.x = element_blank()) +
  scale_fill_manual(values=viridis(5, begin=0.9, end=0.1)) +
  labs(fill="MFWS habitat", x=NULL, y=NULL) 
print(nvis_mfws)
```

```{r}
# Read in den coordinates data and prepare for plotting
den <- read_excel("input/raw data.xlsx", sheet="dens") %>% 
  clean_names() %>% 
  #na.omit() %>%
  mutate(easting=as.numeric(easting), 
         northing=as.numeric(northing), 
         date=as.POSIXct(date), 
         day=as.numeric(yday(date)))

# Convert to spatial points dataframe
den_sp <- SpatialPointsDataFrame( 
  data.frame(den$easting, den$northing), den,
  proj4string=CRS("EPSG:32755")) %>% 
  st_as_sf()

# Read in NVIS layer
nvis <- st_read("input/mfws_nvis_vegetation_groups.shp", quiet=TRUE) %>%
  st_transform("EPSG:32755")

# Join den coordinates and NVIS layer
den_nvis <- st_join(den_sp, nvis, join=st_intersects) %>%
  rename(mfws_habitat=MFWS_NAME) %>% 
  st_drop_geometry()

# Read in NVIS layer (warnings okay)
nvis_shp <- readOGR(dsn="input/mfws_nvis_vegetation_groups.shp") %>%
  spTransform(CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")) 
nvis_shp$area <- raster::area(nvis_shp)/1000000 #library(raster) in km2

# Calculate area per habitat type
mfws_group <- as.data.frame(nvis_shp) %>% 
  clean_names() %>%
  mutate(mfws_habitat=mfws_name) %>%
  group_by(mfws_habitat) %>%
  summarise(area_km2=sum(area))

# Assign MFWS minus aquatic area value
mfws_area_minus_aquatic <- sum(mfws_group$area_km2)-mfws_group$area_km2[1]

# Extract habitat type for each den coordinate and calculate frequencies
den_hab_adj <- den_nvis %>%
  group_by(mfws_habitat) %>% 
  count() %>% 
  # Adjusted by frequency of dens
  mutate(den_freq=n) %>%
  select(-n) %>%
  mutate(den_perc=(den_freq/sum(den_hab_adj$den_freq))*100) %>%
  # Remove dens outside MFWS
  filter(!is.na(mfws_habitat))

# Adjust frequencies and percentages to account for habitat availability
den_hab_plot <- den_hab_adj %>%
  # Join with df to get area_km2 column
  left_join(mfws_group, by="mfws_habitat") %>%
  # Calculate frequency of dens adjusted by available area 
  mutate(den_freq_adj=((den_freq/area_km2)*
                         mfws_area_minus_aquatic/nrow(den_hab_adj))) %>%
  # Calculate percentage of dens adjusted by available area
  mutate(den_perc_adj=(den_freq_adj/sum(den_hab_plot$den_freq_adj))*100, 
         mfws_habitat=factor(mfws_habitat, 
                             levels = c("Grassland",
                                        "Eucalypt forest",
                                        "Eucalypt woodland",
                                        "Regrowth",
                                        "Aquatic")))
# Reorder levels of mfws_habitat
den_hab_plot$mfws_habitat <- factor(den_hab_plot$mfws_habitat, 
                                 levels = c("Grassland", "Eucalypt forest", "Eucalypt woodland", "Regrowth", "Aquatic"))

# Display chi-square test results
den_nvis %>%
  # Remove any dens that were in habitat NA (i.e., outside MFWS)
  filter(!is.na(mfws_habitat)) %>%
  freq_table(mfws_habitat) %>%
  mutate(n=den_hab_plot$den_freq_adj, 
         percent=den_hab_plot$den_perc_adj) %>%
  freq_test() %>%
  # Select the statistics we're interested in
  select(mfws_habitat=cat, n, percent, t_crit, 
         chi2_pearson, df, p=p_chi2_pearson)
```

Bar plot

```{r}
# Plot raw den percentages per habitat type
hab_raw_barplot <- ggplot(data=den_hab_plot, 
  aes(x= mfws_habitat, y=den_perc, fill=mfws_habitat)) + 
  geom_col() +
  annotate("text", x=4.4, y=60, label = parse(text = "bold('a')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"),
        axis.line.y = element_line(colour="grey")) +
  xlab("") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=0.1, end=0.9) +
  scale_y_continuous(limits=c(0, 60), breaks=seq(0, 60, 10))

# Display the plot
print(hab_raw_barplot)

# Plot adjusted den percentages per habitat type
hab_adj_barplot <- ggplot(data=den_hab_plot, 
  aes(x= mfws_habitat, y=den_perc_adj, fill=mfws_habitat)) +
  geom_col() +
  annotate("text", x=4.4, y=60, label = parse(text = "bold('b')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"),
        axis.line.y = element_line(colour="grey")) +
  xlab("") + ylab("Percentage adjusted for availability (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=0.1, end=0.9) +
  scale_y_continuous(limits=c(0, 60), breaks=seq(0, 60, 10)) 

# Display the plot
print(hab_adj_barplot)

# Combine the plots
hab_barplots <- ggarrange(hab_raw_barplot, hab_adj_barplot, ncol= 2, nrow= 1, 
                                   widths= c(1,1), heights = c(1,1)) %>%
  annotate_figure(bottom=text_grob("Habitat type", hjust=0.5, vjust=0))


# Display the combined plot
print(hab_barplots)
```

```{r, eval=FALSE, include=FALSE}
ggsave(hab_barplots, 
       filename="output/habitat type barplots.jpeg", 
       width=225, height=100, units="mm")
```

## Activity 

  a) Plots

```{r}
# Transfer from wide to long format, and calculate percent for each level
long <- data %>%
  group_by(activity) %>%
  summarise(percent = 100/nrow(plot_data)*length(activity), 
            n = length(activity)) %>%
    na.omit() %>%
  mutate(activity = factor(activity, levels=c("Active", "Inactive", "Cannot find")))

# Create histogram plot
activity_hist <- ggplot(long, aes(x= activity, y=percent, fill=activity)) +
  geom_col() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"),
        axis.line.y = element_line(colour="grey")) +
  xlab("Activity") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=1, end=0.7) +
  scale_y_continuous(limits=c(0, 100), breaks=seq(0, 100, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Histogram for publication
activity_hist_pub <- ggplot(long, aes(x= activity, y=percent, fill=activity)) +
  geom_col() +
  annotate("text", x=3.5, y=80, label = parse(text = "bold('e')"), size=5) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.border = element_blank(), 
        plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
        legend.position = "none", 
        axis.line.x = element_line(colour="grey"),
        axis.line.y = element_line(colour="grey")) +
  xlab("Activity") + ylab("Percentage (%)") +
  scale_fill_viridis(discrete=T, option="D", begin=0.1, end=0.9) +
  scale_y_continuous(limits=c(0, 80), breaks=seq(0, 80, 10)) +
  geom_text(aes(label=paste("n = ", n), vjust=-0.5))

# Display the plot
print(activity_hist_pub)

# Subset to those records with valid activity types
plot_data <- subset(data, activity!="NA")

# Create map plot
activity_map <- ggplot() + 
  geom_path(mfws, mapping=aes(x=long, y=lat, group=group), 
            col="grey15") +
  geom_point(plot_data, mapping=aes(x=easting, y=northing, 
                                    col=activity), alpha=0.8) +
  coord_sf(xlim=c(695750, 699250), ylim=c(6104250, 6107750)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.70, 0.25),
        legend.title = element_text(size=9), 
        legend.text = element_text(size=9),
        legend.key.size = unit(0.5, 'cm'),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(-0.5, -1, -1, -1), "cm"),
        legend.key = element_blank()) +
  scale_colour_manual(values=viridis(length(unique(data$activity)), 
                                     begin=1, end=0), name="Den activity",
                      breaks=c('Active', 'Inactive', 'Cannot find')) +
  xlab("") + ylab("")

# Combine histogram and map into single plot
activity_plots <- ggarrange(activity_hist, activity_map, 
                            ncol=2, nrow=1, widths=c(3,3))
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(activity_plots, 
       filename="output/den activity histogram and map.jpeg", 
       width=200, height=75, units="mm")

# Export the plot to a jpeg file
ggsave(activity_hist_pub, 
       filename="output/activity histogram publication.jpeg", 
       width=150, height=75, units="mm")
```

  b) Models
  
```{r}
# Chi-square test
activity_mod <- data %>%
  subset(activity!="Cannot find") %>%
  freq_table(activity) %>%
  freq_test() %>%
  select(c(1:3, 5, 7, 12:14))
```

## Summary plot for publication

```{r}
summary_plot_dens <- ggarrange(den_type_hist_pub, 
                  den_substrate_hist_pub + rremove("ylab"),
                  hjust=-44, vjust=2,
                  ncol=2, nrow=1, widths=c(1,1))

print(summary_plot_dens)
summary_plot_den2 <- ggarrange(den_cover_hist_pub, 
                  hjust=-100, vjust=2, ncol=1, nrow=1, widths=c(1)) 

summary_plot_den3 <- ggarrange(entrances_hist_pub,activity_hist_pub + 
                               rremove("ylab"), 
                               hjust=-60, vjust=2, ncol=2, nrow=1, 
                               widths=c(1,1), heights = 2)

den_summary_plot <- ggarrange(summary_plot_dens,summary_plot_den2,
                               summary_plot_den3,ncol= 1, nrow= 3, 
                               widths= c(2,0.5,1), heights = c(2.5,3,2))
print(den_summary_plot)

# Combine grazing zone and habitat type plots for publication
summary_barplots <- ggarrange(hab_barplots, zone_barplots, ncol= 1, nrow= 2, 
                                   widths= c(1,1), heights = c(1,1)) 

print(summary_barplots)
```

```{r, echo=FALSE, eval=FALSE}
# Export the plot to a jpeg file
ggsave(summary_plot_dens, filename="output/summary plot dens.jpeg", 
       width=275, height=225, units="mm")

ggsave(summary_plot_den3, filename="output/summary plot dens3.jpeg", 
       width=275, height=225, units="mm")

ggsave(den_summary_plot, filename="output/den summary plot.jpeg", 
       width=250, height=225, units="mm")

ggsave(summary_barplots, filename="output/grazing zone habitat summary plot.jpeg", 
       width=225, height=175, units="mm")
```

# **Session information**

```{r}
sessionInfo()
```