analyses.Rmd

---
title: "Inducing conditioned taste aversion in red foxes to bush stone-curlew eggs using synthetic scent"
author: "Ellie Lambden, Belinda Wilson"
date: "May 6, 2024"
output:
  html_document:
    toc: true
    number_sections: true
    toc_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=TRUE, warning = FALSE, message = FALSE)
```

# Background

Major questions:

  1. Did we induce CTA in red foxes?
  2. How long did it take to induce CTA?
  3. Did we induce CTA in ravens?

# Data preparation

We used the pacman [Package Management Tool](https://cran.r-project.org/web/packages/pacman/index.html) to install and load packages in a condensed and efficient way.

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

```{r, results='hide', warning=FALSE, message=FALSE}
# Install and load required packages
pacman::p_load(beepr, dplyr, ggmap, ggspatial, ggplot2, ggpubr, janitor, lme4, lmerTest, lwgeom, MuMIn, ncdf4, raster, readxl, sf, sp, tidyr, tidyverse, viridis)

# Assign beepr sound (options are ping, coin, fanfare, complete, treasure, ready, shotgun, mario, wilhelm, facebook, sword)
sound <- "coin"
```

## Bait take

Here we prepare our bait take data to get it ready for analyses and plotting.

```{r}
# Read in data
data <- read_excel(path="input/raw data.xlsx", sheet="data") %>% 
  clean_names() %>% 
  mutate(date = as.Date(date)) %>% 
  mutate(bait_take = as.numeric(bait_take)) %>% 
  mutate(phase = as.factor(phase)) %>% 
  drop_na(bait_take) %>% 
  mutate(bait_yes = ifelse(bait_take == 1, "Yes", "No"), 
         prints_count = as.numeric(prints_count)) %>%
  suppressWarnings() 
```

We exclude Friday's data from the dataset here. The cameras were inconsistent so, if an egg was missing when baits were checked on Monday and the camera had failed, we would not be able to say with certainty that it had been taken within 24 hours of deployment. This means the data from that site would be excluded for that day. If an egg was not taken, we would still be able to use that data. As a result, this skewed our results towards egg rejection. To be consistent, it is better to reject all of Fridays data.

```{r}
# Create a column for day of the week from the date variable
data <- data %>%
  mutate(day_of_week = weekdays(date)) %>%
# Subset to any day except Fridays
  subset(day_of_week != "Friday")
```

## Weather

The weather data were acquired from the [Australian Bureau of Meteorology weather station 70351](http://www.bom.gov.au/climate/dwo/202402/html/IDCJDW2801.202402.shtml).

```{r}
# Read in weather data
weather <- read_excel(path="input/raw data.xlsx", sheet="weather") %>%
  clean_names() %>%
  mutate(date = as.Date(date))

# Combine bait take and weather data into a single df
data <- left_join(data, weather, by=c("date"))
```

## Visits

```{r}
# Subset to actual instances of bait take by a vertebrate
data <- data %>%
  subset(species != "Lizard") %>%
  # Change columns to numeric
  mutate_at(7:22, as.numeric) %>%
  filter(!is.na(fox_visit_1_pa),
         !is.na(fox_visit_2_pa),
         !is.na(fox_visit_3_pa), 
         !is.na(fox_visit_4_pa), 
         !is.na(fox_visit_5_pa), 
         
         !is.na(raven_visit_1_pa), 
         !is.na(raven_visit_2_pa), 
         !is.na(raven_visit_3_pa)) %>%
  
  mutate(fox_visit_pa = ifelse(fox_visit_1_pa !=0, 1, 
                        ifelse(fox_visit_2_pa !=0, 1,
                        ifelse(fox_visit_3_pa !=0, 1, 
                        ifelse(fox_visit_4_pa !=0, 1, 
                        ifelse(fox_visit_5_pa !=0, 1, 0)))))) %>%
  
  mutate(fox_visit_n = ifelse(fox_visit_5_pa !=0, 5, 
                       ifelse(fox_visit_4_pa !=0, 4, 
                       ifelse(fox_visit_3_pa !=0, 3,
                       ifelse(fox_visit_2_pa !=0, 2, 
                       ifelse(fox_visit_1_pa !=0, 1, 0)))))) %>%

  mutate(raven_visit_pa = ifelse(raven_visit_1_pa !=0, 1, 
                          ifelse(raven_visit_2_pa !=0, 1,
                          ifelse(raven_visit_3_pa !=0, 1, 0)))) %>%
  
  mutate(raven_visit_n = ifelse(raven_visit_3_pa !=0, 3, 
                         ifelse(raven_visit_2_pa !=0, 2,
                         ifelse(raven_visit_1_pa !=0, 1, 0))))

# Read in processed data
data <- data %>%
  clean_names() %>%
  mutate(date = as.Date(date)) %>%
  # Create study day variable
  mutate(study_day = yday(date)-14)

# Save the processed data
write.csv(data, "output/data processed.csv")
```

## GIS

Here we read in and reproject polylines representing the anthropogenic boundary (urban and roads) and natural boundary (river corridor) at the Ginninderry Conservation Corridor (GCC) and calculate the distance from each bait station to these boundaries.

What we need to do for each layer is check its inherent Coordinate Reference System (CRS, using `st_transform()` to find the last 4-digit number of its output), remind the layer of that number (using `st_set_crs()`), and then reproject it to the CRS we wish to work with (using `st_transform("EPSG:4326")`). 

### Bait station points

```{r}
# Read in bait station points
stations <- st_read("input/bait_stations_egg_cta.shp") %>%
  as.data.frame() %>%
  fortify() %>%
  clean_names()

# Transform to Spatial Points Dataframe
stations_spdf <- SpatialPointsDataFrame( 
  data.frame(stations$long, stations$lat), 
  stations, proj4string=CRS("EPSG:4326")) %>% 
  st_as_sf()
```

### Anthropogenic boundary

```{r}
# Read in the layer
gcc_anth <- st_read("input/anthropogenic_boundary.shp")

# Remind the layer of its CRS then reproject to our desired CRS
gcc_anth <- gcc_anth %>%
  # Set projection using sf
  st_set_crs(st_crs(gcc_anth)) %>%
  st_transform("EPSG:4326")

gcc_anth_geo <- st_geometry(obj = gcc_anth) %>%
  st_cast(to = 'LINESTRING')

# Combine the linestrings into a single geometry
gcc_anth_com <- st_union(gcc_anth_geo)

# Calculate distance of bait stations to nearest edge of polygon
dist_anth <- gcc_anth_com %>%
  st_distance(y=stations_spdf) %>%
  as.data.frame() %>%
  # Transpose the df from columns to rows
  t()

# Join distances to bait stations df 
stations_spdf <- cbind(stations_spdf, dist_anth)
```

### Natural boundary

```{r}
# Read in the layer
gcc_nat <- st_read("input/natural_boundary.shp")

# Remind the layer of its CRS then reproject to our desired CRS
gcc_nat <- gcc_nat %>%
  # Set projection using sf
  st_set_crs(st_crs(gcc_nat)) %>%
  st_transform("EPSG:4326")

# Calculate distance of bait stations to nearest edge of polygon
dist_nat <- st_geometry(obj = gcc_nat) %>%
  st_cast(to = 'MULTILINESTRING') %>%
  st_distance(y=stations_spdf) %>%
  as.data.frame() %>%
  # Transpose the df from columns to rows
  t()

# Join distances to bait stations df 
stations_spdf <- cbind(stations_spdf, dist_nat)
```

Here we combine the GIS dataframe with the previously prepared data.

```{r}
# Rename name to site
stations <- stations_spdf %>%
  dplyr::rename(site = name)

# Join response df to spatial df
data_sp <- left_join(data, stations, by=c("site", "treatment")) %>%
  filter(!is.na(bait_take))

# Subset to only variables in the subsequent models and plots
data_mod <- data_sp %>%
  dplyr::select(treatment, site, date, study_day, 
                bait_take, fox_visit_n, raven_visit_n, prints_count, 
                maximum_temperature_c, minimum_temperature_c, 
                rainfall_mm, speed_of_maximum_wind_gust_km_h, 
                dist_anth, dist_nat, elevation) %>%
  na.omit()
```

# Models

Here we model the relationship between our response variables (i.e., `bait_take`, `fox_vists`,  `raven_visit_n`, and `prints_count`) and:

  - Treatment (treatment or control, `treatment`)
  - Study day (1-23, `study_day`)
  - Maximum temperature (`maximum_temperature_c`)
  - Minimum temperature (`minimum_temperature_c`)
  - Rainfall (mm, `rainfall_mm`)
  - Maximum wind gust (km/hr, ``)
  - Distance to anthropogenic boundary (urban area and roads) to the east (`gcc_anth`) 
  - Distance to natural boundary (river corridor) to the west (`gcc_nat`)
  - Elevation (m, `elevation`)


## Responses

Here we pose the methodological question, are our response variables (i.e., `bait_take`, `fox_vists`,  `raven_visit_n`, and `prints_count`) are affected by each other.

```{r}
# Fit a GLM
summary(glm(bait_take ~ fox_visit_n + raven_visit_n + 
              prints_count, family=gaussian, data=data_mod))
```

## Correlation tests

Before running GLMs with our continuous predictors, we generate a correlation matrix (values are R-squared) to check whether any are significantly correlated with each other (i.e., R^2 above 0.7 is considered significant). 

```{r}
# Subset data for correlation test
data_corr <- data_mod %>%
  dplyr::select(maximum_temperature_c, minimum_temperature_c, 
                rainfall_mm, speed_of_maximum_wind_gust_km_h, 
                dist_nat, dist_anth, elevation)

# Generate correlation matrix
cor(data_corr)
```

Because `dist_nat` and `elevation` are significantly correlated, we drop `elevation` from our models.

## Bait take

```{r}
# Experimental and climatic model
summary(glm(bait_take ~ treatment + study_day + rainfall_mm + 
             maximum_temperature_c + minimum_temperature_c + 
             speed_of_maximum_wind_gust_km_h, 
           family=gaussian, data=data_mod))

# Geographic model
summary(glm(bait_take ~ dist_anth + dist_nat, 
            family=gaussian, data=data_mod))

# Competing model
mod <- glm(bait_take ~ treatment + dist_anth + dist_nat, 
            data=data_mod, na.action = "na.fail")

# Display summary statistics
summary(mod)

# Perform model selection using dredge() with AICc
dredge(mod, rank = "AICc") %>%
  filter(delta <2)
```

Bait take was significantly influenced by treatment, distance to anthropogenic boundary, and distance to natural boundary.

## Fox visits

```{r}
# Experimental and climatic model
summary(glm(fox_visit_n ~ treatment + study_day + rainfall_mm + 
              maximum_temperature_c + minimum_temperature_c + 
              speed_of_maximum_wind_gust_km_h, 
            family=gaussian, data=data_mod))

# Geographic model
summary(glm(fox_visit_n ~ dist_anth + dist_nat, 
            family=gaussian, data=data_mod))

# Prints model
summary(glm(fox_visit_n ~ prints_count, 
            family=gaussian, data=data_mod))
```

Fox visits was significantly influenced by distance to natural boundary.

## Raven visits

```{r}
# Experimental and climatic model
summary(glm(raven_visit_n ~ treatment + study_day + rainfall_mm + 
              maximum_temperature_c + minimum_temperature_c + 
              speed_of_maximum_wind_gust_km_h, 
            family=gaussian, data=data_mod))

# Geographic model
summary(lm(raven_visit_n ~ dist_anth + dist_nat, data=data_mod))
```

Raven visits was not significantly influenced by any of our predictors.

# Plots

## Bait take

Percentage bait taken by each species by phase

```{r}
# Subset to only potential foxes
taken <- data %>%
  subset(bait_take != 0) %>%
  subset(species != "Not taken") %>%
  group_by(species, treatment, phase) %>%
  summarise(perc = n() / nrow(data) * 100)

# Plot by treatment type
taken_species_phase_hist <- ggplot(data=taken, aes(x = species, 
                                                   y = perc, 
                                                   col=treatment,
                                                   fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("All species") +
  xlab("Species") + 
  ylab("Percentage of taken baits (%)") + 
  labs(fill="Treatment", col="Treatment") +
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1)) + 
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1))

# Display the plot
print(taken_species_phase_hist)
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave("output/taken_species_phase_hist.jpeg", 
       height=1500, width=2000, unit="px")
```

```{r}
# Plot bait take by treatment across time
bait_treat_day_all_sp <- ggplot(data_mod, aes(x = study_day, 
                                       y = bait_take, 
                                       col=treatment, 
                                       fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position="right") +
  xlab("Study Day") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment") +
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1), name=NULL) +
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1), name=NULL) 

# Plot bait take by treatment
bait_treat <- ggplot(data_mod, aes(x = treatment,
                                   y = bait_take, 
                                   col=treatment, 
                                   fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("none") +
  xlab("Treatment") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment") +
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1)) +
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1)) 

# Display plot
print(bait_treat_day_all_sp)

ggsave("output/bait_treat_day_all_sp.jpeg", 
       height=1000, width=2000, unit="px")
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave(plot=bait_treat_both, "output/bait take by treatment.jpeg", 
       height=50, width=100, unit="mm")
```


### Influences
```{r}
# Bait take plot
bait_take_prints <- ggplot(data_mod, aes(x = prints_count, 
                                       y = bait_take, 
                                       col=treatment, 
                                       fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position="none") +
  ggtitle("Bait-take by prints") +
  xlab("Prints") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment") +
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1)) +
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1)) 

# Display the plot
print(bait_take_prints)

# Fox visits plot
bait_take_fox_visits <- ggplot(data_mod, aes(x = fox_visit_n, 
                                             y = bait_take, 
                                             col = treatment, 
                                             fill = treatment)) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position = "none") +
  ggtitle("Bait-take by Fox Visits") +
  xlab("Fox Visits") + 
  ylab("Bait Take") + 
  labs(fill = "Treatment", col = "Treatment") +
  scale_colour_manual(values = viridis(2, begin = 0.8, end = 0.1)) +
  scale_fill_manual(values = viridis(2, begin = 0.8, end = 0.1)) +
  facet_wrap(~ treatment)

# Display plot
print(bait_take_fox_visits)

# Raven visits plot
bait_take_raven_visits <- ggplot(data_mod, aes(x = raven_visit_n, 
                                               y = bait_take, 
                                               col = treatment, 
                                               fill = treatment)) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position = "none") +
  ggtitle("Bait-take by Raven Visits") +
  xlab("Raven Visits") + 
  ylab("Bait Take") + 
  labs(fill = "Treatment", col = "Treatment") +
  scale_colour_manual(values = viridis(2, begin = 0.8, end = 0.1)) +
  scale_fill_manual(values = viridis(2, begin = 0.8, end = 0.1)) +
  facet_wrap(~ treatment)

# Display plot
print(bait_take_raven_visits)
```

## Fox visits

Percentage fox visits by treatment and phase

```{r}
# Subset to only potential foxes
fox_stats <- data %>%
  subset(species == "Fox") %>%
  group_by(treatment) %>%
  summarise(perc = n() / nrow(data) * 100)
  
# Plot by treatment type
fox_vis_hist <- ggplot(data=fox_stats, 
                       aes(y = perc, x = treatment, 
                           col=treatment, fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("Potential fox visits (including unknowns)") +
  xlab("Species") + 
  ylab("Percentage of fox visits (%)") + 
  labs(fill="Treatment", col="Treatment") +
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1)) +
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1)) 

# Display the plot
print(fox_vis_hist)
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave(plot=fox_vis_hist, "output/fox visits histogram.jpeg", 
       height=150, width=300, unit="mm")
```

```{r}
# Subset to only potential foxes
foxes_only <- data %>%
  subset(species != "Raven") %>%
  subset(species != "Unknown")
  
# Plot by treatment type
bait_treat_phase <- ggplot(data=foxes_only, aes(x=date, 
                                                y=bait_take, 
                                                col=treatment, 
                                                fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("Foxes vs no take (excluding unknowns)") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment") +
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1)) +
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1)) 

# Display the plot
print(bait_treat_phase)
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave("output/fox take by treatment.jpeg", 
       height=1500, width=3000, unit="px")
```

Prints as a representation of fox visits 

```{r}
# Plot by treatment type
visit_print <- ggplot(data, aes(x = prints_count, 
                                y = fox_visit_n)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  xlab("Print Count") + 
  ylab("Fox visits") +
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1)) +
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1)) 

# Display the plot
print(visit_print)
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave("output/fox visits by print count.jpeg", 
       height=1000, width=2000, unit="px")
```

```{r}
# Plot the distances against bait take
dist_anth_plot <- ggplot(data = data_mod) +
  geom_smooth(aes(x=dist_anth, y=bait_take, 
                  col=treatment, fill=treatment)) + 
  xlab("Distance to nearest part of anthropogenic boundary (m)") + 
  ylab("Bait take (proportion)") + 
  theme_minimal() + 
  theme(legend.position="none") + 
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1)) +
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1)) +
  ylim(0,1)

dist_nat_plot <- ggplot(data = data_mod) +
  geom_smooth(aes(x=dist_nat, y=bait_take, 
                  col=treatment, fill=treatment)) + 
  xlab("Distance to nearest part of natural boundary (m)") + 
  ylab("Bait take (proportion)") + 
  theme_minimal() +
  scale_colour_manual(values=viridis(2, begin=0.8, end=0.1), name=NULL) +
  scale_fill_manual(values=viridis(2, begin=0.8, end=0.1), name=NULL) +
  ylim(0,1)

dist_plots <- ggarrange(dist_anth_plot, dist_nat_plot, ncol=2, labels=c("(a)", "(b)"), label.x=0.05, label.y=0.9,font.label=list(size=10))

# Display the plot
print(dist_plots)
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave(plot=dist_plots, 
       "output/bait take by treatment and distances.jpeg", 
       height=3, width=9, dpi=300)
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave(plot=mod_map, "output/fox visits map.jpeg", 
       height=4, width=6, dpi=300)
```

## Raven visits

Percentage raven visits by treatment and phase

```{r}
# Subset to only potential ravens
raven_stats <- data %>%
  subset(species == "Raven") %>%
  group_by(treatment, phase) %>%
  summarise(perc = n() / nrow(data) * 100)
  
# Plot by treatment type
visits_raven_hist <- ggplot(data=raven_stats, 
                              aes(y = perc, x = treatment, 
                                  col=treatment, fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  facet_wrap(~phase) +
  ggtitle("Potential raven visits (including unknowns)") +
  xlab("Species") + 
  ylab("Percentage of raven visits (%)") + 
  labs(fill="Treatment", col="Treatment") +
  scale_colour_manual(values=viridis(2, begin=0.8, 
                                    end=0.1), name=NULL) +    
  scale_fill_manual(values=viridis(2, begin=0.8, 
                                   end=0.1), name=NULL)

# Display the plot
print(visits_raven_hist)
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave("output/visits_raven_hist.jpeg", 
       height=1500, width=3000, unit="px")
```

```{r}
# Subset to only potential ravens
ravens_only <- data %>%
  subset(species != "Fox") %>%
  subset(species != "Unknown")
  
# Plot by treatment type
bait_raven_treat_phase <- ggplot(data=ravens_only, 
                                 aes(x = as.Date(date), 
                                     y = bait_take, 
                                     col=treatment, 
                                     fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("Ravens vs no take (excluding unknowns)") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment") +
  scale_colour_manual(values=viridis(2, begin=0.8, 
                                     end=0.1), name=NULL) +    
  scale_fill_manual(values=viridis(2, begin=0.8, 
                                   end=0.1), name=NULL)

# Display the plot
print(bait_raven_treat_phase)
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave("output/bait_ravens_treat_phase.jpeg", 
       height=1500, width=3000, unit="px")
```

```{r}
# Model for raven visits by study day
summary(glm(raven_visit_n ~ treatment + study_day, 
            data=data))
```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave(plot=mod_map, "output/raven visits map.jpeg", 
       height=4, width=6, dpi=300)
```

## Bait take by visits

```{r}
# Plot bait take by species visits

bait_visits <- ggplot(data, mapping=aes(x = study_day)) +
  geom_smooth(data, mapping=aes(y = as.numeric(bait_take), 
                                color = "Bait Take", 
                                fill = "Bait Take")) +
  geom_smooth(data, mapping=aes(y = as.numeric(fox_visit_n), 
                                color = "Fox Visits", 
                                fill = "Fox Visits"))  +
  geom_smooth(data, mapping=aes(y = as.numeric(raven_visit_n), 
                                color = "Raven Visits", 
                                fill = "Raven Visits"))  +
  scale_color_manual(values = c("Bait Take" = "red", 
                                "Fox Visits" = "darkorange", 
                                "Raven Visits" = "black"),
                     labels = c("Bait Take", 
                                "Fox Visits", 
                                "Raven Visits"), name = NULL) +
  scale_fill_manual(values = c("Bait Take" = "red", 
                               "Fox Visits" = "darkorange", 
                               "Raven Visits" = "black"),
                    labels = c("Bait Take", "Fox Visits", 
                               "Raven Visits"), name = NULL) +
  theme_minimal() + 
  facet_wrap(~treatment) +
  scale_y_continuous(name = "Bait take (%)", limits = c(0, 1), 
                     sec.axis = sec_axis(~., 
                                         name = "Fox and raven visits", 
                                         breaks = c(0, 1))) +
  xlab("Study Day")

# Display the plot
print(bait_visits)

```

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave("output/bait_visits.jpeg", 
       height=800, width=2000, unit="px")
```

## Maps

```{r, eval=FALSE, include=FALSE}
# Save the plot as a jpeg
ggsave(plot=gcc_map, "output/gcc features and bait stations map.jpeg", 
       height=4, width=6, dpi=300)
```

Satellite map

```{r}
# Use Google API to fetch a base map
# ggmap::register_google(key="[add API key here]", write=TRUE)
map <- get_map(location=c(lon=148.9811, lat=-35.2350),
               zoom=15, source="google", maptype="satellite", crop=FALSE)
```

```{r}
# Read in shapefile
gis <- st_read("input/bait_stations_egg_cta.shp") %>%
  fortify() %>%
  clean_names() %>%
  rename(site=name)

# Map the bait stations
gin_map <- ggmap(map) + 
  geom_point(data=gis, size=2, 
             aes(x=long, y=lat, col=treatment)) + 
  theme_minimal() +
  xlab("") + ylab("") + labs(col="Treatment")

# Display the plot
print(gin_map)
```

# Session information

```{r}
# Display version information for R, OS, and packages
sessionInfo()
```