# INPUT: - Cleaned, combined arthropod ID and measuring data: combined_id_meas.csv

# OUTPUT: - JWM manuscript Figure 3: Arthropod abundance, biomass, and family 
# richness measures (untransformed data) through time at reference and focal 
# study sites. Error bars represent standard errors. Asterisks indicate significant 
# differences based on model results (step 2; βsite type × timing).

# Load libraries
library(dplyr)
library(tidyr)
library(ggplot2)
library(forcats)
library(tidyselect)
library(ggpubr)
theme_set(theme_pubr())

# Read in insect data
dat = read.csv("data/20230303_combined_id_meas_manually_cleaned.csv")


################# Calculate arthropod biomass ##################################

##### Functions: biomass by order - Straus & Aviles 2018

# Coleoptera
bio.col <- function(x) {
  mg.col = 10^(-1.960 + 2.966*(log10(x)))
  return(mg.col)
}

# Orthoptera
bio.ort <- function(x) {
  mg.ort = 10^(-1.928 + 2.680*(log10(x)))
  return(mg.ort)
}

# Hymenoptera
bio.hym <- function(x) {
  mg.hym = 10^(-0.611 + 1.141*(log10(x)))
  return(mg.hym)
}

# Lepidoptera
bio.lep <- function(x) {
  mg.lep = 10^(-3.169 + 3.624*(log10(x)))
  return(mg.lep)
}

# Hemiptera
bio.hem <- function(x) {
  mg.hem = 10^(-1.749 + 2.529*(log10(x)))
  return(mg.hem)
}

# Diptera
bio.dip <- function(x) {
  mg.dip = 10^(-1.361 + 2.026*(log10(x)))
  return(mg.dip)
}

# Araneae
bio.ara <- function(x) {
  mg.ara = 10^(-2.643 + 4.011*(log10(x)))
  return(mg.ara)
}

# All groups
bio.all <- function(x) {
  mg.all = 10^(-1.412 + 2.085*(log10(x)))
  return(mg.all)
}

# Test functions
y = 20
bio.ort(y)
bio.col(y)
bio.hym(y)
bio.lep(y)
bio.hem(y)
bio.dip(y)
bio.ara(y)
bio.all(y)

# Make all columns character
dat[] = lapply(dat, as.character)

# Adjust data types
dat = dat %>%
  mutate_at(c("key",
              "code",
              "year",
              "site", 
              "timing",
              "transect",
              "distance",
              "ID",
              "ID_family",
              "ID_order"), 
            as.factor) %>% # columns to be factors
  mutate_at(c("mm_1",
              "mm_2",
              "mm_3",
              "mm_4",
              "mm_5",
              "mm_6",
              "mm_7",
              "mm_8",
              "mm_9",
              "mm_10",
              "ct_1",
              "ct_2",
              "ct_3",
              "ct_4",
              "ct_5",
              "ct_6",
              "ct_7",
              "ct_8",
              "ct_9",
              "ct_10"), 
            as.numeric) # Columns to be numeric. Warnings: "." changed to NA.

# Make family & order "." = NA
dat$ID_family[dat$ID_family == "."] = NA
dat$ID_order[dat$ID_order == "."] = NA

# Biomass Coleoptera
dat$bio_1= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_1), NA)
dat$bio_2= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_2), NA)
dat$bio_3= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_3), NA)
dat$bio_4= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_4), NA)
dat$bio_5= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_5), NA)
dat$bio_6= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_6), NA)
dat$bio_7= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_7), NA)
dat$bio_8= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_8), NA)
dat$bio_9= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_9), NA)
dat$bio_10= ifelse(dat$ID_order == "Coleoptera", bio.col(dat$mm_10), NA)

# Biomass Orthoptera
dat$bio_1= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(dat$ID_order == "Orthoptera", bio.ort(dat$mm_10), dat$bio_10)

# Biomass Hymenoptera
dat$bio_1= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(dat$ID_order == "Hymenoptera", bio.hym(dat$mm_10), dat$bio_10)

# Biomass Lepidoptera
dat$bio_1= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(dat$ID_order == "Lepidoptera", bio.lep(dat$mm_10), dat$bio_10)

# Biomass Hemiptera
dat$bio_1= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(dat$ID_order == "Hemiptera", bio.hem(dat$mm_10), dat$bio_10)

# Biomass Diptera
dat$bio_1= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(dat$ID_order == "Diptera", bio.dip(dat$mm_10), dat$bio_10)

# Biomass Araneae
dat$bio_1= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(dat$ID_order == "Araneae", bio.ara(dat$mm_10), dat$bio_10)

# Biomass other: for orders other than Coleoptera, Orthoptera, Hymenoptera,
# Lepidoptera, Hemiptera, Diptera, or Araneae
dat$bio_1= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(dat$ID_order != "Coleoptera" & dat$ID_order != "Orthoptera" & dat$ID_order != "Hymenoptera" & dat$ID_order != "Lepidoptera" & dat$ID_order != "Hemiptera" & dat$ID_order != "Diptera" & dat$ID_order != "Araneae", bio.all(dat$mm_10), dat$bio_10)

# Biomass other: for order = Unknown
dat$bio_1= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(dat$ID_order == "Unknown", bio.all(dat$mm_10), dat$bio_10)

# Biomass other: for order = NA
dat$bio_1= ifelse(is.na(dat$ID_order), bio.all(dat$mm_1), dat$bio_1)
dat$bio_2= ifelse(is.na(dat$ID_order), bio.all(dat$mm_2), dat$bio_2)
dat$bio_3= ifelse(is.na(dat$ID_order), bio.all(dat$mm_3), dat$bio_3)
dat$bio_4= ifelse(is.na(dat$ID_order), bio.all(dat$mm_4), dat$bio_4)
dat$bio_5= ifelse(is.na(dat$ID_order), bio.all(dat$mm_5), dat$bio_5)
dat$bio_6= ifelse(is.na(dat$ID_order), bio.all(dat$mm_6), dat$bio_6)
dat$bio_7= ifelse(is.na(dat$ID_order), bio.all(dat$mm_7), dat$bio_7)
dat$bio_8= ifelse(is.na(dat$ID_order), bio.all(dat$mm_8), dat$bio_8)
dat$bio_9= ifelse(is.na(dat$ID_order), bio.all(dat$mm_9), dat$bio_9)
dat$bio_10= ifelse(is.na(dat$ID_order), bio.all(dat$mm_10), dat$bio_10)

# Make columns numeric
dat = dat %>%
  mutate_at(c("bio_1",
              "bio_2",
              "bio_3",
              "bio_4",
              "bio_5",
              "bio_6",
              "bio_7",
              "bio_8",
              "bio_9",
              "bio_10"), 
            as.numeric)

# Multiply each biomass by count into new column mbio_#
dat = dat %>% 
  mutate(mbio_1 = bio_1*ct_1,
         mbio_2 = bio_2*ct_2,
         mbio_3 = bio_3*ct_3,
         mbio_4 = bio_4*ct_4,
         mbio_5 = bio_5*ct_5,
         mbio_6 = bio_6*ct_6,
         mbio_7 = bio_7*ct_7,
         mbio_8 = bio_8*ct_8,
         mbio_9 = bio_9*ct_9,
         mbio_10 = bio_10*ct_10)

# New column mbio_sum: sum of mbio_1-10
dat = dat %>% 
  mutate(mbio_sum = rowSums(dat[,c("mbio_1", "mbio_2", "mbio_3", "mbio_4", "mbio_5", "mbio_6", "mbio_7", "mbio_8", "mbio_9", "mbio_10")], 
                            na.rm = TRUE)) # exclude NAs

# New column count_sum: sum of ct_1-10
dat = dat %>% 
  mutate(count_sum = rowSums(dat[,c("ct_1", "ct_2", "ct_3", "ct_4", "ct_5", "ct_6", "ct_7", "ct_8", "ct_9", "ct_10")], 
                             na.rm = TRUE)) # exclude NAs
# Ideally, column "count_sum" should = column "count"

######################## Count number of vials in dataset #################################

dat %>% 
  distinct(site,timing,transect,distance,ID) %>% 
  tally() # 7640 vials

######################### Count vials with no calcs #########################################

# Count vials with no mbio_sum measurement
dat %>% 
  filter(mbio_sum == 0) %>% 
  tally() # 3 vials

# Count vials with no count_sum measurement
dat %>% 
  filter(count_sum == 0) %>% 
  tally () # 3 vials

##################### Summarize data by vials #############################################

# Each vial own row
vialsum = dat %>% 
  group_by(site, timing, transect, distance, ID, ID_family, ID_order, ID_notes) %>% 
  summarize(totalbio = sum(mbio_sum, na.rm = TRUE),
            totalct = sum(count_sum, na.rm = TRUE))

# Add column of site type (control/treatment)
vialsum$type = ifelse(vialsum$site == "HL", "Treatment", "Control")
vialsum$type = ifelse(vialsum$site == "LM", "Treatment", vialsum$type)
vialsum$type = ifelse(vialsum$site == "MS", "Treatment", vialsum$type)
vialsum$type = ifelse(vialsum$site == "SM", "Treatment", vialsum$type)
vialsum$type = ifelse(vialsum$site == "WW", "Treatment", vialsum$type)


############ Add Larvae/Adult to Lepidoptera ###############

# New column larva: "Y" if ID_notes contain words larva/larvae
vialsum = vialsum %>% 
  ungroup() %>% 
  mutate(larva = ifelse(grepl("larva|larvae", vialsum$ID_notes, ignore.case = TRUE), "Y", "N"))

# New column ID_order.new with new "Lepidoptera Larva" value
vialsum = vialsum %>% 
  mutate(ID_order.new = ifelse(ID_order == "Lepidoptera" & larva == "Y", "Lepidoptera Larva", as.character(ID_order)))

# Have vialsum only include new ID_order column
vialsum = vialsum %>% 
  dplyr::select(-ID_order, -ID_notes)

# Rename ID_order.new to ID_order
names(vialsum)[names(vialsum) == "ID_order.new"] <- "ID_order"

# Rename "Lepidoptera" to "Lepidoptera Adult"
vialsum$ID_order[vialsum$ID_order == 'Lepidoptera'] <- 'Lepidoptera Adult'



################################## Additional cleaning ###########################

# Exclude 1 observation with timing = "."
vialsum = vialsum %>% 
  filter(timing != ".")


# Exclude 1 observation from HL P T 25 - did not collect on transect T pre-spraying
vialsum = vialsum %>% 
  filter(!(timing == "P" & transect == "T"))


############################## Group ############################################

# Group data by site/timing/transect/distance/order/family/type - some arthropods with same family & order are split into separate vials
ins = vialsum %>% 
  group_by(site, timing, transect, distance, ID_order, ID_family, type) %>% 
  summarize(bio = sum(totalbio, na.rm = TRUE),
            ct = sum(totalct, na.rm = TRUE))


############ Correct for HL pre-spraying transects being longer than all other transects ######################

# HL pre-spraying transects were 60 m long while all other transects were 40 m long

# Multiply bio & ct by (2/3) for all HL pre-spraying samples
ins = ins %>% 
  mutate(bio = ifelse(site == "HL" & timing == "P", bio*(2/3), bio),
         ct = ifelse(site == "HL" & timing == "P", ct*(2/3), ct))



########### Address NAs #################

# Change order & family to character
ins = ins %>%
  mutate_at(c("ID_family",
              "ID_order"), 
            as.character)

# Change order & family "NA" to "Unidentified"
ins = ins %>%
  replace_na(list(ID_family = 'Unidentified', ID_order = 'Unidentified')) 



############ Model response variables over time at trt/control sites - scatterplots with error bars ##############

#### Response variables: total abundance, total biomas, bird prey abundance, bird prey biomass, family richness (Araneae, Coleoptera, Hemiptera, & Orthoptera) ####

### Total abundance and biomass at treatment and control sites at 0 & 25 m, primary
### (X/Y/Z) transects

# Subset data
dat.total = ins %>% 
  filter(transect != "T") %>% 
  filter(distance == 0 | distance == 25) %>% 
  group_by(site, type, timing, transect, distance) %>% 
  summarize(ct.sum= sum(ct),
         bio.sum = sum(bio)) %>% 
  arrange(factor(timing, levels = c("P", "3", "20")))

# Summary stats
dat.summary.total = dat.total %>% 
  group_by(type, timing) %>%
  summarize(bio.mean = mean(bio.sum, na.rm = TRUE),
           bio.sd = sd(bio.sum, na.rm = TRUE),
           ct.mean = mean(ct.sum, na.rm = TRUE),
           ct.sd = sd(ct.sum, na.rm = TRUE),
           ss = n(),
           bio.se= bio.sd/(sqrt(ss)),
           ct.se = ct.sd/sqrt(ss)) %>% 
  arrange(factor(timing, levels = c("P", "3", "20")))
  

### Bird prey abundance and biomass at treatment and control sites at 0 & 25 m
# Subset data
dat.bp = ins %>% 
  filter((ID_order == "Coleoptera" | ID_order == "Lepidoptera Larva" | ID_order == "Orthoptera" | ID_order == "Araneae") &
           transect != "T") %>%  
  filter(distance == 0 | distance == 25) %>% 
  group_by(site, type, timing, transect, distance) %>%
  summarize(ct.sum= sum(ct),
            bio.sum = sum(bio)) %>% 
  arrange(factor(timing, levels = c("P", "3", "20")))

# Summary stats
dat.summary.bp = dat.bp %>% 
  group_by(type, timing) %>%
  summarize(bio.bp.mean = mean(bio.sum, na.rm = TRUE),
            bio.bp.sd = sd(bio.sum, na.rm = TRUE),
            ct.bp.mean = mean(ct.sum, na.rm = TRUE),
            ct.bp.sd = sd(ct.sum, na.rm = TRUE),
            ss.bp = n(),
            bio.bp.se= bio.bp.sd/(sqrt(ss.bp)),
            ct.bp.se = ct.bp.sd/sqrt(ss.bp)) %>% 
  arrange(factor(timing, levels = c("P", "3", "20")))

### Araneae family richness at treatment and control sites at 0 & 25 m
# Subset data
dat.ara = ins %>% 
  filter((ID_order == "Araneae") &
           transect != "T" & ID_family != "(Missing)") %>% 
  filter(distance == 0 | distance == 25) %>% 
  group_by(site, type, timing, transect, distance) %>% 
  summarize(fam = length(unique(ID_family)))

# Summary stats
dat.summary.ara = dat.ara %>% 
  group_by(type, timing) %>%
  summarize(fam.ara.mean = mean(fam, na.rm = TRUE),
            fam.ara.sd = sd(fam, na.rm = TRUE),
            ss.ara = n(),
            fam.ara.se= fam.ara.sd/(sqrt(ss.ara))) %>% 
  arrange(factor(timing, levels = c("P", "3", "20")))

### Coleoptera family richness at treatment and control sites at 0 & 25 m
# Subset data
dat.col = ins %>% 
  filter((ID_order == "Coleoptera") &
           transect != "T" & ID_family != "(Missing)") %>%
  filter(distance == 0 | distance == 25) %>% 
  group_by(site, type, timing, transect, distance) %>% 
  summarize(fam = length(unique(ID_family)))

# Summary stats
dat.summary.col = dat.col %>% 
  group_by(type, timing) %>%
  summarize(fam.col.mean = mean(fam, na.rm = TRUE),
            fam.col.sd = sd(fam, na.rm = TRUE),
            ss.col = n(),
            fam.col.se= fam.col.sd/(sqrt(ss.col))) %>% 
  arrange(factor(timing, levels = c("P", "3", "20")))

### Hemiptera family richness at treatment and control sites at 0 & 25 m
# Subset data
dat.hem = ins %>% 
  filter((ID_order == "Hemiptera") &
           transect != "T" & ID_family != "(Missing)") %>%
  filter(distance == 0 | distance == 25) %>% 
  group_by(site, type, timing, transect, distance) %>% 
  summarize(fam = length(unique(ID_family)))

# Summary stats
dat.summary.hem = dat.hem %>% 
  group_by(type, timing) %>%
  summarize(fam.hem.mean = mean(fam, na.rm = TRUE),
            fam.hem.sd = sd(fam, na.rm = TRUE),
            ss.hem = n(),
            fam.hem.se= fam.hem.sd/(sqrt(ss.hem))) %>% 
  arrange(factor(timing, levels = c("P", "3", "20")))

### Orthoptera family richness at treatment and control sites at 0 & 25 m
# Subset data
dat.ort = ins %>% 
  filter((ID_order == "Orthoptera") &
           transect != "T" & ID_family != "(Missing)") %>% 
  filter(distance == 0 | distance == 25) %>% 
  group_by(site, type, timing, transect, distance) %>% 
  summarize(fam = length(unique(ID_family)))

# Summary stats
dat.summary.ort = dat.ort %>% 
  group_by(type, timing) %>%
  summarize(fam.ort.mean = mean(fam, na.rm = TRUE),
            fam.ort.sd = sd(fam, na.rm = TRUE),
            ss.ort = n(),
            fam.ort.se= fam.ort.sd/(sqrt(ss.ort))) %>% 
  arrange(factor(timing, levels = c("P", "3", "20")))


### Subset data to combine into one dataframe

dat.combine.total = dat.summary.total %>% 
  select(type, timing, bio.mean, bio.se, ct.mean, ct.se)

dat.combine.bp = dat.summary.bp %>% 
  select(type, timing, bio.bp.mean, bio.bp.se, ct.bp.mean, ct.bp.se)

dat.combine.ara = dat.summary.ara %>% 
  select(type, timing, fam.ara.mean, fam.ara.se)

dat.combine.col = dat.summary.col %>% 
  select(type, timing, fam.col.mean, fam.col.se)

dat.combine.hem = dat.summary.hem %>% 
  select(type, timing, fam.hem.mean, fam.hem.se)

dat.combine.ort = dat.summary.ort %>% 
  select(type, timing, fam.ort.mean, fam.ort.se)


### Combine tables

alldat = dat.combine.total %>% 
  left_join(dat.combine.bp,  by = c("type", "timing")) %>% 
  left_join(dat.combine.ara,  by = c("type", "timing")) %>%
  left_join(dat.combine.col,  by = c("type", "timing")) %>%
  left_join(dat.combine.hem,  by = c("type", "timing")) %>%
  left_join(dat.combine.ort,  by = c("type", "timing"))

####################################### Individual plots #######################################

### Plot: total abundance
t.ct = ggplot(alldat) + 
  geom_point(aes(timing, ct.mean, shape = type),
             size = 5, position=position_dodge(width=0.4)) +
  scale_x_discrete(limits=c("P", "3", "20"),
                   labels=c("P" = "Pre-spraying",
                            "3" = "3-5 days \n post-spraying",
                            "20" = "19-21 days \n post-spraying")) + 
  xlab("Timing") + ylab("Mean total abundance\n ") + # labels
  geom_errorbar(aes(timing, ymin = (ct.mean - ct.se), ymax = (ct.mean + ct.se), # xyz error bars
                    width = .2, linetype = type), position=position_dodge(width=0.4)) +
  scale_shape_discrete(name = "Site type", labels = c("Reference", "Focal")) +
  scale_linetype_manual(name = "Site type", labels = c("Reference", "Focal"),
                        values=c("Control"=1,"Treatment"=1)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  theme(text=element_text(size=20),   
        legend.text = element_text(size = 20), 
        legend.title = element_text(size = 20), 
        legend.key.size = unit(2, 'cm')) +
  geom_text(aes(x=2, y=151, label=c("*"), hjust=0.5, vjust=0), size=12) +
  geom_segment(aes(x=1.8, xend=2.2, y=155, yend=155), size=1)

t.ct

### Plot: total biomass
t.bio = ggplot(alldat) + 
  geom_point(aes(timing, bio.mean, shape = type),
      size = 5, position=position_dodge(width=0.4)) +
  scale_x_discrete(limits=c("P", "3", "20"),
                   labels=c("P" = "Pre-spraying",
                            "3" = "3-5 days \n post-spraying",
                            "20" = "19-21 days \n post-spraying")) + 
  xlab("Timing") + ylab("Mean total consumable \ndry biomass (mg)") + # labels
  geom_errorbar(aes(timing, ymin = (bio.mean - bio.se), ymax = (bio.mean + bio.se), # xyz error bars
                    width = .2, linetype = type), position=position_dodge(width=0.4)) +
  scale_shape_discrete(name = "Site type", labels = c("Reference", "Focal")) +
  scale_linetype_manual(name = "Site type", labels = c("Reference", "Focal"),
                        values=c("Control"=1,"Treatment"=1)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  theme(text=element_text(size=20),   
        legend.text = element_text(size = 20), 
        legend.title = element_text(size = 20), 
        legend.key.size = unit(2, 'cm'))

t.bio

### Plot: bird prey abundance
t.ct.bp = ggplot(alldat) + 
  geom_point(aes(timing, ct.bp.mean, shape = type),
             size = 5, position=position_dodge(width=0.4)) +
  scale_x_discrete(limits=c("P", "3", "20"),
                   labels=c("P" = "Pre-spraying",
                            "3" = "3-5 days \n post-spraying",
                            "20" = "19-21 days \n post-spraying")) + 
  xlab("Timing") + ylab("Mean bird prey abundance\n ") + # labels
  geom_errorbar(aes(timing, ymin = (ct.bp.mean - ct.bp.se), ymax = (ct.bp.mean + ct.bp.se), # xyz error bars
                    width = .2, linetype = type), position=position_dodge(width=0.4)) +
  scale_shape_discrete(name = "Site type", labels = c("Reference", "Focal")) +
  scale_linetype_manual(name = "Site type", labels = c("Reference", "Focal"),
                        values=c("Control"=1,"Treatment"=1)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  theme(text=element_text(size=20),   
        legend.text = element_text(size = 20), 
        legend.title = element_text(size = 20), 
        legend.key.size = unit(2, 'cm'))

t.ct.bp

### Plot: bird prey biomass
t.bio.bp = ggplot(alldat) + 
  geom_point(aes(timing, bio.bp.mean, shape = type),
             size = 5, position=position_dodge(width=0.4)) +
  scale_x_discrete(limits=c("P", "3", "20"),
                   labels=c("P" = "Pre-spraying",
                            "3" = "3-5 days \n post-spraying",
                            "20" = "19-21 days \n post-spraying")) + 
  xlab("Timing") + ylab("Mean bird prey consumable \ndry biomass (mg)") + # labels
  geom_errorbar(aes(timing, ymin = (bio.bp.mean - bio.bp.se), ymax = (bio.bp.mean + bio.bp.se), # xyz error bars
                    width = .2, linetype = type), position=position_dodge(width=0.4)) +
  scale_shape_discrete(name = "Site type", labels = c("Reference", "Focal")) +
  scale_linetype_manual(name = "Site type", labels = c("Reference", "Focal"),
                        values=c("Control"=1,"Treatment"=1)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  theme(text=element_text(size=20),   
        legend.text = element_text(size = 20), 
        legend.title = element_text(size = 20), 
        legend.key.size = unit(2, 'cm'))

t.bio.bp

### Plot: Araneae family richness
t.fam.ara = ggplot(alldat) + 
  geom_point(aes(timing, fam.ara.mean, shape = type),
             size = 5, position=position_dodge(width=0.4)) +
  scale_x_discrete(limits=c("P", "3", "20"),
                   labels=c("P" = "Pre-spraying",
                            "3" = "3-5 days \n post-spraying",
                            "20" = "19-21 days \n post-spraying")) + 
  xlab("Timing") + ylab("Araneae mean \nfamily richness") + # labels
  geom_errorbar(aes(timing, ymin = (fam.ara.mean - fam.ara.se), ymax = (fam.ara.mean + fam.ara.se), # xyz error bars
                    width = .2, linetype = type), position=position_dodge(width=0.4)) +
  scale_shape_discrete(name = "Site type", labels = c("Reference", "Focal")) +
  scale_linetype_manual(name = "Site type", labels = c("Reference", "Focal"),
                        values=c("Control"=1,"Treatment"=1)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  theme(text=element_text(size=20),   
        legend.text = element_text(size = 20), 
        legend.title = element_text(size = 20), 
        legend.key.size = unit(2, 'cm'))

t.fam.ara

### Plot: Coleoptera family richness
t.fam.col = ggplot(alldat) + 
  geom_point(aes(timing, fam.col.mean, shape = type),
             size = 5, position=position_dodge(width=0.4)) +
  scale_x_discrete(limits=c("P", "3", "20"),
                   labels=c("P" = "Pre-spraying",
                            "3" = "3-5 days \n post-spraying",
                            "20" = "19-21 days \n post-spraying")) + 
  xlab("Timing") + ylab("Coleoptera mean \nfamily richness") + # label
  geom_errorbar(aes(timing, ymin = (fam.col.mean - fam.col.se), ymax = (fam.col.mean + fam.col.se), # xyz error bars
                    width = .2, linetype = type), position=position_dodge(width=0.4)) +
  scale_shape_discrete(name = "Site type", labels = c("Reference", "Focal")) +
  scale_linetype_manual(name = "Site type", labels = c("Reference", "Focal"),
                        values=c("Control"=1,"Treatment"=1)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  theme(text=element_text(size=20),   
        legend.text = element_text(size = 20), 
        legend.title = element_text(size = 20), 
        legend.key.size = unit(2, 'cm')) +
  geom_text(aes(x=2, y=3.85, label=c("*"), hjust=0.5, vjust=0), size=12) +
  geom_segment(aes(x=1.8, xend=2.2, y=3.89, yend=3.89), size=1)

t.fam.col

### Plot: Hemiptera family richness
t.fam.hem = ggplot(alldat) + 
  geom_point(aes(timing, fam.hem.mean, shape = type),
             size = 5, position=position_dodge(width=0.4)) +
  scale_x_discrete(limits=c("P", "3", "20"),
                   labels=c("P" = "Pre-spraying",
                            "3" = "3-5 days \n post-spraying",
                            "20" = "19-21 days \n post-spraying")) + 
  ylab("Hemiptera mean \nfamily richness") + # labels
  geom_errorbar(aes(timing, ymin = (fam.hem.mean - fam.hem.se), ymax = (fam.hem.mean + fam.hem.se), # xyz error bars
                    width = .2, linetype = type), position=position_dodge(width=0.4)) +
  scale_shape_discrete(name = "Site type", labels = c("Reference", "Focal")) +
  scale_linetype_manual(name = "Site type", labels = c("Reference", "Focal"),
                        values=c("Control"=1,"Treatment"=1)) +
  theme(text=element_text(size=20),   
        legend.text = element_text(size = 20), 
        legend.title = element_text(size = 20), 
        legend.key.size = unit(2, 'cm')) +
  theme(axis.title.x=element_blank())

t.fam.hem

### Plot: Orthoptera family richness
t.fam.ort = ggplot(alldat) + 
  geom_point(aes(timing, fam.ort.mean, shape = type),
             size = 5, position=position_dodge(width=0.4)) +
  scale_x_discrete(limits=c("P", "3", "20"),
                   labels=c("P" = "Pre-spraying",
                            "3" = "3-5 days \n post-spraying",
                            "20" = "19-21 days \n post-spraying")) + 
  ylab("Orthoptera mean \nfamily richness") + # labels
  geom_errorbar(aes(timing, ymin = (fam.ort.mean - fam.ort.se), ymax = (fam.ort.mean + fam.ort.se), # xyz error bars
                    width = .2, linetype = type), position=position_dodge(width=0.4)) +
  scale_shape_discrete(name = "Site type", labels = c("Reference", "Focal")) +
  scale_linetype_manual(name = "Site type", labels = c("Reference", "Focal"),
                        values=c("Control"=1,"Treatment"=1)) +
  theme(text=element_text(size=20),   
        legend.text = element_text(size = 20), 
        legend.title = element_text(size = 20), 
        legend.key.size = unit(2, 'cm')) +
  theme(axis.title.x=element_blank())


t.fam.ort

####################################### Combine plots #######################################

figure = ggarrange(t.ct, t.bio, t.ct.bp, t.bio.bp, t.fam.ara, t.fam.col, t.fam.hem, t.fam.ort,
                    labels = c("(A)", "(B)", "(C)", "(D)", "(E)", "(F)", "(G)", "(H)"),
                    ncol = 2, nrow = 4,
                   common.legend = TRUE, legend = "bottom",
                   font.label = list(size = 24))

figure

#######################################################################################################################

###### Export #######
# Datestamp
datestamp=as.Date(Sys.Date())
datestamp=format(datestamp, format="%Y%m%d")

### Combined figure
jpeg(filename = paste0("figures/", datestamp, "_JWM_abund_bio_richness.jpeg"), # file name
     width = 1000, height = 1400, # size
     quality = 100) # image quality (%)
plot(figure) # what to save
dev.off() # close