library(tidyverse)
library(multcomp)
library(tibble)
library(lubridate)
library(ggplot2)
library(ggpubr)
library(glmmTMB)
library(DHARMa)

# Data preparation ####
carabids <- c("BEMLAM","BEMQUA","HARRUF","PLADOR","POECUP","PTEMEL")
carabids_autumn <- c("BEMLAM","BEMQUA","HARRUF","NEBBRE","PLADOR","POECUP","PTEMEL")
carabid_labels = c("BEMLAM" = "Bembidion\nlampros",
                   "BEMQUA" = "Bembidion\nquadrimaculatum",
                   "HARRUF" = "Harpalus\nrufipes",
                   "NEBBRE" = "Nebria\nbrevicollis",
                   "PLADOR" = "Anchomenus\ndorsalis",
                   "POECUP" = "Poecilus\ncupreus",
                   "PTEMEL" = "Pterostichus\nmelanarius")

breaks_2021 <- as.POSIXct(c("2021-06-07", "2021-06-14","2021-06-21", "2021-06-28","2021-07-05","2021-07-12"))
breaks_2022 <- as.POSIXct(c("2022-06-07", "2022-06-14","2022-06-21", "2022-06-28","2022-07-05","2022-07-12"))


# 2021a (RS / RundS)
c_UKA_RS_21<-read.csv("2021a carabids.csv",sep=";",header=T)
c_UKA_RS_21 <- c_UKA_RS_21 %>%
  mutate(across(c(ROW, DIRECTION, SURROUNDING, BLOCK, TREATMENT), 
                as.factor),
         DATE = as.Date(DATE),
         TOTAL_CARABIDS = rowSums(across(7:last_col()), na.rm = TRUE))

gc_UKA_RS_21 <- c_UKA_RS_21 %>%
  gather("CARABID","COUNT",7:last_col()) %>%
  mutate(CARABID = as.factor(CARABID))

sc_UKA_RS_21 <- gc_UKA_RS_21 %>% 
  filter(DATE < "2021-07-31",
         !(TREATMENT %in% c("WT","S")) | SURROUNDING == "mulch") %>%
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

tc_UKA_RS_21 <- sc_UKA_RS_21 %>% 
  group_by(BLOCK,TREATMENT,CARABID) %>% 
  summarise(sum = sum(t_mean,na.rm = T))

fc_UKA_RS_21 <- tc_UKA_RS_21 %>% 
  filter(CARABID %in% carabids)

ssc_UKA_RS_21 <- fc_UKA_RS_21 %>% 
  group_by(TREATMENT,CARABID) %>% 
  summarise(mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

t_UKA_RS_21 <- sc_UKA_RS_21 %>% 
  group_by(DATE,TREATMENT,CARABID) %>% 
  summarise(mean = mean(t_mean, na.rm = T),
            se = sd(t_mean, na.rm = T)/sqrt(length(t_mean)),
            n = length(t_mean))

tt_UKA_RS_21 <- t_UKA_RS_21 %>% 
  filter(CARABID == "TOTAL_CARABIDS")

a1_UKA_RS_21 <- subset(gc_UKA_RS_21, DATE > "2021-07-31" & DATE < "2021-10-01") %>% 
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

a2_UKA_RS_21 <- subset(gc_UKA_RS_21, DATE > "2021-10-01") %>% 
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

a3_UKA_RS_21 <- subset(gc_UKA_RS_21, DATE > "2021-07-31") %>% 
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

ta1_UKA_RS_21 <- a1_UKA_RS_21 %>% 
  group_by(BLOCK,TREATMENT,CARABID) %>% 
  summarise(sum = sum(t_mean,na.rm = T))

ta2_UKA_RS_21 <- a2_UKA_RS_21 %>% 
  group_by(BLOCK,TREATMENT,CARABID) %>% 
  summarise(sum = sum(t_mean,na.rm = T))

fa1_UKA_RS_21 <- ta1_UKA_RS_21 %>% 
  filter(CARABID %in% carabids_autumn)

fa2_UKA_RS_21 <- ta2_UKA_RS_21 %>% 
  filter(CARABID %in% carabids_autumn)

ssa1_UKA_RS_21 <- fa1_UKA_RS_21 %>% 
  group_by(TREATMENT,CARABID) %>% 
  summarise(mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

ssa2_UKA_RS_21 <- fa2_UKA_RS_21 %>% 
  group_by(TREATMENT,CARABID) %>% 
  summarise(mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

at_UKA_RS_21 <- a3_UKA_RS_21 %>% 
  group_by(DATE,TREATMENT,CARABID) %>% 
  summarise(mean = mean(t_mean, na.rm = T),
            se = sd(t_mean, na.rm = T)/sqrt(length(t_mean)),
            n = length(t_mean))

att_UKA_RS_21 <- at_UKA_RS_21 %>% 
  filter(CARABID == "TOTAL_CARABIDS")

surr_UKA_RS_21 <- subset(gc_UKA_RS_21, DATE < "2021-07-31" & 
                                CARABID != "TOTAL_CARABIDS" &
                                TREATMENT == "WT") %>% 
  group_by(DATE,BLOCK,TREATMENT,SURROUNDING) %>% 
  summarise(sum = sum(COUNT)) %>%
  group_by(TREATMENT,SURROUNDING,BLOCK) %>%
  summarise(mean = mean(sum, na.rm = T)) %>%
  group_by(TREATMENT,SURROUNDING) %>% 
  summarise(t_mean = mean(mean, na.rm = T),
            se = sd(mean, na.rm = T)/sqrt(length(mean)),
            n = length(mean))

surr_spec_UKA_RS_21 <- subset(gc_UKA_RS_21, DATE < "2021-07-31" & 
                                TREATMENT == "WT" &
                                CARABID %in% carabids) %>%
  group_by(BLOCK,TREATMENT,SURROUNDING,CARABID) %>% 
  summarise(mean1 = mean(COUNT, na.rm = T)) %>%
  group_by(TREATMENT,SURROUNDING,CARABID) %>% 
  summarise(t_mean = mean(mean1, na.rm = T),
            se = sd(mean1, na.rm = T)/sqrt(length(mean1)),
            n = length(mean1))

surr_date_UKA_RS_21 <- subset(gc_UKA_RS_21, DATE < "2021-07-31" & 
                                CARABID != "TOTAL_CARABIDS" &
                                TREATMENT == "WT") %>% 
  group_by(DATE,BLOCK,TREATMENT,SURROUNDING) %>% 
  summarise(sum = sum(COUNT)) %>%
  group_by(TREATMENT,SURROUNDING,DATE) %>%
  summarise(t_mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

# 2021b (MS)
c_UKA_MS_21<-read.csv("2021b carabids.csv",sep=";",header=T)
c_UKA_MS_21 <- c_UKA_MS_21 %>%
  mutate(across(c(ROW, DIRECTION, SURROUNDING, BLOCK, TREATMENT), 
                as.factor),
         DATE = as.Date(DATE),
         TOTAL_CARABIDS = rowSums(across(7:last_col()), na.rm = T))

gc_UKA_MS_21 <- c_UKA_MS_21 %>%
  gather("CARABID","COUNT",7:last_col()) %>%
  mutate(CARABID = as.factor(CARABID))

sc_UKA_MS_21 <- gc_UKA_MS_21 %>% 
  filter(DATE < "2021-07-31",
         !(TREATMENT %in% c("WT","S")) | SURROUNDING == "mulch") %>%
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

tc_UKA_MS_21 <- sc_UKA_MS_21 %>% 
  group_by(BLOCK,TREATMENT,CARABID) %>% 
  summarise(sum = sum(t_mean,na.rm = T))

fc_UKA_MS_21 <- tc_UKA_MS_21 %>% 
  filter(CARABID %in% carabids)

ssc_UKA_MS_21 <- fc_UKA_MS_21 %>% 
  group_by(TREATMENT,CARABID) %>% 
  summarise(mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

t_UKA_MS_21 <- sc_UKA_MS_21 %>% 
  group_by(DATE,TREATMENT,CARABID) %>% 
  summarise(mean = mean(t_mean, na.rm = T),
            se = sd(t_mean, na.rm = T)/sqrt(length(t_mean)),
            n = length(t_mean))

tt_UKA_MS_21 <- t_UKA_MS_21 %>% 
  filter(CARABID == "TOTAL_CARABIDS")

a_UKA_MS_21 <- subset(gc_UKA_MS_21, DATE > "2021-07-31") %>% 
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

ta_UKA_MS_21 <- a_UKA_MS_21 %>% 
  group_by(BLOCK,TREATMENT,CARABID) %>% 
  summarise(sum = sum(t_mean,na.rm = T))

fa_UKA_MS_21 <- ta_UKA_MS_21 %>% 
  filter(CARABID %in% carabids_autumn)

ssa_UKA_MS_21 <- fa_UKA_MS_21 %>% 
  group_by(TREATMENT,CARABID) %>% 
  summarise(mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

at_UKA_MS_21 <- a_UKA_MS_21 %>% 
  group_by(DATE,TREATMENT,CARABID) %>% 
  summarise(mean = mean(t_mean, na.rm = T),
            se = sd(t_mean, na.rm = T)/sqrt(length(t_mean)),
            n = length(t_mean))

att_UKA_MS_21 <- at_UKA_MS_21 %>% 
  filter(CARABID == "TOTAL_CARABIDS")

surr_UKA_MS_21 <- subset(gc_UKA_MS_21, DATE < "2021-07-31" & 
                           CARABID != "TOTAL_CARABIDS" &
                           TREATMENT %in% c("WT", "S")) %>%
  group_by(DATE,BLOCK,TREATMENT,SURROUNDING) %>% 
  summarise(sum = sum(COUNT)) %>%
  group_by(TREATMENT,SURROUNDING,BLOCK) %>%
  summarise(mean = mean(sum, na.rm = T)) %>%
  group_by(TREATMENT,SURROUNDING) %>% 
  summarise(t_mean = mean(mean, na.rm = T),
            se = sd(mean, na.rm = T)/sqrt(length(mean)),
            n = length(mean))

surr_spec_UKA_MS_21 <- subset(gc_UKA_MS_21, DATE < "2021-07-31" & 
                                TREATMENT %in% c("WT", "S") &
                                CARABID %in% carabids) %>%
  group_by(BLOCK,TREATMENT,SURROUNDING,CARABID) %>% 
  summarise(mean1 = mean(COUNT, na.rm = T)) %>%
  group_by(TREATMENT,SURROUNDING,CARABID) %>% 
  summarise(t_mean = mean(mean1, na.rm = T),
            se = sd(mean1, na.rm = T)/sqrt(length(mean1)),
            n = length(mean1))

surr_date_UKA_MS_21 <- subset(gc_UKA_MS_21, DATE < "2021-07-31" & 
                           CARABID != "TOTAL_CARABIDS" &
                           TREATMENT %in% c("WT", "S")) %>%
  group_by(DATE,BLOCK,TREATMENT,SURROUNDING) %>% 
  summarise(sum = sum(COUNT)) %>%
  group_by(TREATMENT,SURROUNDING,DATE) %>%
  summarise(t_mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

# 2022 (MS)
c_UKA_MS_22<-read.csv("2022 carabids.csv",sep=";",header=T)
c_UKA_MS_22 <- c_UKA_MS_22 %>%
  mutate(across(c(ROW, BLOCK, TREATMENT), 
                as.factor),
         DATE = as.Date(DATE),
         TOTAL_CARABIDS = rowSums(across(7:last_col()), na.rm = T))

gc_UKA_MS_22 <- c_UKA_MS_22 %>%
  gather("CARABID","COUNT",c(5:last_col())) %>%
  mutate(CARABID = as.factor(CARABID))

sc_UKA_MS_22 <- gc_UKA_MS_22 %>% 
  filter(DATE < "2022-07-31") %>%
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

tc_UKA_MS_22 <- sc_UKA_MS_22 %>% 
  group_by(BLOCK,TREATMENT,CARABID) %>% 
  summarise(sum = sum(t_mean,na.rm = T))

fc_UKA_MS_22 <- tc_UKA_MS_22 %>% 
  filter(CARABID %in% carabids)

ssc_UKA_MS_22 <- fc_UKA_MS_22 %>% 
  group_by(TREATMENT,CARABID) %>% 
  summarise(mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

t_UKA_MS_22 <- sc_UKA_MS_22 %>% 
  group_by(DATE,TREATMENT,CARABID) %>% 
  summarise(mean = mean(t_mean, na.rm = T),
            se = sd(t_mean, na.rm = T)/sqrt(length(t_mean)),
            n = length(t_mean))

tt_UKA_MS_22 <- t_UKA_MS_22 %>% 
  filter(CARABID == "TOTAL_CARABIDS")


# 2021a (RS) following year
c_UKA_RundS_21_fy<-read.csv("2021a carabids - following year.csv",sep=";",header=T)
c_UKA_RundS_21_fy <- c_UKA_RundS_21_fy %>%
  mutate(across(c(DIRECTION, BLOCK, TREATMENT), 
                as.factor),
         DATE = as.Date(DATE),
         TOTAL_CARABIDS = rowSums(across(5:(length(c_UKA_RundS_21_fy)-4)), na.rm = TRUE))

gc_UKA_RundS_21_fy <- c_UKA_RundS_21_fy %>%
  gather("CARABID","COUNT",c(5:24)) %>%
  mutate(CARABID = as.factor(CARABID))

sc_UKA_RundS_21_fy <- gc_UKA_RundS_21_fy %>% 
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

tc_UKA_RundS_21_fy <- sc_UKA_RundS_21_fy %>% 
  group_by(BLOCK,TREATMENT,CARABID) %>% 
  summarise(sum = sum(t_mean,na.rm = T))

fc_UKA_RundS_21_fy <- tc_UKA_RundS_21_fy %>% 
  filter(CARABID %in% carabids)

ssc_UKA_RundS_21_fy <- fc_UKA_RundS_21_fy %>% 
  group_by(TREATMENT,CARABID) %>% 
  summarise(mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

t_UKA_RundS_21_fy <- sc_UKA_RundS_21_fy %>% 
  group_by(DATE,TREATMENT,CARABID) %>% 
  summarise(mean = mean(t_mean, na.rm = T),
            se = sd(t_mean, na.rm = T)/sqrt(length(t_mean)),
            n = length(t_mean))

tt_UKA_RundS_21_fy <- t_UKA_RundS_21_fy %>% 
  filter(CARABID == "TOTAL_CARABIDS")


# 2021b (MS) following year
c_UKA_MS_21_fy<-read.csv("2021b carabids - following year.csv",sep=";",header=T)
c_UKA_MS_21_fy <- c_UKA_MS_21_fy %>%
  mutate(across(c(DIRECTION, BLOCK, TREATMENT), 
                as.factor),
         DATE = as.Date(DATE),
         TOTAL_CARABIDS = rowSums(across(5:last_col()), na.rm = TRUE))

gc_UKA_MS_21_fy <- c_UKA_MS_21_fy %>%
  gather("CARABID","COUNT",c(5:last_col())) %>%
  mutate(CARABID = as.factor(CARABID))

sc_UKA_MS_21_fy <- gc_UKA_MS_21_fy %>% 
  group_by(DATE,BLOCK,TREATMENT,CARABID) %>% 
  summarise(t_mean = mean(COUNT, na.rm = T))

tc_UKA_MS_21_fy <- sc_UKA_MS_21_fy %>% 
  group_by(BLOCK,TREATMENT,CARABID) %>% 
  summarise(sum = sum(t_mean,na.rm = T))

fc_UKA_MS_21_fy <- tc_UKA_MS_21_fy %>% 
  filter(CARABID %in% carabids)

ssc_UKA_MS_21_fy <- fc_UKA_MS_21_fy %>% 
  group_by(TREATMENT,CARABID) %>% 
  summarise(mean = mean(sum, na.rm = T),
            se = sd(sum, na.rm = T)/sqrt(length(sum)),
            n = length(sum))

t_UKA_MS_21_fy <- sc_UKA_MS_21_fy %>% 
  group_by(DATE,TREATMENT,CARABID) %>% 
  summarise(mean = mean(t_mean, na.rm = T),
            se = sd(t_mean, na.rm = T)/sqrt(length(t_mean)),
            n = length(t_mean))

tt_UKA_MS_21_fy <- t_UKA_MS_21_fy %>% 
  filter(CARABID == "TOTAL_CARABIDS")

# Calculations for Table A.2 ####
summarise_gc <- function(df) {
  df %>%
    group_by(CARABID) %>%
    summarise(TOTAL_COUNT = sum(COUNT, na.rm = TRUE)) %>%
    arrange(desc(TOTAL_COUNT))
}
gc_UKA_RS_21_sum      <- summarise_gc(gc_UKA_RS_21)
gc_UKA_MS_21_sum      <- summarise_gc(gc_UKA_MS_21)
gc_UKA_MS_22_sum      <- summarise_gc(gc_UKA_MS_22)
gc_UKA_RundS_21_fy_sum <- summarise_gc(gc_UKA_RundS_21_fy) # Relatively high H. rufipes catches in this dataset
gc_UKA_MS_21_fy_sum   <- summarise_gc(gc_UKA_MS_21_fy) # Relatively high H. rufipes catches in this dataset

all_gc <- bind_rows(
  gc_UKA_RS_21  %>% dplyr::select(CARABID, COUNT),
  gc_UKA_MS_21  %>% dplyr::select(CARABID, COUNT),
  gc_UKA_MS_22  %>% dplyr::select(CARABID, COUNT),
  gc_UKA_RundS_21_fy %>% dplyr::select(CARABID, COUNT),
  gc_UKA_MS_21_fy %>% dplyr::select(CARABID, COUNT)
)
sum_per_species <- all_gc %>%
  group_by(CARABID) %>%
  summarise(TOTAL_COUNT = sum(COUNT, na.rm = TRUE)) %>%
  arrange(desc(TOTAL_COUNT))
sum_per_species

# Calculations for results and GLMM ####
a <- gc_UKA_RS_21 %>%
  filter(DATE < "2021-07-31",
         !(TREATMENT %in% c("WT","S")) | SURROUNDING == "mulch",
         CARABID == "TOTAL_CARABIDS") %>%
  group_by(TREATMENT,BLOCK) %>%
  summarise(mean = mean(COUNT, na.rm = T))
a_sum <- a %>%
  group_by(TREATMENT) %>%
  summarise(mean = mean(mean))
a_sum
11.2/10.7
23.5/10.7
M.p <- glmmTMB(mean ~ TREATMENT + (1 | BLOCK), 
               family = "gaussian", 
               data = a)
summary(M.p)
M.p_simout <- simulateResiduals(M.p, n = 1000)
testDispersion(M.p_simout)
plot(M.p_simout, quantreg = FALSE)
summary(glht(M.p, mcp(TREATMENT = "Dunnett"))) 


b <- gc_UKA_MS_21 %>%
  filter(DATE < "2021-07-31",
         !(TREATMENT %in% c("WT","S")) | SURROUNDING == "mulch",
         CARABID == "TOTAL_CARABIDS") %>%
  group_by(TREATMENT,BLOCK) %>%
  group_by(TREATMENT, BLOCK) %>%
  summarise(mean = mean(COUNT, na.rm = T))
b_sum <- b %>%
  group_by(TREATMENT) %>%
  summarise(mean = mean(mean))
b_sum
10.1/8.19
12.3/8.19
M.p <- glmmTMB(mean ~ TREATMENT + (1 | BLOCK), 
               family = "gaussian", 
               data = b)
summary(M.p)
M.p_simout <- simulateResiduals(M.p, n = 1000)
testDispersion(M.p_simout)
plot(M.p_simout, quantreg = FALSE)
summary(glht(M.p, mcp(TREATMENT = "Dunnett"))) 

c <- gc_UKA_MS_22 %>%
  filter(CARABID == "TOTAL_CARABIDS") %>%
  group_by(TREATMENT,BLOCK) %>%
  group_by(TREATMENT, BLOCK) %>%
  summarise(mean = mean(COUNT, na.rm = T))
c_sum <- c %>%
  group_by(TREATMENT) %>%
  summarise(mean = mean(mean))
c_sum
11.1/13.0
11.5/13.0
M.p <- glmmTMB(mean ~ TREATMENT + (1 | BLOCK), 
               family = "gaussian", 
               data = c)
summary(M.p)


gc_UKA_RS_21 %>%
  filter(DATE < "2021-07-31",
         !(TREATMENT %in% c("WT","S")) | SURROUNDING == "mulch") %>%
  filter(CARABID == "POECUP") %>%
  group_by(TREATMENT) %>%
  summarise(mean = mean(COUNT, na.rm = T))
4.72/3.5
11.0/3.5

gc_UKA_MS_21 %>%
  filter(DATE < "2021-07-31",
         !(TREATMENT %in% c("WT","S")) | SURROUNDING == "mulch") %>%
  filter(CARABID == "POECUP") %>%
  group_by(TREATMENT) %>%
  summarise(mean = mean(COUNT, na.rm = T))
1.24/0.265
1.15/0.265

gc_UKA_MS_22 %>%
  filter(CARABID == "POECUP") %>%
  group_by(TREATMENT) %>%
  summarise(mean = mean(COUNT, na.rm = T))
0.946/0.232
2.09/0.232


gc_UKA_RS_21 %>%
  filter(DATE < "2021-07-31",
         !(TREATMENT %in% c("WT","S")) | SURROUNDING == "mulch") %>%
  filter(CARABID == "BEMLAM" | CARABID == "BEMQUA") %>%
  group_by(TREATMENT, CARABID) %>%
  summarise(mean = mean(COUNT, na.rm = T))
1.42/1.15 #BEMLAM, S/C
4.50/1.15 #BEMLAM, WT/C
1.65/1.24 #BEMQUA, S/C
2.09/1.24 #BEMQUA, WT/C

gc_UKA_MS_21 %>%
  filter(DATE < "2021-07-31",
         !(TREATMENT %in% c("WT","S")) | SURROUNDING == "mulch") %>%
  filter(CARABID == "BEMLAM" | CARABID == "BEMQUA") %>%
  group_by(TREATMENT, CARABID) %>%
  summarise(mean = mean(COUNT, na.rm = T))
3.12/1.71 #BEMLAM, S/C
4.06/1.71 #BEMLAM, WT/C
1.85/1.65 #BEMQUA, S/C
3.35/1.65 #BEMQUA, WT/C

gc_UKA_MS_22 %>%
  filter(CARABID == "BEMLAM" | CARABID == "BEMQUA") %>%
  group_by(TREATMENT, CARABID) %>%
  summarise(mean = mean(COUNT, na.rm = T))
7.00/6.95 #BEMLAM, SIL/C
5.84/6.95 #BEMLAM, WT/C
0.357/1.09 #BEMQUA, SIL/C
0.768/1.09 #BEMQUA, WT/C


# Catch by date plots ####

# 2021a plot
sampling_A <- tribble(
  ~sampling, ~start_date, ~end_date,
  "sampling", "2021-06-04", "2021-06-07",
  "sampling", "2021-06-07", "2021-06-10",
  "sampling", "2021-06-10", "2021-06-14",
  "sampling", "2021-06-14", "2021-06-17",
  "no sampling", "2021-06-17", "2021-07-07",
  "sampling", "2021-07-07", "2021-07-12",
  "sampling", "2021-07-12", "2021-07-16")
sampling_A <- sampling_A %>%
  mutate(start_date = as.POSIXct(start_date),
         end_date = as.POSIXct(end_date))
tt_UKA_RS_21 <- tt_UKA_RS_21 %>%
  mutate(DATE = as.POSIXct(DATE))
tt_UKA_RS_21 <- tt_UKA_RS_21 %>%
  mutate(row_id = row_number())

ggplot() +
  geom_col(data = tt_UKA_RS_21, aes(x=DATE-hours(2), y=mean, fill=TREATMENT),
           position="dodge",width = 140000) +
  geom_segment(data = sampling_A, aes(color = sampling, x=start_date+hours(1), 
                                      xend=end_date-hours(1), y=-1, yend=-1), 
               linewidth = 0.7) +
  geom_errorbar(data = tt_UKA_RS_21, aes(ymin = mean - se, ymax = mean + se,
                                          x=DATE-hours(2), y=mean, fill=TREATMENT), 
                width = 45000, position = position_dodge(140000)) +
  labs(title = "carabids - trial 2021a", 
       x = "", y = "mean catch per trap") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control","straw","triticale/vetch")) +
  scale_color_manual(name = NULL,
                     values = c("cadetblue3","darkgreen")) +
  theme_light() +
  theme(
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = 0.5),
    plot.title = element_text(hjust = 0.98, size = 11, vjust = 0),
    legend.position = "none",
    panel.grid.minor.x = element_blank(),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -5, r=10)) +
  scale_x_datetime(breaks = breaks_2021, 
                   date_labels = "%d %b",
                   limits = as.POSIXct(c("2021-06-04","2021-07-18"))) +
  scale_y_continuous(
    breaks = seq(0, 100, by = 10),
    limits = c(-2, NA),
    expand = expansion(mult = c(0, 0.02)))

ggsave(filename = "Carabids Dates - trial 2021a.pdf",
       width = 17, height = 5.5, units = "cm")


# 2021b plot
sampling_B <- tribble(
  ~sampling, ~start_date, ~end_date,
  "no sampling", "2021-06-04", "2021-06-11",
  "sampling", "2021-06-11", "2021-06-15",
  "sampling", "2021-06-15", "2021-06-18",
  "sampling", "2021-06-18", "2021-06-22",
  "sampling", "2021-06-22", "2021-06-25",
  "no sampling", "2021-06-25", "2021-07-08",
  "sampling", "2021-07-08", "2021-07-12",
  "sampling", "2021-07-12", "2021-07-16")
sampling_B <- sampling_B %>%
  mutate(start_date = as.POSIXct(start_date),
         end_date = as.POSIXct(end_date))
tt_UKA_MS_21 <- tt_UKA_MS_21 %>%
  mutate(DATE = as.POSIXct(DATE))
tt_UKA_MS_21 <- tt_UKA_MS_21 %>%
  mutate(row_id = row_number())

ggplot() +
  geom_col(data = tt_UKA_MS_21, aes(x=DATE-hours(2), y=mean, fill=TREATMENT),
           position="dodge",width = 140000) +
  geom_segment(data = sampling_B, aes(color = sampling, x=start_date+hours(1), 
                                      xend=end_date-hours(1), y=-1, yend=-1), 
               linewidth = 0.7) +
  geom_errorbar(data = tt_UKA_MS_21, aes(ymin = mean - se, ymax = mean + se,
                                          x=DATE-hours(2), y=mean, fill=TREATMENT), 
                width = 45000, position = position_dodge(140000)) +
  labs(title = "carabids - trial 2021b", 
       x = "", y = "mean catch per trap") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control","straw","triticale/vetch")) +
  scale_color_manual(name = NULL,
                     values = c("cadetblue3","darkgreen")) +
  theme_light() +
  theme(
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = 0.5),
    plot.title = element_text(hjust = 0.98, size = 11, vjust = 0),
    legend.position = "none",
    panel.grid.minor.x = element_blank(),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -5, r=10)) +
  scale_x_datetime(breaks = breaks_2021, 
                   date_labels = "%d %b",
                   limits = as.POSIXct(c("2021-06-04","2021-07-18"))) +
  scale_y_continuous(
    breaks = seq(0, 100, by = 10),
    limits = c(-2, NA),
    expand = expansion(mult = c(0, 0.02)))

ggsave(filename = "Carabids Dates - trial 2021b.pdf",
       width = 17, height = 4, units = "cm")


# 2022 plot
sampling_C <- tribble(
  ~sampling, ~start_date, ~end_date,
  "no sampling", "2022-06-04", "2022-06-10",
  "sampling", "2022-06-10", "2022-06-13",
  "sampling", "2022-06-13", "2022-06-15",
  "no sampling", "2022-06-15", "2022-06-27",
  "sampling", "2022-06-27", "2022-06-30",
  "sampling", "2022-06-30", "2022-07-04",
  "sampling", "2022-07-04", "2022-07-07",
  "sampling", "2022-07-07", "2022-07-11",
  "sampling", "2022-07-11", "2022-07-14",
  "no sampling", "2022-07-14", "2022-07-16")
sampling_C <- sampling_C %>%
  mutate(start_date = as.POSIXct(start_date),
         end_date = as.POSIXct(end_date))
tt_UKA_MS_22 <- tt_UKA_MS_22 %>%
  mutate(DATE = as.POSIXct(DATE))

ggplot() +
  geom_col(data = tt_UKA_MS_22, aes(x=DATE-hours(2), y=mean, fill=TREATMENT),
           position="dodge",width = 140000) +
  geom_segment(data = sampling_C, aes(color = sampling, x=start_date+hours(1), 
                                      xend=end_date-hours(1), y=-1, yend=-1), 
               linewidth = 0.7) +
  geom_errorbar(data = tt_UKA_MS_22, aes(ymin = mean - se, ymax = mean + se,
                                                x=DATE-hours(2), y=mean, fill=TREATMENT), 
                width = 45000, position = position_dodge(140000)) +
  labs(title = "carabids - trial 2022", 
       x = "", y = "mean catch per trap") +
  scale_fill_manual(name = NULL,
                    values =c("#EE6677","#228833","#CCBB44","#4477AA"),
                    labels = c("control","silage","triticale/vetch","straw"), 
                    drop = FALSE) +
  scale_color_manual(name = NULL,
                     values = c("cadetblue3","darkgreen")) +
  theme_light() +
  theme(
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    legend.text = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = 0.5),
    plot.title = element_text(hjust = 0.98, size = 11, vjust = 0),
    panel.grid.minor.x = element_blank(),
    legend.position = "none",
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -5, r=10)) +
  scale_x_datetime(breaks = breaks_2022, 
                   date_labels = "%d %b",
                   limits = as.POSIXct(c("2022-06-04","2022-07-18"))) +
  scale_y_continuous(
    breaks = seq(0, 100, by = 10),
    limits = c(-2, NA),
    expand = expansion(mult = c(0, 0.02)))

ggsave(filename = "Carabids Dates - trial 2022.pdf",
       width = 17, height = 3.5, units = "cm")


# Total catch by species plots ####

# 2021a plot
trt_order_car3 <- c('C', 'S', 'WT')
ggplot(ssc_UKA_RS_21, aes(x = CARABID, y = mean, fill = factor(TREATMENT, levels = trt_order_car3))) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1, 
                position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap", 
       title = "carabids - trial 2021a - June/July") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control", "straw","triticale/vetch")) +
  theme_light() +
  scale_x_discrete(labels = carabid_labels) +
  theme(
    axis.text.x = element_text(size = 8.5, angle = 45, hjust = 1, lineheight = 0.75),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.95, size = 11, vjust = 0),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -11, r = 1),
    legend.position = "none") +
  scale_y_continuous(
    limits = c(-0.8, 76),
    breaks = seq(0, 70, by = 10),
    expand = c(0, 0))
ggsave(filename = "Carabid Species - trial 2021a.pdf",
       width = 7.6, height = 7.4, units = "cm")

# 2021b plot
ggplot(ssc_UKA_MS_21, aes(x = CARABID, y = mean, fill = factor(TREATMENT, levels = trt_order_car3))) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1, 
                position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap", 
       title = "carabids - trial 2021b - June/July") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control", "straw","triticale/vetch")) +
  theme_light() +
  scale_x_discrete(labels = carabid_labels) +
  theme(
    axis.text.x = element_text(size = 8.5, angle = 45, hjust = 1, lineheight = 0.75),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.98, size = 11, vjust = 0),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -11, r = 1),
    legend.position = "none") +
  scale_y_continuous(
    limits = c(-0.8, 76),
    breaks = seq(0, 70, by = 10),
    expand = c(0, 0))
ggsave(filename = "Carabid Species - trial 2021b.pdf",
       width = 7.6, height = 7.4, units = "cm")

# 2022 plot
trt_order_car4 <- c('C', 'SIL', 'WT')
ggplot(ssc_UKA_MS_22, aes(x = CARABID, y = mean, fill = factor(TREATMENT, levels = trt_order_car4))) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1, 
                position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap", 
       title = "carabids - trial 2022 - June/July") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#228833","#CCBB44"),
                    labels = c("control", "silage", "triticale/vetch")) +
  theme_light() +
  scale_x_discrete(labels = carabid_labels) +
  theme(
    axis.text.x = element_text(size = 8.5, angle = 45, hjust = 1, lineheight = 0.75),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.98, size = 11, vjust = 0),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -11, r = 1),
    legend.position = "none") +
  scale_y_continuous(
    limits = c(-0.8, 76),
    breaks = seq(0, 70, by = 10),
    expand = c(0, 0))
ggsave(filename = "Carabid Species - trial 2022.pdf",
       width = 7.6, height = 7.4, units = "cm")

# 2021a September plot
ggplot(ssa1_UKA_RS_21, aes(x = CARABID, y = mean, fill = factor(TREATMENT, levels = trt_order_car3))) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1, 
                position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap", 
       title = "carabids - trial 2021a - early September") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control", "straw","triticale/vetch")) +
  theme_light() +
  scale_x_discrete(labels = carabid_labels) +
  theme(
    axis.text.x = element_text(size = 8.5, angle = 45, hjust = 1, lineheight = 0.75),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.98, size = 10, vjust = 0),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -11, r = 1),
    legend.position = "none") +
  scale_y_continuous(
    limits = c(-0.4, 38),
    breaks = seq(0, 35, by = 5),
    expand = c(0, 0))
ggsave(filename = "Carabid Species - trial 2021a - September.pdf",
       width = 7.6, height = 6.5, units = "cm")


# 2021a October
ggplot(ssa2_UKA_RS_21, aes(x = CARABID, y = mean, fill = factor(TREATMENT, levels = trt_order_car3))) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1, 
                position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap", 
       title = "carabids - trial 2021a - October (after harvest)") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control", "straw","triticale/vetch")) +
  theme_light() +
  scale_x_discrete(labels = carabid_labels) +
  theme(
    axis.text.x = element_text(size = 8.5, angle = 45, hjust = 1, lineheight = 0.75),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.98, size = 10, vjust = 0),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -11, r = 1),
    legend.position = "none") +
  scale_y_continuous(
    limits = c(-0.4, 38),
    breaks = seq(0, 35, by = 5),
    expand = c(0, 0))
ggsave(filename = "Carabid Species - trial 2021a - October.pdf",
       width = 7.6, height = 6.5, units = "cm")


# 2021b September plot
ggplot(ssa_UKA_MS_21, aes(x = CARABID, y = mean, fill = factor(TREATMENT, levels = trt_order_car3))) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1, 
                position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap", 
       title = "carabids - trial 2021b - late August") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control", "straw", "triticale/vetch")) +
  theme_light() +
  scale_x_discrete(labels = carabid_labels) +
  theme(
    axis.text.x = element_text(size = 8.5, angle = 45, hjust = 1, lineheight = 0.75),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.98, size = 10, vjust = 0),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -11, r = 1),
    legend.position = "none") +
  scale_y_continuous(
    limits = c(-0.4, 38),
    breaks = seq(0, 35, by = 5),
    expand = c(0, 0))
ggsave(filename = "Carabid Species - trial 2021b - August.pdf",
       width = 7.6, height = 6.5, units = "cm")


# 2021a following year plot
ggplot(ssc_UKA_RundS_21_fy, aes(x = CARABID, y = mean, fill = factor(TREATMENT, levels = trt_order_car3))) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1, 
                position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap", 
       title = "carabids - trial 2021a - following year") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control", "straw", "triticale/vetch")) +
  theme_light() +
  scale_x_discrete(labels = carabid_labels) +
  theme(
    axis.text.x = element_text(size = 8.5, angle = 45, hjust = 1, lineheight = 0.75),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.98, size = 11, vjust = 0),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -11, r = 1),
    legend.position = "none") +
  scale_y_continuous(
    limits = c(-0.6, 62),
    breaks = seq(0, 60, by = 10),
    expand = c(0, 0))
ggsave(filename = "Carabid Species - trial 2021a - following year.pdf",
       width = 7.6, height = 7.6, units = "cm")


# 2021b following year plot
ggplot(ssc_UKA_MS_21_fy, aes(x = CARABID, y = mean, fill = factor(TREATMENT, levels = trt_order_car3))) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1, 
                position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap", 
       title = "carabids - trial 2021b - following year") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control", "straw", "triticale/vetch")) +
  theme_light() +
  scale_x_discrete(labels = carabid_labels) +
  theme(
    axis.text.x = element_text(size = 8.5, angle = 45, hjust = 1, lineheight = 0.75),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.98, size = 11, vjust = 0),
    panel.border = element_rect(colour = "lightgray", size = 0.5),
    plot.margin = margin(b = -11, r = 1),
    legend.position = "none") +
  scale_y_continuous(
    limits = c(-0.6, 62),
    breaks = seq(0, 60, by = 10),
    expand = c(0, 0))
ggsave(filename = "Carabid Species - trial 2021b - following year.pdf",
       width = 7.6, height = 7.6, units = "cm")


# Legend 1
dummy_plot <- ggplot(data.frame(x = 1:4, y = 1:4, fill = c("control", "triticale/vetch", "silage", "straw")), 
                     aes(x, y, fill = fill)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44","#228833"),
                    labels = c("control","straw","triticale/vetch","silage")) +  
  theme_void() +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10))
leg <- get_legend(dummy_plot)
as_ggplot(leg)
ggsave(filename = "legend.pdf",
       width = 4.25, height = 0.3)

# Legend 2
dummy_plot <- ggplot(data.frame(x = 1:3, y = 1:3, fill = c("control", "triticale/vetch", "straw")), 
                     aes(x, y, fill = fill)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(name = "treatments:",
                    values =c("#EE6677","#4477AA","#CCBB44"),
                    labels = c("control", "straw", "triticale/vetch")) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10))
dummy_plot
leg <- get_legend(dummy_plot)
as_ggplot(leg)
ggsave(filename = "legend_wo_SIL.pdf",
       width = 3.5, height = 0.3)


# Surroundings ####

combined_surr <- bind_rows(surr_spec_UKA_RS_21, surr_spec_UKA_MS_21, 
                           .id = "TRIAL") %>%
  mutate(TRIAL = case_when(
    TRIAL == "1" ~ "A",
    TRIAL == "2" ~ "B")) %>%
  mutate(TRIAL_TREATMENT = paste(TRIAL,TREATMENT,sep = "_"))
names2 <- c(
  `A_WT` = "carabids\ntrial 2021a - triticale/vetch",
  `B_S` = "carabids\ntrial 2021b - straw",
  `B_WT` = "carabids\ntrial 2021b - triticale/vetch"
)

ggplot(combined_surr, aes(x = CARABID, y = t_mean, fill = SURROUNDING)) +
  geom_bar(position = "dodge", stat = "identity", width = 0.6) +
  geom_errorbar(aes(ymin = t_mean - se, ymax = t_mean + se), 
                width = 0.1, position = position_dodge(0.6)) +
  labs(x = "", y = "mean catch sums per trap") +
  scale_fill_manual(name = "surrounding:",
                    values =c("#1B9E77","#E7298A"),
                    labels = c("mulch","no mulch"))+
  scale_x_discrete(labels = carabid_labels) +
  theme_light() +
  theme(
    axis.text.x = element_text(size = 7.5, angle = 45, hjust = 1, lineheight = 0.7),
    axis.text.y = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90),
    plot.title = element_text(hjust = 0.98, size = 9, vjust = 0),
    plot.margin = margin(b = -11, r = 1),
    strip.text = element_text(
      size = 10,
      color = "black"),
    strip.background = element_rect(
      fill = "gray90",
      color = "gray70"),
    legend.position = "top") +
  facet_grid(.~TRIAL_TREATMENT, labeller = as_labeller(names2))

ggsave(filename = "Surrounding - carabid species.pdf",
       width = 17, height = 10, units = "cm")