#  Section 1: Data description ----

# This data set contains data from a 3 year study on parasitism rates in diversified cropping systems.

# Measured variables are: 
# - Presence of caterpillars (10 caterpillars were placed, but X returned)
# - Parasitism of caterpillar (1 = parasitized by any parasitoid, 0 = not parasitized)
# - Number of eggs of C. glomerata (Cg), C. rubecula (Cr) and H. ebininus (He). 
#   * In 2019 and 2020 H. ebininus was not recognized as such, so Cr in that year contains both species.
#   * In 2021 H. ebininus was seperated from C. rubecula.

# Environmental variables:
# - Treatments (6 levels): 
#   * REF_SPACE
#   * STRIP
#   * STRIP_VAR
#   * STRIP_ADD
#   * STRIP_ROT
#   * PIXEL
# - Cabbage cultivar (3 levels):
#   * Rivera
#   * Christmas Drumhead (only measured in STRIP_VAR)
#   * Expect (was used instead of Rivera in Lelystad in 2021)
# - Location (2 levels):
#   * Droevendaal (Wageningen)
#   * Broekemahoeve (Lelystad)
# - Fields (9 levels):
#   * Droevendaal:
#     + 1
#     + 2
#     + 3
#     + 4 (REF_SPACE and STRIP in 2021)
#     + 10 (REF_SPACE and STRIP in 2019 and 2020)
#   * Broekemahoeve:
#     + J8_4 (including REF_SPACE in 2019)
#     + J10_1
#     + J9_1 (REF_SPACE and STRIP in 2020)
#     + J8_1 (REF_SPACE and STRIP in 2021)
# - Years (3 levels):
#   * 2019
#   * 2020
#   * 2021
# - Dates (rounds, dates on which caterpillars were recaptured, many levels)


#  Section 2: Import the data and load the packages----

#Create clean R workspace
rm(list=ls())
graphics.off()
gc(FALSE)
cat("\014") 
Sys.setenv(LANG = "en")

#' We need the following packages:

Pkgs2Load <- c("readxl", "car", "MuMIn", "DHARMa", "emmeans", "multcomp",
               "lattice", "ggplot2", "mgcv", "plyr", "MASS", 
               "pscl", "glmmTMB", "cowplot", "rgl", "pals", "ggtext")

invisible(lapply(Pkgs2Load, library, character.only = TRUE))

#' Set the working directory.
setwd("[FILL IN]")

#' Import the data.
PR <- read_xlsx("Par_data.xlsx")


# Section 3: Prepare the data----

#' Subsection 3.1: Create variables
PR$year <- format(PR$date, format = "%Y")
PR$wheat <- ifelse(PR$location == "Droevendaal", 
                   ifelse(PR$round < 4 | 
                            PR$round > 8 & PR$round < 13 |
                            PR$round > 15 & PR$round < 21, "before", "after"),
                   ifelse(PR$round < 5 | 
                            PR$round > 8 & PR$round < 14 |
                            PR$round > 15 & PR$round < 20, "before", "after"))
PR$wheat <- factor(PR$wheat, levels = c("before", "after"))

#' Subsection 3.2: Change characters to factors
cols <- c("location","treatment","cultivar","field","ID","round")
PR[cols] <- lapply(PR[cols],factor)

#' Check the imported data.
names(PR)
str(PR)

#' Subsection 3.3: Change levels treatments
PR$treatment <- factor(PR$treatment, levels = c("REF_SPACE", "STRIP", "STRIP_VAR", "STRIP_ADD", "STRIP_ROT", "PIXEL"))
PR$location <- factor(PR$location, levels = c("Droevendaal", "Broekemahoeve"))
PR$cultivar[PR$cultivar == "Expect"] <- "Rivera"

#' Subsection 3.4: Change reference field name and remove ref_space 2019
#PR <- subset(PR[PR$field != "J8_3",])

PR$field <- factor(PR$field,
                   levels = c("1","2","3","4","10","J8_1","J8_3","J8_4","J9_3","J10_1"),
                   labels = c("1", "2", "3", "4", "4", "J9","J9" , "J8", "J9", "J10"))
PR[cols] <- lapply(PR[cols],factor)

PR <- na.omit(PR)

#' Subsection 3.5: Remove non-recaptured caterpillars and unparasitized caterpillars
CgE <- subset(PR[PR$cg > 0,])

#' Subsection 3.6: Make dummy variable system
#CgE$strip <- ifelse(CgE$treatment == "REF_SPACE" | CgE$treatment == "PIXEL", 0, 1)
#CgE$pixel <- ifelse(CgE$treatment == "PIXEL", 1, 0)
#CgE$nectar <- ifelse(CgE$treatment == "STRIP_ADD" | CgE$treatment == "STRIP_ROT" | CgE$treatment == "PIXEL", 1, 0)
#CgE$intravar <- ifelse(CgE$treatment == "STRIP_VAR" | CgE$treatment == "STRIP_ROT" | CgE$treatment == "PIXEL", 1, 0)
#CgE$diverse <- ifelse(CgE$treatment == "STRIP_ROT" | CgE$treatment == "PIXEL", 1, 0)

#CgP$crops <- factor(CgP$treatment,
#                    levels = c("REF_SPACE", "STRIP", "STRIP_VAR", "STRIP_ADD", "STRIP_ROT", "PIXEL"),
#                    labels = c(1, 2, 2, 3, 6, 6))
#CgP$crops <- as.numeric(CgP$crops)
#CgP$crops <- ifelse(CgP$crops == 4, 6, CgP$crops)

#' Subsection 3.7: Split CD vs Rivera from treatment effects. 
CgET <- CgE[CgE$cultivar == "Rivera",]
CgET[cols] <- lapply(CgET[cols],factor)

CgEC <- CgE[CgE$treatment == "STRIP_VAR",]
CgEC <- CgEC[CgEC$year > 2019 | CgEC$year == 2019 & CgEC$round == "7"| CgEC$year == 2019 & CgEC$round == "8",]  #remove rounds without samples on Christmas Drumhead
CgEC$cultivar <- factor(CgEC$cultivar, levels = c("Rivera", "Christmas drumhead"))
CgEC[cols] <- lapply(CgEC[cols],factor)

#' Subsection 3.8: Split Droevendaal and Broekemahoeve
CgET.droef <- CgET[CgET$location == "Droevendaal",]
CgET.droef[cols] <- lapply(CgET.droef[cols], factor)

CgET.broek <- CgET[CgET$location == "Broekemahoeve",]
CgET.broek[cols] <- lapply(CgET.broek[cols], factor)

CgET.droef.nopix <- CgET.droef[CgET.droef$treatment != "PIXEL",]
CgET.droef.nopix[cols] <- lapply(CgET.droef.nopix[cols],factor)

#' Subsection 3.9: Remove strip_rot and pixel from both locations dataset
CgET <- CgET[CgET$treatment != "STRIP_ROT",] 
CgET <- CgET[CgET$treatment != "PIXEL",]
CgET[cols] <- lapply(CgET[cols],factor)



# Section 4: Apply the models in GLM----
# Subsection 4.1: both locations (excl strip_rot and pixel)
TM1 <- glmmTMB(cg ~ 1 + treatment +
                 location + year + wheat + 
                 (1|round) + (1|field) + (1|ID) +
                 location*year +
                 year*wheat,
               family = nbinom2,
               ziformula =~ 0,
               data = CgET, na.action = na.fail)


# Subsection 4.2: all treatments / Droevendaal
TM2 <- glmmTMB(cg ~ 1 + treatment +
                 year + wheat +
                 (1|round) + (1|field) + (1|ID) +
                 treatment*year,
               family = nbinom2,
               data = CgET.droef, na.action = na.fail)

# Subsection 4.3: all treatments / Broekemahoeve
TM3 <- glmmTMB(cg ~ 1 + treatment +
                 year + wheat +
                 (1|round) + (1|field) + (1|ID) +
                 treatment*wheat +
                 year*wheat,
               family = nbinom2,
               data = CgET.broek, na.action = na.fail)

# Subsection 4.4: cultivar effects
CM1 <- glmmTMB(cg ~ 1 + cultivar +
                 location + wheat + year + 
                 (1|round) + (1|field) + (1|ID) +
                 cultivar*wheat,
               family = nbinom2,
               data = CgEC, na.action = na.fail)


# Section 5: Results in words----
# Subsection 5.1: both locations (excl strip_rot and pixel)
Anova(TM1)
summary(TM1)
HSD.TM1.1 <- emmeans(ref_grid(TM1, transform = "response"), pairwise ~ treatment|wheat)
CLD.TM1.1 <- cld(HSD.TM1.1$emmeans, Letters = letters)
CLD.TM1.1$.group <- gsub(" ", "", CLD.TM1.1$.group)

HSD.TM1.2 <- emmeans(ref_grid(TM1, transform = "response"), pairwise ~ treatment|year)
CLD.TM1.2 <- cld(HSD.TM1.2$emmeans, Letters = letters)
CLD.TM1.2$.group <- gsub(" ", "", CLD.TM1.2$.group)

# Subsection 5.2: all treatments / Droevendaal
Anova(TM2)
summary(TM2)
HSD.TM2.1 <- emmeans(ref_grid(TM2, transform = "response"), pairwise ~ treatment|wheat)
CLD.TM2.1 <- cld(HSD.TM2.1$emmeans, Letters = letters)
CLD.TM2.1$.group <- gsub(" ", "", CLD.TM2.1$.group)

HSD.TM2.2 <- emmeans(ref_grid(TM2, transform = "response"), pairwise ~ treatment|year)
CLD.TM2.2 <- cld(HSD.TM2.2$emmeans, Letters = letters)
CLD.TM2.2$.group <- gsub(" ", "", CLD.TM2.2$.group)

# Subsection 5.3: all treatments / Broekemahoeve
Anova(TM3)
summary(TM3)
HSD.TM3.1 <- emmeans(ref_grid(TM3, transform = "response"), pairwise ~ treatment|wheat)
CLD.TM3.1 <- cld(HSD.TM3.1$emmeans, Letters = letters)
CLD.TM3.1$.group <- gsub(" ", "", CLD.TM3.1$.group)

HSD.TM3.2 <- emmeans(ref_grid(TM3, transform = "response"), pairwise ~ treatment|year)
CLD.TM3.2 <- cld(HSD.TM3.2$emmeans, Letters = letters)
CLD.TM3.2$.group <- gsub(" ", "", CLD.TM3.2$.group)

# Subsection 5.4: Cultivar effects
Anova(CM1)
summary(CM1)
HSD.CM1.1 <- emmeans(ref_grid(CM1, transform = "response"), pairwise ~ cultivar|wheat)
CLD.CM1.1 <- cld(HSD.CM1.1$emmeans, Letters = letters)
CLD.CM1.1$.group <- gsub(" ", "", CLD.CM1.1$.group)

HSD.CM1.2 <- emmeans(ref_grid(CM1, transform = "response"), pairwise ~ cultivar|year)
CLD.CM1.2 <- cld(HSD.CM1.2$emmeans, Letters = letters)
CLD.CM1.2$.group <- gsub(" ", "", CLD.CM1.2$.group)


# Section 6: Figures----
# Subsection 6.0: Preparations
Droef.colors <- c("#895F39", "#B6D663", "#27a6b0", "#213F9A", "#9671B2", "#B32E5F")
Broek.colors <- c("#895F39", "#B6D663", "#27a6b0", "#213F9A")
Cultivar.colors <- c("#FF9997", "#619CFF")

# Subsection 6.1: both locations (excl strip_rot and pixel)
png(filename = "CgE - Both locations - Wheat.png", width = 2060, height = 2000, res = 400)
ggplot(data = CLD.TM1.1, aes(x = wheat, y = response)) +
  geom_pointrange(aes(x = wheat, y = response, ymin = lower.CL, ymax = upper.CL, 
                      size = 1,
                      color = treatment, group = treatment), shape = 16, position = position_dodge(width = 0.75)) +
  geom_text(aes(label = .group, group = treatment, y = upper.CL + 2, size = 2) , position = position_dodge(width = 0.75)) +
  scale_color_manual(values = Broek.colors) +
  theme_minimal() +
  theme(title = element_text(size = 14), axis.title = element_text(size = 17), 
        axis.text.x = element_text(size = 16, color = "black"), axis.title.x = element_text(margin = margin(t = 15)),  
        axis.text.y = element_text(size = 14, color = "black"), axis.title.y = element_markdown(margin = margin(r = 15)), panel.grid.major.x = element_blank(),
        legend.position = "none") +
  labs(title = "Both locations", x = "Wheat harvest", y = "Number of *Cotesia glomerata* eggs per parasitized caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  ylim(c(0,48))
dev.off()

png(filename = "CgE - Both locations - Year.png", width = 3100, height = 2000, res = 400)
ggplot(data = CLD.TM1.2, aes(x = year, y = response)) +
  geom_pointrange(aes(x = year, y = response, ymin = lower.CL, ymax = upper.CL, 
                      size = 1,
                      color = treatment, group = treatment), shape = 16, position = position_dodge(width = 0.75)) +
  geom_text(aes(label = .group, group = treatment, y = upper.CL + 2, size = 2) , position = position_dodge(width = 0.75)) +
  scale_color_manual(values = Broek.colors) +
  theme_minimal() +
  theme(title = element_text(size = 14), axis.title = element_text(size = 17), 
        axis.text.x = element_text(size = 16, color = "black"), axis.title.x = element_text(margin = margin(t = 15)),  
        axis.text.y = element_text(size = 14, color = "black"), axis.title.y = element_markdown(margin = margin(r = 15)), panel.grid.major.x = element_blank(),
        legend.position = "none") +
  labs(title = "Both locations", x = "Year", y = "Number of *Cotesia glomerata* eggs per parasitized caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  geom_vline(xintercept = 2.5, linetype = 2) + 
  ylim(c(0,48))
dev.off()

# Subsection 6.2: All treatments / Droevendaal
png(filename = "CgE - Droevendaal - Wheat.png", width = 2060, height = 2000, res = 400)
ggplot(data = CLD.TM2.1, aes(x = wheat, y = response)) +
  geom_pointrange(aes(x = wheat, y = response, ymin = lower.CL, ymax = upper.CL, 
                      size = 1,
                      color = treatment, group = treatment), shape = 16, position = position_dodge(width = 0.75)) +
  geom_text(aes(label = .group, group = treatment, y = upper.CL + 2, size = 2) , position = position_dodge(width = 0.75)) +
  scale_color_manual(values = Droef.colors) +
  theme_minimal() +
  theme(title = element_text(size = 14), axis.title = element_text(size = 17), 
        axis.text.x = element_text(size = 16, color = "black"), axis.title.x = element_text(margin = margin(t = 15)),  
        axis.text.y = element_text(size = 14, color = "black"), axis.title.y = element_markdown(margin = margin(r = 15)), panel.grid.major.x = element_blank(),
        legend.position = "none") +
  labs(title = "Droevendaal", x = "Wheat harvest", y = "Number of *Cotesia glomerata* eggs per parasitized caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  ylim(c(0,48))
dev.off()

png(filename = "CgE - Droevendaal - Year.png", width = 3100, height = 2000, res = 400)
ggplot(data = CLD.TM2.2, aes(x = year, y = response)) +
  geom_pointrange(aes(x = year, y = response, ymin = lower.CL, ymax = upper.CL, 
                      size = 1,
                      color = treatment, group = treatment), shape = 16, position = position_dodge(width = 0.75)) +
  geom_text(aes(label = .group, group = treatment, y = upper.CL + 2, size = 2) , position = position_dodge(width = 0.75)) +
  scale_color_manual(values = Droef.colors) +
  theme_minimal() +
  theme(title = element_text(size = 14), axis.title = element_text(size = 17), 
        axis.text.x = element_text(size = 16, color = "black"), axis.title.x = element_text(margin = margin(t = 15)),  
        axis.text.y = element_text(size = 14, color = "black"), axis.title.y = element_markdown(margin = margin(r = 15)), panel.grid.major.x = element_blank(),
        legend.position = "none") +
  labs(title = "Droevendaal", x = "Year", y = "Number of *Cotesia glomerata* eggs per parasitized caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  geom_vline(xintercept = 2.5, linetype = 2) + 
  ylim(c(0,48))
dev.off()


# Subsection 6.3: All treatments / Broekemahoeve
png(filename = "CgE - Broekemahoeve - Wheat.png", width = 2060, height = 2000, res = 400)
ggplot(data = CLD.TM3.1, aes(x = wheat, y = response)) +
  geom_pointrange(aes(x = wheat, y = response, ymin = lower.CL, ymax = upper.CL, 
                      size = 1,
                      color = treatment, group = treatment), shape = 16, position = position_dodge(width = 0.75)) +
  geom_text(aes(label = .group, group = treatment, y = upper.CL + 2, size = 2) , position = position_dodge(width = 0.75)) +
  scale_color_manual(values = Broek.colors) +
  theme_minimal() +
  theme(title = element_text(size = 14), axis.title = element_text(size = 17), 
        axis.text.x = element_text(size = 16, color = "black"), axis.title.x = element_text(margin = margin(t = 15)),  
        axis.text.y = element_text(size = 14, color = "black"), axis.title.y = element_markdown(margin = margin(r = 15)), panel.grid.major.x = element_blank(),
        legend.position = "none") +
  labs(title = "Broekemahoeve", x = "Wheat harvest", y = "Number of *Cotesia glomerata* eggs per parasitized caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  ylim(c(0,48))
dev.off()

png(filename = "CgE - Broekemahoeve - Year.png", width = 3100, height = 2000, res = 400)
ggplot(data = CLD.TM3.2, aes(x = year, y = response)) +
  geom_pointrange(aes(x = year, y = response, ymin = lower.CL, ymax = upper.CL, 
                      size = 1,
                      color = treatment, group = treatment), shape = 16, position = position_dodge(width = 0.75)) +
  geom_text(aes(label = .group, group = treatment, y = upper.CL + 2, size = 2) , position = position_dodge(width = 0.75)) +
  scale_color_manual(values = Broek.colors) +
  theme_minimal() +
  theme(title = element_text(size = 14), axis.title = element_text(size = 17), 
        axis.text.x = element_text(size = 16, color = "black"), axis.title.x = element_text(margin = margin(t = 15)),  
        axis.text.y = element_text(size = 14, color = "black"), axis.title.y = element_markdown(margin = margin(r = 15)), panel.grid.major.x = element_blank(),
        legend.position = "none") +
  labs(title = "Broekemahoeve", x = "Year", y = "Number of *Cotesia glomerata* eggs per parasitized caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  geom_vline(xintercept = 2.5, linetype = 2) + 
  ylim(c(0,48))
dev.off()


# Subsection 6.4: Rivera vs Christmas Drumhead
png(filename = "CgE - Cultivars - Wheat.png", width = 2060, height = 2000, res = 400)
ggplot(data = CLD.CM1.1, aes(x = wheat, y = response)) +
  geom_pointrange(aes(x = wheat, y = response, ymin = lower.CL, ymax = upper.CL, 
                      size = 1,
                      color = cultivar, group = cultivar), shape = 16, position = position_dodge(width = 0.75)) +
  geom_text(aes(label = .group, group = cultivar, y = upper.CL + 2, size = 2) , position = position_dodge(width = 0.75)) +
  scale_color_manual(values = Cultivar.colors) +
  theme_minimal() +
  theme(title = element_text(size = 14), axis.title = element_text(size = 16), 
        axis.text.x = element_text(size = 16, color = "black"), axis.title.x = element_text(margin = margin(t = 15)),  
        axis.text.y = element_text(size = 14, color = "black"), axis.title.y = element_markdown(margin = margin(r = 15)), panel.grid.major.x = element_blank(),
        legend.position = "none") +
  labs(x = "Wheat harvest", y = "Number of *Cotesia glomerata* eggs per parasitized caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  ylim(c(0,48))
dev.off()

png(filename = "CgE - Cultivars - Year.png", width = 3100, height = 2000, res = 400)
ggplot(data = CLD.CM1.2, aes(x = year, y = response)) +
  geom_pointrange(aes(x = year, y = response, ymin = lower.CL, ymax = upper.CL, 
                      size = 1,
                      color = cultivar, group = cultivar), shape = 16, position = position_dodge(width = 0.75)) +
  geom_text(aes(label = .group, group = cultivar, y = upper.CL + 2, size = 2) , position = position_dodge(width = 0.75)) +
  scale_color_manual(values = Cultivar.colors) +
  theme_minimal() +
  theme(title = element_text(size = 14), axis.title = element_text(size = 16), 
        axis.text.x = element_text(size = 16, color = "black"), axis.title.x = element_text(margin = margin(t = 15)),  
        axis.text.y = element_text(size = 14, color = "black"), axis.title.y = element_markdown(margin = margin(r = 15)), panel.grid.major.x = element_blank(),
        legend.position = "none") +
  labs(x = "Year", y = "Number of *Cotesia glomerata* eggs per parasitized caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  geom_vline(xintercept = 2.5, linetype = 2) + 
  ylim(c(0,48))
dev.off()

