#  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", "DHARMa", "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","wheat","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$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 unfound caterpillars
PR$cg <- ifelse(PR$cg > 1, 1, PR$cg)

CgP <- aggregate(cg ~ ID + round + year + wheat + location + treatment + cultivar + field, data = PR, sum)
Pres <- aggregate(presence ~ ID + round, data = PR, sum)
CgP <- merge(CgP, Pres)

CgP <- CgP[CgP$presence > 0,]
CgP$nocg <- CgP$presence - CgP$cg


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

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

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

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

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

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



# Section 4: Apply the models in GLM----
yT <- cbind(CgPT$cg, CgPT$nocg)
yT.droef <- cbind(CgPT.droef$cg, CgPT.droef$nocg)
yT.broek <- cbind(CgPT.broek$cg, CgPT.broek$nocg)
yC <- cbind(CgPC$cg, CgPC$nocg)

# Subsection 4.1: both locations (excl strip_rot and pixel)
TM1 <- glmmTMB(yT ~ 1 + treatment +
                 location + wheat + year + 
                 (1|round*location) + (1|field) +
                 treatment*(location + year) +
                 location*(wheat),
               family = binomial,
               ziformula =~ 1,
               data = CgPT, na.action = na.fail)

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

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

# Subsection 4.4: Rivera vs Christmas Drumhead 
CM1 <- glmmTMB(yC ~ 1 + cultivar +
                 location + wheat + year + 
                 (1|round) + (1|field) +
                 cultivar*(year + location) + 
                 location*(year),
               family = binomial,
               ziformula =~ 1,
               data = CgPC, 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: Rivera vs Christmas Drumhead 
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 = "CgPR.cat - Both locations - Wheat.png", width = 2060, height = 2000, res = 400)
ggplot(data = CLD.TM1.1, aes(x = wheat, y = prob)) +
  geom_pointrange(aes(x = wheat, y = prob, 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 + 0.05, 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 = 18), 
        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 = "*Cotesia glomerata* parasitism rate per caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  ylim(c(0,1.05))
dev.off()

png(filename = "CgPR.cat - Both locations - Year.png", width = 3100, height = 2000, res = 400)
ggplot(data = CLD.TM1.2, aes(x = year, y = prob)) +
  geom_pointrange(aes(x = year, y = prob, 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 + 0.05, 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 = 18), 
        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 = "*Cotesia glomerata* parasitism rate per 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,1.05))
dev.off()

# Subsection 6.2: All treatments / Droevendaal
png(filename = "CgPR.cat - Droevendaal - Wheat.png", width = 2060, height = 2000, res = 400)
ggplot(data = CLD.TM2.1, aes(x = wheat, y = prob)) +
  geom_pointrange(aes(x = wheat, y = prob, 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 + 0.05, 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 = 18), 
        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 = "*Cotesia glomerata* parasitism rate per caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  ylim(c(0,1.05))
dev.off()

png(filename = "CgPR.cat - Droevendaal - Year.png", width = 3100, height = 2000, res = 400)
ggplot(data = CLD.TM2.2, aes(x = year, y = prob)) +
  geom_pointrange(aes(x = year, y = prob, 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 + 0.05, 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 = 18), 
        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 = "*Cotesia glomerata* parasitism rate per 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.01,1.07))
dev.off()


# Subsection 6.3: All treatments / Broekemahoeve
png(filename = "CgPR.cat - Broekemahoeve - Wheat.png", width = 2060, height = 2000, res = 400)
ggplot(data = CLD.TM3.1, aes(x = wheat, y = prob)) +
  geom_pointrange(aes(x = wheat, y = prob, 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 + 0.05, 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 = 18), 
        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 = "*Cotesia glomerata* parasitism rate per caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  ylim(c(0,1.12))
dev.off()

png(filename = "CgPR.cat - Broekemahoeve - Year.png", width = 3100, height = 2000, res = 400)
ggplot(data = CLD.TM3.2, aes(x = year, y = prob)) +
  geom_pointrange(aes(x = year, y = prob, 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 + 0.05, 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 = 18), 
        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 = "*Cotesia glomerata* parasitism rate per 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.14,1.23))
dev.off()


# Subsection 6.4: Rivera vs Christmas Drumhead
png(filename = "CgPR.cat - Cultivars - Wheat.png", width = 2060, height = 2000, res = 400)
ggplot(data = CLD.CM1.1, aes(x = wheat, y = prob)) +
  geom_pointrange(aes(x = wheat, y = prob, 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 + 0.05, 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 = 18), 
        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 = "*Cotesia glomerata* parasitism rate per caterpillar") +
  geom_vline(xintercept = 0.4, linetype = 1) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1.5, linetype = 2) +
  ylim(c(0,1.05))
dev.off()

png(filename = "CgPR.cat - Cultivars - Year.png", width = 3100, height = 2000, res = 400)
ggplot(data = CLD.CM1.2, aes(x = year, y = prob)) +
  geom_pointrange(aes(x = year, y = prob, 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 + 0.05, 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 = 18), 
        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 = "*Cotesia glomerata* parasitism rate per 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,1.14))
dev.off()
