# GENERAL INFORMATION  ##########

# This R script contains all statistical analyses presented in the manuscript
#"Service crops for protection against multiple insect pests in sugar beet during crop establishment"
# by Neus Rodriguez-Gasol, Göran Bergkvist & Ola Lundin. 

#For a detailed description of sampling protocols and inclusion criteria, please see the main text.

# All analyses were performed in R version 4.4.2 (2024-10-31 ucrt) on x86_64, mingw32 system.

# Creative Commons CC0 1.0 Universal license


# 1. Set working directory ----
#setwd()
#getwd()



# 2. Load packages ----
library(dplyr)
library(lmerTest)   # for LM
library(car)        # for ANOVA
library(emmeans)    # for posthoc
library(multcomp)   # for cld
library(DHARMa)     # Check model assumptions
library(glmmTMB)    # GLMM
library(ggplot2)    # for plots
library(patchwork)  # to create combined figure
library(purrr)
library(rstatix)
library(tidyr)




# 3. Load datasets & data handling ----
##    3.1. Leaf pest assessment 2024 ---- 

##Load dataset
pd24 <- read.csv("Leaf_pest_assessment_2024.csv", header=T, sep=";", na.strings="NA") 
pd24$Year       <- as.factor(pd24$Year)
pd24$Site       <- as.factor(pd24$Site)
pd24$Date <- as.factor(pd24$Date)
pd24$Survey <- as.factor(pd24$Survey)
pd24$Block <- as.factor(pd24$Block)
pd24$Treatment <- as.factor(pd24$Treatment)
pd24$Treatment_name <- as.factor(pd24$Treatment_name)

##exclude 10 extra plants for dates in which root assessment was also done (plant 21 to 30)
pd24 <- pd24 %>%
  filter(!(Plant >= 21 & Plant <= 30))

View(pd24)
str(pd24)


##    3.2. Leaf pest assessment 2025 ---- 

##Load dataset
pd25 <- read.csv("Leaf_pest_assessment_2025.csv", header=T, sep=";", na.strings="NA") 
pd25$Year       <- as.factor(pd25$Year)
pd25$Site       <- as.factor(pd25$Site)
pd25$Date <- as.factor(pd25$Date)
pd25$Survey <- as.factor(pd25$Survey)
pd25$Block <- as.factor(pd25$Block)
pd25$Treatment <- as.factor(pd25$Treatment)
pd25$Treatment_name <- as.factor(pd25$Treatment_name)

View(pd25)
str(pd25)


###    3.3. Merge both files into a single one for overall assessments ---- 

##Create new variables for 2025 dataset (overall plant assessment instead of true leaves or cotyledons)
pd25$Thrips_damage_nb <- rowSums(pd25[, c("Thrips_damage_cotyledon_nb", "Thrips_damage_leaf_nb")], na.rm = TRUE)
pd25$Fleabeetle_damage_nb <- rowSums(pd25[, c("Fleabeetle_damage_cotyledon_nb", "Fleabeetle_damage_leaf_nb")], na.rm = TRUE)
pd25$Mangold_beetle_damage_nb <- rowSums(pd25[, c("Mangold_beetle_damage_cotyledon_nb", "Mangold_beetle_damage_leaf_nb")], na.rm = TRUE)
pd25$Fly_eggs_nb <- rowSums(pd25[, c("Fly_eggs_cotyledon_nb", "Fly_eggs_leaf_nb")], na.rm = TRUE)
pd25$Vertebrate_damage_nb <- rowSums(pd25[, c("Vertebrate_damage_cotyledon_nb", "Vertebrate_damage_leaf_nb")], na.rm = TRUE)


##Combine 2024 and 2025 files
pd <- bind_rows(
  pd24 %>% mutate(.source = "pd24"),
  pd25 %>% mutate(.source = "pd25")
)


View(pd)
str(pd)




##    3.4. Plant assessment 2024 ---- 

##Load dataset
pa24 <- read.csv("Plant_assessment_2024.csv", header=T, sep=";", na.strings="NA") 
pa24$Year       <- as.factor(pa24$Year)
pa24$Site       <- as.factor(pa24$Site)
pa24$Date <- as.factor(pa24$Date)
pa24$Block <- as.factor(pa24$Block)
pa24$Treatment <- as.factor(pa24$Treatment)
pa24$Treatment_name <- as.factor(pa24$Treatment_name)
pa24$Quadrat <- as.factor(pa24$Quadrat)

View(pa24)
str(pa24)





##    3.5. Plant assessment 2025 ---- 

##Load dataset
pa25 <- read.csv("Plant_assessment_2025.csv", header=T, sep=";", na.strings="NA") 
pa25$Year       <- as.factor(pa25$Year)
pa25$Site       <- as.factor(pa25$Site)
pa25$Date <- as.factor(pa25$Date)
pa25$Block <- as.factor(pa25$Block)
pa25$Treatment <- as.factor(pa25$Treatment)
pa25$Treatment_name <- as.factor(pa25$Treatment_name)
pa25$Quadrat <- as.factor(pa25$Quadrat)

View(pa25)
str(pa25)


###    3.6. Merge both files into a single one for overall assessments ---- 

##Combine 2024 and 2025 files
pa <- bind_rows(
  pa24 %>% mutate(.source = "pa24"),
  pa25 %>% mutate(.source = "pa25")
)


View(pa)
str(pa)




##    3.7. Stem pest assessment 2024---- 
##Load dataset
rd24 <- read.csv("Stem_pest_assessment_2024.csv", header=T, sep=";", na.strings="NA") 
rd24$Year       <- as.factor(rd24$Year)
rd24$Site       <- as.factor(rd24$Site)
rd24$Date <- as.factor(rd24$Date)
rd24$Week <- as.factor(rd24$Week)
rd24$Block <- as.factor(rd24$Block)
rd24$Treatment <- as.factor(rd24$Treatment)
rd24$Treatment_name <- as.factor(rd24$Treatment_name)
rd24$Quadrat <- as.factor(rd24$Plant)

View(rd24)
str(rd24)






##    3.8. Stem damage 2025---- 
##Load dataset
rd25 <- read.csv("Stem_pest_assessment_2025.csv", header=T, sep=";", na.strings="NA") 
rd25$Year       <- as.factor(rd25$Year)
rd25$Date <- as.factor(rd25$Date)
rd25$Site       <- as.factor(rd25$Site)
rd25$Block <- as.factor(rd25$Block)
rd25$Treatment <- as.factor(rd25$Treatment)

View(rd25)
str(rd25)







# 4. Analysis - Leaf pest assessment ----
## 4.1 GLMM - summs ---- 
### Calculate overall sums (across all survey dates and plants)---- 
pd_summed <- pd %>%
  group_by(Year, Site, Treatment, Treatment_name, Block) %>%
  summarise(
    n_plants = n(),  # number of plants contributing to each group
    Thrips_damage_nb = sum(Thrips_damage_nb, na.rm = TRUE),
    Fleabeetle_damage_nb = sum(Fleabeetle_damage_nb, na.rm = TRUE),
    Mangold_beetle_damage_nb = sum(Mangold_beetle_damage_nb, na.rm = TRUE),
    Fly_eggs_nb = sum(Fly_eggs_nb, na.rm = TRUE),
    Vertebrate_damage_nb = sum(Vertebrate_damage_nb, na.rm = TRUE),
    .groups = "drop"
  )

View(pd_summed)
str(pd_summed)


### Exclude treatments 2 and 3 in 2024 due to low biomass---- 
pd_summed_filtered <- pd_summed %>%
  filter(!(Year == 2024 & Treatment %in% c("2", "3")))
View(pd_summed_filtered)



###Proper models---- 
### A. Thrips---- 
options(contrasts = c("contr.sum", "contr.poly"))
thrips_poisson <- glmmTMB(Thrips_damage_nb ~ Treatment * Year + (1 | Site/Block) + offset(log(n_plants)), 
                          data = pd_summed_filtered, 
                          family = poisson(),
                          offset = log(n_plants))

options(contrasts = c("contr.sum", "contr.poly"))
thrips_nbm <- glmmTMB(Thrips_damage_nb ~ Treatment * Year + (1 | Site/Block) + offset(log(n_plants)),  
                      data = pd_summed_filtered, 
                      family = nbinom2())


#### Check model assumptions with DHARMa
#poisson
sim <- simulateResiduals(thrips_poisson, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)

#negative binomial
sim <- simulateResiduals(thrips_nbm, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)



##Model signification
summary(thrips_nbm)
Anova(thrips_nbm, type = "II")


## Emmeans by treatment and year
emm_thrips_nbm <- emmeans(thrips_nbm, ~ Treatment | Year, type = "response")
summary(emm_thrips_nbm)
pairs(emm_thrips_nbm, adjust = "tukey")

cld_thrips_nbm <- multcomp::cld(emm_thrips_nbm, adjust = "tukey", Letters = letters, reversed = TRUE)
cld_thrips_nbm


#To get descriptions by year
pd_summed_filtered |>
  dplyr::group_by(Year) |>
  dplyr::summarise(mean_thrips = mean(Thrips_damage_nb / n_plants, na.rm = TRUE))



### B. Flea beetle---- 
options(contrasts = c("contr.sum", "contr.poly"))
Fleabeetle_poisson <- glmmTMB(Fleabeetle_damage_nb ~ Treatment + Year + (1 | Site/Block) + offset(log(n_plants)), 
                          data = pd_summed_filtered, 
                          family = poisson())


options(contrasts = c("contr.sum", "contr.poly"))
Fleabeetle_nbm <- glmmTMB(Fleabeetle_damage_nb ~ Treatment + Year + (1 | Site/Block) + offset(log(n_plants)), 
                      data = pd_summed_filtered, 
                      family = nbinom2())


#### Check model assumptions with DHARMa
#poisson
sim <- simulateResiduals(Fleabeetle_poisson, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)

#negative binomial
sim <- simulateResiduals(Fleabeetle_nbm, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)



##Model signification
summary(Fleabeetle_nbm)
Anova(Fleabeetle_nbm, type = "II")

## Emmeans by treatment
emm_Fleabeetle_nbm <- emmeans(Fleabeetle_nbm, ~ Treatment, type = "response")
summary(emm_Fleabeetle_nbm)
pairs(emm_Fleabeetle_nbm, adjust = "tukey")

cld_Fleabeetle_nbm <- multcomp::cld(emm_Fleabeetle_nbm, adjust = "tukey", Letters = letters, reversed = TRUE)
cld_Fleabeetle_nbm





### C. Fly eggs---- 
options(contrasts = c("contr.sum", "contr.poly"))
fly_eggs_poisson <- glmmTMB(Fly_eggs_nb ~ Treatment + Year + (1 | Site/Block) + offset(log(n_plants)), 
                         data = pd_summed_filtered, 
                         family = poisson())

options(contrasts = c("contr.sum", "contr.poly"))
fly_eggs_nbm <- glmmTMB(Fly_eggs_nb ~ Treatment + Year + (1 | Site/Block) + offset(log(n_plants)), 
                            data = pd_summed_filtered, 
                            family = nbinom2())



#### Check model assumptions with DHARMa
#poisson
sim <- simulateResiduals(fly_eggs_poisson, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)

#negative binomial
sim <- simulateResiduals(fly_eggs_nbm, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)



##Model signification
summary(fly_eggs_nbm)
Anova(fly_eggs_nbm, type = "II")

## Emmeans by treatment
emm_fly_eggs_nbm <- emmeans(fly_eggs_nbm, ~ Treatment, type = "response")
summary(emm_fly_eggs_nbm)
pairs(emm_fly_eggs_nbm, adjust = "tukey")

cld_fly_eggs_nbm <- multcomp::cld(emm_fly_eggs_nbm, adjust = "tukey", Letters = letters, reversed = TRUE)
cld_fly_eggs_nbm


## To get descriptions by year
emm_fly_eggs_nbm_year <- emmeans(fly_eggs_nbm, ~ Year, type = "response")
summary(emm_fly_eggs_nbm_year)
pairs(emm_fly_eggs_nbm_year, adjust = "tukey")



### D. Vertebrate---- 
options(contrasts = c("contr.sum", "contr.poly"))
Vertebrate_poisson <- glmmTMB(Vertebrate_damage_nb ~ Treatment * Year + (1 | Site/Block) + offset(log(n_plants)), 
                      data = pd_summed_filtered, 
                      family = poisson())

options(contrasts = c("contr.sum", "contr.poly"))
Vertebrate_nbm <- glmmTMB(Vertebrate_damage_nb ~ Treatment * Year + (1 | Site/Block) + offset(log(n_plants)), 
                      data = pd_summed_filtered, 
                      family = nbinom2())

#### Check model assumptions with DHARMa
#poisson
sim <- simulateResiduals(Vertebrate_poisson, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)

#negative binomial
sim <- simulateResiduals(Vertebrate_nbm, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)



##Model signification
summary(Vertebrate_nbm)
Anova(Vertebrate_nbm, type = "II")

## Emmeans by treatment and year
emm_Vertebrate_nbm <- emmeans(Vertebrate_nbm, ~ Treatment | Year, type = "response" )
summary(emm_thrips_nbm)
pairs(emm_Vertebrate_nbm, adjust = "tukey")

cld_Vertebrate_nbm <- multcomp::cld(emm_Vertebrate_nbm, adjust = "tukey", Letters = letters, reversed = TRUE)
cld_Vertebrate_nbm






## 4.2 Plots (EMM + CI and partial residuals) ----

#### Define labels, colors and order before plotting---- 
custom_colors <- c(
  "1" = "#007439",
  "2" = "purple4",
  "3" = "#E2D200FF",
  "4" = "#E58601FF",
  "5" = "firebrick4"
)

custom_labels <- c(
  "1" = "Barley",
  "2" = "Clover",
  "3" = "Mustard",
  "4" = "Control\nband",
  "5" = "Control\nbroad"
)

order_levels <- c("1", "2", "3", "4", "5")



### A. Thrips---- 

# 1) Model frame + extract numeric offset (then back to n_plants)
md_thrips   <- model.frame(thrips_nbm)        # data actually used by model
off_thrips  <- model.offset(md_thrips)        # numeric vector of log(n_plants)
stopifnot(!is.null(off_thrips))
np_thrips   <- exp(off_thrips)                # recover n_plants from offset

# 2) Estimated marginal means BY YEAR (per plant, offset = 0)
rg_thrips <- ref_grid(thrips_nbm, offset = 0)

emm_thrips_by_year <- emmeans(
  rg_thrips,
  ~ Treatment | Year,
  type = "response",
  weights = "proportional"
)

emm_thrips_df <- as.data.frame(emm_thrips_by_year)

# Normalize CI column names if emmeans returned asymp.LCL/UCL
if (!("lower.CL" %in% names(emm_thrips_df)) && "asymp.LCL" %in% names(emm_thrips_df)) {
  emm_thrips_df <- emm_thrips_df |>
    dplyr::rename(lower.CL = asymp.LCL, upper.CL = asymp.UCL)
}

# Enforce treatment order and factor Year
emm_thrips_df$Treatment <- factor(as.character(emm_thrips_df$Treatment),
                                  levels = order_levels)
emm_thrips_df$Year <- as.factor(emm_thrips_df$Year)

# 3) Partial residuals per plant
fitted_tot_thrips <- predict(thrips_nbm, type = "response")                  # includes offset + RE
fixed_tot_thrips  <- exp(predict(thrips_nbm, type = "link", re.form = NA))   # includes offset, no RE
obs_tot_thrips    <- md_thrips$Thrips_damage_nb

stopifnot(length(obs_tot_thrips) == length(fitted_tot_thrips),
          length(fixed_tot_thrips) == length(fitted_tot_thrips),
          length(np_thrips) == length(fitted_tot_thrips),
          all(np_thrips > 0))

partial_resid_thrips_pp <- (obs_tot_thrips / np_thrips) -
  (fitted_tot_thrips / np_thrips) +
  (fixed_tot_thrips / np_thrips)

pd_plot_thrips <- md_thrips |>
  dplyr::mutate(
    Treatment = factor(as.character(Treatment), levels = order_levels),
    Year      = as.factor(Year),
    partial_resid_thrips = partial_resid_thrips_pp
  )



# Order within the dodge: 2024 left, 2025 right (if present)
emm_thrips_df$Year <- forcats::fct_relevel(emm_thrips_df$Year, "2024", "2025")
pd_plot_thrips$Year <- forcats::fct_relevel(pd_plot_thrips$Year, "2024", "2025")

# 5) Plot (Treatment on x; dodged by Year)
dodge_w <- 0.60
pos_dodge <- position_dodge(width = dodge_w)
shape_vals <- c("2024" = 22, "2025" = 24)  # both are filled shapes


thrips_plot_by_year <- ggplot() +
  # Partial residuals (dodged by Year; shape aids legend)
  geom_jitter(
    data = pd_plot_thrips,
    aes(x = Treatment,
        y = partial_resid_thrips,
        color = Treatment,
        shape = Year,
        group = Year),
    position = position_jitterdodge(
      jitter.width = 0.10, jitter.height = 0,
      dodge.width = dodge_w
    ),
    size = 2, stroke = 0.3, fill = "white", alpha = 0.7
  ) +
  # EMM 95% CI (dodged by Year)
  geom_errorbar(
    data = dplyr::filter(emm_thrips_df, is.finite(response)),
    aes(x = Treatment, ymin = lower.CL, ymax = upper.CL,
        color = Treatment, group = Year),
    position = pos_dodge,
    width = 0.15, linewidth = 0.6
  ) +
  # EMM points (dodged by Year)
  geom_point(
    data = dplyr::filter(emm_thrips_df, is.finite(response)),
    aes(x = Treatment, y = response,
        color = Treatment, fill = Treatment,  # <-- fill mapped
        shape = Year, group = Year),
    position = pos_dodge,
    size = 3.8, stroke = 0.6
  ) +
  # Aesthetics
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_color_manual(values = custom_colors, breaks = order_levels,
                     labels = custom_labels, guide = "none") +
  scale_fill_manual(values = custom_colors, breaks = order_levels,  # <-- fill palette
                    labels = custom_labels, guide = "none") +
  scale_shape_manual(values = shape_vals) +
  labs(x = NULL, y = "Thrips injuries per plant") +
  coord_cartesian(clip = "off") +
  theme_classic(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    legend.position = "none",
    plot.margin = margin(t = 14, r = 10, b = 10, l = 10)
  )

print(thrips_plot_by_year)






### B. Flea beetle---- 
# 1) Model frame + extract numeric offset (then back to n_plants)
md_Fleabeetle   <- model.frame(Fleabeetle_nbm)        # data actually used by model
off_Fleabeetle  <- model.offset(md_Fleabeetle)        # numeric vector of log(n_plants)
stopifnot(!is.null(off_Fleabeetle))
np_Fleabeetle   <- exp(off_Fleabeetle)                # recover n_plants from offset

# 2) Estimated marginal means (per plant, offset = 0)
rg_Fleabeetle <- ref_grid(Fleabeetle_nbm, offset = 0)
emm_Fleabeetle_df <- as.data.frame(emmeans(rg_Fleabeetle, ~ Treatment, type = "response", weights = "proportional"))

# Rename columns for consistency (GLMMs usually have asymp.LCL/UCL)
if (!("lower.CL" %in% names(emm_Fleabeetle_df)) && "asymp.LCL" %in% names(emm_Fleabeetle_df)) {
  emm_Fleabeetle_df <- emm_Fleabeetle_df %>%
    rename(lower.CL = asymp.LCL, upper.CL = asymp.UCL)
}


# enforce treatment order in all data frames
emm_Fleabeetle_df$Treatment <- factor(as.character(emm_Fleabeetle_df$Treatment),
                                      levels = order_levels)
md_Fleabeetle$Treatment     <- factor(as.character(md_Fleabeetle$Treatment),
                                      levels = order_levels)



# 3) Partial residuals per plant
fitted_tot_Fleabeetle <- predict(Fleabeetle_nbm, type = "response")                  # includes offset + RE
fixed_tot_Fleabeetle  <- exp(predict(Fleabeetle_nbm, type = "link", re.form = NA))   # includes offset, no RE
obs_tot_Fleabeetle    <- md_Fleabeetle$Fleabeetle_damage_nb

stopifnot(length(obs_tot_Fleabeetle) == length(fitted_tot_Fleabeetle),
          length(fixed_tot_Fleabeetle) == length(fitted_tot_Fleabeetle),
          length(np_Fleabeetle) == length(fitted_tot_Fleabeetle),
          all(np_Fleabeetle > 0))

partial_resid_Fleabeetle_pp <- (obs_tot_Fleabeetle / np_Fleabeetle) -
  (fitted_tot_Fleabeetle / np_Fleabeetle) +
  (fixed_tot_Fleabeetle / np_Fleabeetle)

pd_plot_Fleabeetle <- md_Fleabeetle %>%
  mutate(
    Treatment = factor(Treatment),
    partial_resid_Fleabeetle = partial_resid_Fleabeetle_pp
  )



# 4) Plot (per-plant scale for everything)
Fleabeetle_plot <- ggplot() +
  # partial residuals (per plant)
  geom_jitter(
    data = pd_plot_Fleabeetle,
    aes(x = Treatment, y = partial_resid_Fleabeetle, color = Treatment),
    width = 0.1, size = 2, shape = 21, fill = "white", alpha = 0.7
  ) +
  # EMMs + 95% CI (per plant)
  geom_errorbar(
    data = emm_Fleabeetle_df,
    aes(x = Treatment, ymin = lower.CL, ymax = upper.CL, color = Treatment),
    width = 0.15, linewidth = 0.6
  ) +
  geom_point(
    data = emm_Fleabeetle_df,
    aes(x = Treatment, y = response, color = Treatment),
    size = 4
  ) +
  # your order & display names for the x-axis
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_color_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  scale_fill_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  labs(x = NULL, y = "Flea beetle injuries per plant") +
  coord_cartesian(clip = "off") +
  theme_classic(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    legend.position = "none",
    plot.margin = margin(t = 14, r = 10, b = 10, l = 10)
  )

print(Fleabeetle_plot)






### C. Fly eggs---- 
# 1) Model frame + extract numeric offset (then back to n_plants)
md_Fly_eggs   <- model.frame(fly_eggs_nbm)        # data actually used by model
off_Fly_eggs  <- model.offset(md_Fly_eggs)        # numeric vector of log(n_plants)
stopifnot(!is.null(off_Fly_eggs))
np_Fly_eggs   <- exp(off_Fly_eggs)                # recover n_plants from offset

# 2) Estimated marginal means (per plant, offset = 0)
rg_Fly_eggs <- ref_grid(fly_eggs_nbm, offset = 0)
emm_Fly_eggs_df <- as.data.frame(emmeans(rg_Fly_eggs, ~ Treatment, type = "response", weights = "proportional"))

# Rename columns for consistency (GLMMs usually have asymp.LCL/UCL)
if (!("lower.CL" %in% names(emm_Fly_eggs_df)) && "asymp.LCL" %in% names(emm_Fly_eggs_df)) {
  emm_Fly_eggs_df <- emm_Fly_eggs_df %>%
    rename(lower.CL = asymp.LCL, upper.CL = asymp.UCL)
}


# enforce treatment order in all data frames
emm_Fly_eggs_df$Treatment <- factor(as.character(emm_Fly_eggs_df$Treatment),
                                  levels = order_levels)
md_Fly_eggs$Treatment     <- factor(as.character(md_Fly_eggs$Treatment),
                                  levels = order_levels)



# 3) Partial residuals per plant
fitted_tot_Fly_eggs <- predict(fly_eggs_nbm, type = "response")                  # includes offset + RE
fixed_tot_Fly_eggs  <- exp(predict(fly_eggs_nbm, type = "link", re.form = NA))   # includes offset, no RE
obs_tot_Fly_eggs    <- md_Fly_eggs$Fly_eggs_nb

stopifnot(length(obs_tot_Fly_eggs) == length(fitted_tot_Fly_eggs),
          length(fixed_tot_Fly_eggs) == length(fitted_tot_Fly_eggs),
          length(np_Fly_eggs) == length(fitted_tot_Fly_eggs),
          all(np_Fly_eggs > 0))

partial_resid_Fly_eggs_pp <- (obs_tot_Fly_eggs / np_Fly_eggs) -
  (fitted_tot_Fly_eggs / np_Fly_eggs) +
  (fixed_tot_Fly_eggs / np_Fly_eggs)

pd_plot_Fly_eggs <- md_Fly_eggs %>%
  mutate(
    Treatment = factor(Treatment),
    partial_resid_Fly_eggs = partial_resid_Fly_eggs_pp
  )


# 4) Plot (per-plant scale for everything)
Fly_eggs_plot <- ggplot() +
  # partial residuals (per plant)
  geom_jitter(
    data = pd_plot_Fly_eggs,
    aes(x = Treatment, y = partial_resid_Fly_eggs, color = Treatment),
    width = 0.1, size = 2, shape = 21, fill = "white", alpha = 0.7
  ) +
  # EMMs + 95% CI (per plant)
  geom_errorbar(
    data = emm_Fly_eggs_df,
    aes(x = Treatment, ymin = lower.CL, ymax = upper.CL, color = Treatment),
    width = 0.15, linewidth = 0.6
  ) +
  geom_point(
    data = emm_Fly_eggs_df,
    aes(x = Treatment, y = response, color = Treatment),
    size = 4
  ) +
  # your order & display names for the x-axis
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_color_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  scale_fill_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  labs(x = NULL, y = "Beet leafminer eggs per plant") +
  coord_cartesian(clip = "off") +
  theme_classic(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    legend.position = "none",
    plot.margin = margin(t = 14, r = 10, b = 10, l = 10)
  )

print(Fly_eggs_plot)




### D. Vertebrate---- 

# 1) Model frame + extract numeric offset (then back to n_plants)
md_Vertebrate   <- model.frame(Vertebrate_nbm)        # data actually used by model
off_Vertebrate  <- model.offset(md_Vertebrate)        # numeric vector of log(n_plants)
stopifnot(!is.null(off_Vertebrate))
np_Vertebrate   <- exp(off_Vertebrate)                # recover n_plants from offset

# 2) Estimated marginal means BY YEAR (per plant, offset = 0)
rg_Vertebrate <- ref_grid(Vertebrate_nbm, offset = 0)

emm_Vertebrate_by_year <- emmeans(
  rg_Vertebrate,
  ~ Treatment | Year,
  type = "response",
  weights = "proportional"
)

emm_Vertebrate_df <- as.data.frame(emm_Vertebrate_by_year)

# Normalize CI column names if emmeans returned asymp.LCL/UCL
if (!("lower.CL" %in% names(emm_Vertebrate_df)) && "asymp.LCL" %in% names(emm_Vertebrate_df)) {
  emm_Vertebrate_df <- emm_Vertebrate_df |>
    dplyr::rename(lower.CL = asymp.LCL, upper.CL = asymp.UCL)
}

# Enforce treatment order and factor Year
emm_Vertebrate_df$Treatment <- factor(as.character(emm_Vertebrate_df$Treatment),
                                  levels = order_levels)
emm_Vertebrate_df$Year <- as.factor(emm_Vertebrate_df$Year)

# 3) Partial residuals per plant
fitted_tot_Vertebrate <- predict(Vertebrate_nbm, type = "response")                  # includes offset + RE
fixed_tot_Vertebrate  <- exp(predict(Vertebrate_nbm, type = "link", re.form = NA))   # includes offset, no RE
obs_tot_Vertebrate    <- md_Vertebrate$Vertebrate_damage_nb

stopifnot(length(obs_tot_Vertebrate) == length(fitted_tot_Vertebrate),
          length(fixed_tot_Vertebrate) == length(fitted_tot_Vertebrate),
          length(np_Vertebrate) == length(fitted_tot_Vertebrate),
          all(np_Vertebrate > 0))

partial_resid_Vertebrate_pp <- (obs_tot_Vertebrate / np_Vertebrate) -
  (fitted_tot_Vertebrate / np_Vertebrate) +
  (fixed_tot_Vertebrate / np_Vertebrate)

pd_plot_Vertebrate <- md_Vertebrate |>
  dplyr::mutate(
    Treatment = factor(as.character(Treatment), levels = order_levels),
    Year      = as.factor(Year),
    partial_resid_Vertebrate = partial_resid_Vertebrate_pp
  )


# Order within the dodge: 2024 left, 2025 right (if present)
emm_Vertebrate_df$Year <- forcats::fct_relevel(emm_Vertebrate_df$Year, "2024", "2025")
pd_plot_Vertebrate$Year <- forcats::fct_relevel(pd_plot_Vertebrate$Year, "2024", "2025")

# 4) Plot (Treatment on x; dodged by Year)
dodge_w <- 0.60
pos_dodge <- position_dodge(width = dodge_w)
shape_vals <- c("2024" = 22, "2025" = 24)  # both are filled shapes


Vertebrate_plot_by_year <- ggplot() +
  # Partial residuals (dodged by Year; shape aids legend)
  geom_jitter(
    data = pd_plot_Vertebrate,
    aes(x = Treatment,
        y = partial_resid_Vertebrate,
        color = Treatment,
        shape = Year,
        group = Year),
    position = position_jitterdodge(
      jitter.width = 0.10, jitter.height = 0,
      dodge.width = dodge_w
    ),
    size = 2, stroke = 0.3, fill = "white", alpha = 0.7
  ) +
  # EMM 95% CI (dodged by Year)
  geom_errorbar(
    data = dplyr::filter(emm_Vertebrate_df, is.finite(response)),
    aes(x = Treatment, ymin = lower.CL, ymax = upper.CL,
        color = Treatment, group = Year),
    position = pos_dodge,
    width = 0.15, linewidth = 0.6
  ) +
  # EMM points (dodged by Year)
  geom_point(
    data = dplyr::filter(emm_Vertebrate_df, is.finite(response)),
    aes(x = Treatment, y = response,
        color = Treatment, fill = Treatment,  # <-- fill mapped
        shape = Year, group = Year),
    position = pos_dodge,
    size = 3.8, stroke = 0.6
  ) +
  # Aesthetics
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_color_manual(values = custom_colors, breaks = order_levels,
                     labels = custom_labels, guide = "none") +
  scale_fill_manual(values = custom_colors, breaks = order_levels,  # <-- fill palette
                    labels = custom_labels, guide = "none") +
  scale_shape_manual(values = shape_vals) +
  labs(x = NULL, y = "Vertebrate injuries per plant") +
  coord_cartesian(clip = "off") +
  theme_classic(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    legend.position = "none",
    plot.margin = margin(t = 14, r = 10, b = 10, l = 10)
  )

print(Vertebrate_plot_by_year)








### Combined figure---- 
Leaf_pest_assessment <- 
  (thrips_plot_by_year + Fleabeetle_plot + Fly_eggs_plot + Vertebrate_plot_by_year)+
  plot_layout(ncol = 2, nrow= 2) +
  plot_annotation(tag_levels = "A")

Leaf_pest_assessment


#Save figure
ggsave("Leaf_pest_assessment.png", Leaf_pest_assessment, width = 10, height = 8, dpi = 300, bg = "white")








# 5. Analysis - Stem damage ----

## 5.1 2024 Presence/absence AGGREGATED per plot---- 
## binary response 
rd24_pa <- rd24 %>%
  mutate(
    damage_pa = case_when(
      is.na(Mangold_beetle_Stem) ~ NA_real_,
      Mangold_beetle_Stem > 0 ~ 1,
      TRUE ~ 0
    )
  )

## exclude treatments 2 and 3
rd24_pa_filtered <- rd24_pa %>%
  filter(!Treatment %in% c("2", "3")) %>%
  filter(!is.na(damage_pa))

##data checks
rd24_pa_filtered %>%
  count(Site, Block, Treatment, Date)

rd24_pa_filtered %>%
  count(Site, Block, Treatment) %>%
  count(n)


##aggregate plants to the experimental plot
rd24_plot <- rd24_pa_filtered %>%
  group_by(Site, Block, Treatment) %>%
  summarise(
    injured = sum(damage_pa),
    assessed = n(),
    uninjured = assessed - injured,
    .groups = "drop"
  )


View(rd24_plot)
rd24_plot

##verify
rd24_plot %>%
  count(assessed)


##make treatment a clean factor
rd24_plot <- rd24_plot %>%
  mutate(
    Treatment = droplevels(Treatment),
    Site = droplevels(Site),
    Block = droplevels(Block)
  )

## model
Stem_pest_assessment_24 <- glmmTMB(
  cbind(injured, uninjured) ~ Treatment + (1 | Site:/Block),
  data = rd24_plot,
  family = binomial(link = "logit")
)

summary(Stem_pest_assessment_24)
Anova(Stem_pest_assessment_24, type = "II")


# Diagnostics
sim_24 <- simulateResiduals(Stem_pest_assessment_24, n = 1000)
plot(sim_24)
testDispersion(sim_24)

# Treatment estimates, averaged across sites/blocks
emm_24 <- emmeans(Stem_pest_assessment_24, ~ Treatment,  type = "response")
summary(emm_24)
pairs(emm_24, adjust = "tukey")

cld_24 <- multcomp::cld(emm_24, adjust = "tukey", Letters = letters,reversed = TRUE)
cld_24





## 5.2 2025 Presence/absence AGGREGATED ALL SITES -------------------------

## Transform damage to presence/absence
rd25_pa <- rd25 %>%
  mutate(
    damage_pa = case_when(
      is.na(Damage_percentage) ~ NA_real_,
      Damage_percentage > 0 ~ 1,
      TRUE ~ 0
    )
  ) %>%
  filter(!is.na(damage_pa))


## Data checks
# Check number of plants per plot
rd25_pa %>%
  count(Site, Block, Treatment)

# Check that all plots have the expected number of plants
rd25_pa %>%
  count(Site, Block, Treatment) %>%
  count(n)


## Aggregate plants to the experimental plot
rd25_plot <- rd25_pa %>%
  group_by(Site, Block, Treatment) %>%
  summarise(
    injured = sum(damage_pa),
    assessed = n(),
    uninjured = assessed - injured,
    injury_proportion = injured / assessed,
    .groups = "drop"
  )


## Inspect data
View(rd25_plot)
rd25_plot


## Binomial GLMM – all three sites
Stem_pest_assessment_25 <- glmmTMB(
  cbind(injured, uninjured) ~ Treatment + (1 | Site/Block),
  data = rd25_plot,
  family = binomial(link = "logit")
)

## Model summary
summary(Stem_pest_assessment_25)
Anova(Stem_pest_assessment_25,type = "II")


## Diagnostics
sim_25 <- simulateResiduals(Stem_pest_assessment_25,  n = 1000)
plot(sim_25)
testDispersion(sim_25)
testZeroInflation(sim_25)

## Treatment estimates
emm_25 <- emmeans(Stem_pest_assessment_25, ~ Treatment, type = "response")
summary(emm_25)
pairs(  emm_25,  adjust = "tukey")
cld_25 <- multcomp::cld(  emm_25,  adjust = "tukey",  Letters = letters,  reversed = TRUE)
cld_25




### Plot (EMM + CI) ----

##### Define labels, colors and order before plotting
custom_colors <- c(
  "1" = "#007439",
  "2" = "purple4",
  "3" = "#E2D200FF",
  "4" = "#E58601FF",
  "5" = "firebrick4"
)

custom_labels <- c(
  "1" = "Barley",
  "2" = "Clover",
  "3" = "Mustard",
  "4" = "Control\nband",
  "5" = "Control\nbroad"
)

order_levels <- c("1", "2", "3", "4", "5")


# detect no-event groups (all 0s)
counts <- rd25_plot %>%
  dplyr::group_by(Treatment) %>%
  dplyr::summarise(
    n_total  = dplyr::n(),
    n_damage = sum(injured, na.rm = TRUE),
    .groups  = "drop"
  )

## 2) EMMs table (data frame + harmonize names + flags)
# Get probabilities (type = "response") from emmeans
emm_bin_df <- summary(emm_25, type = "response") %>%  # <- response scale
  as.data.frame()

# Standardize columns for plotting
# emmeans on response scale for binomial uses 'prob', 'asymp.LCL', 'asymp.UCL'
if ("prob" %in% names(emm_bin_df)) {
  emm_bin_df <- dplyr::rename(emm_bin_df, response = prob)
}
if ("asymp.LCL" %in% names(emm_bin_df)) {
  emm_bin_df <- dplyr::rename(emm_bin_df, lower.CL = asymp.LCL, upper.CL = asymp.UCL)
}

# Enforce Treatment order and join counts; flag where to draw a CI
emm_bin_df <- emm_bin_df %>%
  dplyr::mutate(Treatment = factor(as.character(Treatment), levels = order_levels)) %>%
  dplyr::left_join(counts, by = "Treatment") %>%
  dplyr::mutate(draw_CI = !is.na(n_damage) & n_damage > 0)

# Data for the CI layer
emm_bin_ci <- emm_bin_df %>%
  dplyr::filter(draw_CI, is.finite(lower.CL), is.finite(upper.CL))

## 3) Plot
Stem_pest_assessment_plot <- ggplot() +
  # CIs (only where at least one event was observed)
  geom_errorbar(
    data = emm_bin_ci,
    aes(x = Treatment, ymin = lower.CL, ymax = upper.CL, color = Treatment),
    width = 0.15, linewidth = 0.7, na.rm = TRUE
  ) +
  # Points (always)
  geom_point(
    data = emm_bin_df,
    aes(x = Treatment, y = response, color = Treatment),
    size = 4, na.rm = TRUE
  ) +
  scale_x_discrete(limits = order_levels, labels = c("Barley","Clover","Mustard","Control\nband","Control\nbroad")) +
  scale_color_manual(values = custom_colors, breaks = order_levels, guide = "none") +
  scale_y_continuous(limits = c(0, 0.5), expand = expansion(mult = c(0.02, 0.10))) +
  labs(x = NULL, y = "Probability of stem injury") +
  coord_cartesian(clip = "off") +
  theme_classic(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    legend.position = "none",
    plot.margin = margin(t = 14, r = 10, b = 10, l = 10)
  )

print(Stem_pest_assessment_plot)


#Save figure
ggsave("Stem_pest_assessment.png", Stem_pest_assessment_plot, width = 5, height = 4, dpi = 300, bg = "white")








# 6. Analysis - Plant assessment ----

## 6.1 GLMM - sums ---- 

###  Calculate overall sums (across replicates)---- 
sum_or_na <- function(x) {
  if (all(is.na(x))) NA_real_ else sum(x, na.rm = TRUE) # important step to keep the NA as they are
}

pa_sum <- pa %>%
  group_by(Year, Site, Treatment, Treatment_name, Block) %>%
  summarise(across(where(is.numeric), sum_or_na), .groups = "drop")

View(pa_sum)


### Exclude treatments 2 and 3 in 2024 due to low biomass---- 
pa_sum_filtered <- pa_sum %>%
  filter(!(Year == 2024 & Treatment %in% c("2", "3")))
View(pa_sum_filtered)
str(pa_sum_filtered)




###Proper models---- 
### A. Sugar beet number----
options(contrasts = c("contr.sum", "contr.poly"))
Sugar_beet_nb_poisson <- glmmTMB(Sugar_beet_nb ~ Treatment + Year + (1 | Site/Block) + offset(log(Survey_area_sugar_beet_sqm)), 
                                 data = pa_sum_filtered,
                                 family = poisson())

options(contrasts = c("contr.sum", "contr.poly"))
Sugar_beet_nb_nbm <- glmmTMB(Sugar_beet_nb ~ Treatment + Year  + (1 | Site/Block) + offset(log(Survey_area_sugar_beet_sqm)), 
                             data = pa_sum_filtered,
                             family = nbinom2())
  

#### Check model assumptions with DHARMa
#poisson
sim <- simulateResiduals(Sugar_beet_nb_poisson, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)

#negative binomial
sim <- simulateResiduals(Sugar_beet_nb_nbm, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)


##Model signification
summary(Sugar_beet_nb_poisson)
Anova(Sugar_beet_nb_poisson, type = "II")

## Emmeans by treatment
emm_Sugar_beet_nb_poisson <- emmeans(Sugar_beet_nb_poisson, ~ Treatment, type = "response")
summary(emm_Sugar_beet_nb_poisson)
pairs(emm_Sugar_beet_nb_poisson, adjust = "tukey")

cld_Sugar_beet_nb_poisson <- multcomp::cld(emm_Sugar_beet_nb_poisson, adjust = "tukey", Letters = letters, reversed = TRUE)
cld_Sugar_beet_nb_poisson







### B. Sugar beet aboveground biomass (2025)----
## Exclude 2024
pa_sum_filtered_2025 <- pa_sum_filtered %>%
  filter(!(Year == 2024))
View(pa_sum_filtered_2025)


##models
#gaussian
pa_sum_filtered_2025 <- pa_sum_filtered_2025 |>
  dplyr::mutate(Sugar_beet_weight_per_m2 = Sugar_beet_weight / Survey_area_sugar_beet_sqm)

options(contrasts = c("contr.sum", "contr.poly"))
Sugar_beet_weight_gaussian <- glmmTMB(Sugar_beet_weight_per_m2 ~ Treatment + (1 | Site/Block), 
                                      data = pa_sum_filtered_2025,
                                      family = gaussian())

#gamma
options(contrasts = c("contr.sum", "contr.poly"))
Sugar_beet_weight_gamma <- glmmTMB(Sugar_beet_weight ~ Treatment  + (1 | Site/Block) + offset(log(Survey_area_sugar_beet_sqm)), 
                                   data = pa_sum_filtered_2025,
                                   family = Gamma(link = "log"))



#### Check model assumptions with DHARMa
#gaussian
sim <- simulateResiduals(Sugar_beet_weight_gaussian, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)


#gamma
sim <- simulateResiduals(Sugar_beet_weight_gamma, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)



##Model signification
summary(Sugar_beet_weight_gaussian)
Anova(Sugar_beet_weight_gaussian, type = "II")

## Emmeans by treatment
emm_Sugar_beet_weight_gaussian <- emmeans(Sugar_beet_weight_gaussian, ~ Treatment, type = "response")
pairs(emm_Sugar_beet_weight_gaussian, adjust = "tukey")
summary(emm_Sugar_beet_weight_gaussian)

cld_Sugar_beet_weight_gaussian <- multcomp::cld(emm_Sugar_beet_weight_gaussian, adjust = "tukey", Letters = letters, reversed = TRUE)
cld_Sugar_beet_weight_gaussian





### C. Service crop aboveground biomass---- 
#gaussian
pa_sum_filtered <- pa_sum_filtered |>
  dplyr::mutate(Service_crop_weight_per_m2 = Service_crop_weight / Survey_area_service_crop_sqm)


options(contrasts = c("contr.sum", "contr.poly"))
Service_crop_weight_gaussian <- glmmTMB(Service_crop_weight_per_m2 ~ Treatment + Year + (1 | Site/Block), 
                                        data = pa_sum_filtered,
                                        family = gaussian())

#gamma
options(contrasts = c("contr.sum", "contr.poly"))
Service_crop_weight_gamma <- glmmTMB(Service_crop_weight ~ Treatment + Year + (1 | Site/Block) + offset(log(Survey_area_service_crop_sqm)), 
                                     data = pa_sum_filtered,
                                     family = Gamma(link = "log"))



#### Check model assumptions with DHARMa
#gaussian
sim <- simulateResiduals(Service_crop_weight_gaussian, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)


#gamma
sim <- simulateResiduals(Service_crop_weight_gamma, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)


##Model signification
summary(Service_crop_weight_gamma)
Anova(Service_crop_weight_gamma, type = "II")

## Emmeans by treatment
emm_Service_crop_weight_gamma <- emmeans(Service_crop_weight_gamma, ~ Treatment, type = "response")
pairs(emm_Service_crop_weight_gamma, adjust = "tukey")
summary(emm_Service_crop_weight_gamma)

cld_Service_crop_weight_gamma <- multcomp::cld(emm_Service_crop_weight_gamma, adjust = "tukey", Letters = letters, reversed = TRUE)
cld_Service_crop_weight_gamma





### D. Weed aboveground biomass---- 
#gaussian
pa_sum_filtered <- pa_sum_filtered |>
  dplyr::mutate(Weed_weight_per_m2 = Weed_weight / Survey_area_weed_sqm)


options(contrasts = c("contr.sum", "contr.poly"))
Weed_weight_gaussian <- glmmTMB(Weed_weight_per_m2 ~ Treatment * Year + (1 | Site/Block), 
                                        data = pa_sum_filtered,
                                        family = gaussian())

#gaussian with log transformed variable
options(contrasts = c("contr.sum", "contr.poly"))
Weed_weight_gaussian_loga <- glmmTMB(log(Weed_weight_per_m2) ~ Treatment * Year + (1 | Site/Block), 
                                    data = pa_sum_filtered,
                                    family = gaussian())

#gamma
options(contrasts = c("contr.sum", "contr.poly"))
Weed_weight_gamma <- glmmTMB(Weed_weight ~ Treatment * Year + (1 | Site/Block) + offset(log(Survey_area_weed_sqm)), 
                             data = pa_sum_filtered,
                             family = Gamma(link = "log"))

#### Check model assumptions with DHARMa
#gaussian
sim <- simulateResiduals(Weed_weight_gaussian, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)

#gaussian with log transformed variable
sim <- simulateResiduals(Weed_weight_gaussian_loga, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)

#gamma
sim <- simulateResiduals(Weed_weight_gamma, n = 1000)
plot(sim)                       # Uniform QQ, residuals vs predicted, etc.
testDispersion(sim)             # Over/under-dispersion
testZeroInflation(sim)          # Zero inflation (for count models)


##Model signification
summary(Weed_weight_gaussian_loga)
Anova(Weed_weight_gaussian_loga, type = "II")

## Emmeans by treatment and year
emm_Weed_weight_gaussian_loga <- emmeans(Weed_weight_gaussian_loga, ~ Treatment | Year, type = "response")
summary(emm_Weed_weight_gaussian_loga)
pairs(emm_Weed_weight_gaussian_loga, adjust = "tukey")

cld_Weed_weight_gaussian_loga <- multcomp::cld(emm_Weed_weight_gaussian_loga, adjust = "tukey", Letters = letters, reversed = TRUE)
cld_Weed_weight_gaussian_loga





## 6.2 Plots (EMM + CI and partial residuals)----

##### Define labels, colors and order before plotting---- 
custom_colors <- c(
  "1" = "#007439",
  "2" = "purple4",
  "3" = "#E2D200FF",
  "4" = "#E58601FF",
  "5" = "firebrick4"
)

custom_labels <- c(
  "1" = "Barley",
  "2" = "Clover",
  "3" = "Mustard",
  "4" = "Control\nband",
  "5" = "Control\nbroad"
)

order_levels <- c("1", "2", "3", "4", "5")


### A. Sugar beet number---- 

# 1) Model frame + extract numeric offset (then back to Survey_area_sugar_beet_sqm)
md_sugar_beet_nb  <- model.frame(Sugar_beet_nb_poisson)      # data actually used by model
off_sugar_beet_nb <- model.offset(md_sugar_beet_nb)          # numeric vector of log(survey area)
stopifnot(!is.null(off_sugar_beet_nb))
np_sugar_beet_nb  <- exp(off_sugar_beet_nb)                  # recover Survey_area_sugar_beet_sqm
stopifnot(length(np_sugar_beet_nb) == nrow(md_sugar_beet_nb), all(np_sugar_beet_nb > 0, na.rm = TRUE))

# 2) Estimated marginal means (per plant, offset = 0)
rg_sugar_beet_nb <- ref_grid(Sugar_beet_nb_poisson, offset = 0)
emm_sugar_beet_nb_df <- as.data.frame(
  emmeans(rg_sugar_beet_nb, ~ Treatment, type = "response", weights = "proportional")
)

# Standardize CI column names if emmeans returns asymp.LCL/UCL
if (!("lower.CL" %in% names(emm_sugar_beet_nb_df)) && "asymp.LCL" %in% names(emm_sugar_beet_nb_df)) {
  emm_sugar_beet_nb_df <- emm_sugar_beet_nb_df %>%
    rename(lower.CL = asymp.LCL, upper.CL = asymp.UCL)
}

# Standardize the predicted column name to `response` (handles rate/prob/emmean/etc.)
resp_candidates <- c("response", "rate", "prob", "p", "emmean")
resp_col <- intersect(resp_candidates, names(emm_sugar_beet_nb_df))
if (length(resp_col) == 0) {
  stop("No response-like column found in emmeans output. Columns are: ",
       paste(names(emm_sugar_beet_nb_df), collapse = ", "))
}
emm_sugar_beet_nb_df <- emm_sugar_beet_nb_df %>%
  rename(response = !!resp_col[1])


# enforce treatment order in all data frames
emm_sugar_beet_nb_df$Treatment <- factor(as.character(emm_sugar_beet_nb_df$Treatment),
                                         levels = order_levels)
md_sugar_beet_nb$Treatment     <- factor(as.character(md_sugar_beet_nb$Treatment),
                                         levels = order_levels)



# 3) Partial residuals per plant
# fitted with RE + offset on response scale (counts)
fitted_tot_sugar_beet_nb <- predict(Sugar_beet_nb_poisson, type = "response")                  
# fixed effects only, includes offset, on response scale
fixed_tot_sugar_beet_nb  <- exp(predict(Sugar_beet_nb_poisson, type = "link", re.form = NA))   
obs_tot_sugar_beet_nb    <- md_sugar_beet_nb$Sugar_beet_nb

stopifnot(length(obs_tot_sugar_beet_nb) == length(fitted_tot_sugar_beet_nb),
          length(fixed_tot_sugar_beet_nb) == length(fitted_tot_sugar_beet_nb),
          length(np_sugar_beet_nb) == length(fitted_tot_sugar_beet_nb))

# per m² partial residuals (divide each component by survey area)
partial_resid_sugar_beet_nb_pp <- (obs_tot_sugar_beet_nb / np_sugar_beet_nb) -
  (fitted_tot_sugar_beet_nb / np_sugar_beet_nb) +
  (fixed_tot_sugar_beet_nb / np_sugar_beet_nb)

pd_plot_sugar_beet_nb <- md_sugar_beet_nb %>%
  mutate(
    Treatment = factor(Treatment, levels = order_levels),
    partial_resid_sugar_beet_nb = partial_resid_sugar_beet_nb_pp
  )

# 4) Plot (per-plant scale for everything)
sugar_beet_nb_plot <- ggplot() +
  # partial residuals (per plant)
  geom_jitter(
    data = pd_plot_sugar_beet_nb,
    aes(x = Treatment, y = partial_resid_sugar_beet_nb, color = Treatment),
    width = 0.1, size = 2, shape = 21, fill = "white", alpha = 0.7
  ) +
  # EMMs + 95% CI (per plant)
  geom_errorbar(
    data = emm_sugar_beet_nb_df,
    aes(x = Treatment, ymin = lower.CL, ymax = upper.CL, color = Treatment),
    width = 0.12, linewidth = 0.9
  ) +
  geom_point(
    data = emm_sugar_beet_nb_df,
    aes(x = Treatment, y = response, color = Treatment),
    size = 4
  ) +
  # your order & display names for the x-axis
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_color_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  scale_fill_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  labs(x = NULL, y = expression(Sugar~beet~plants~(m^{-2}))) +
  coord_cartesian(clip = "off") +
  theme_classic(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    legend.position = "none",
    plot.margin = margin(t = 14, r = 10, b = 10, l = 10)
  )

print(sugar_beet_nb_plot)






### B. Sugar beet aboveground biomass (2025)----

## 1) EMMs
emm_sugar_beet_weight_df <- as.data.frame(
  emmeans(Sugar_beet_weight_gaussian, ~ Treatment, type = "response", weights = "proportional")
)

# harmonize CI names if needed
if (!"response" %in% names(emm_sugar_beet_weight_df) && "emmean" %in% names(emm_sugar_beet_weight_df)) {
  emm_sugar_beet_weight_df$response <- emm_sugar_beet_weight_df$emmean
}

# enforce order
emm_sugar_beet_weight_df$Treatment <- factor(as.character(emm_sugar_beet_weight_df$Treatment), levels = order_levels)

## 2) Partial residuals on the response scale (identity link)
md_sugar_beet_weight <- model.frame(Sugar_beet_weight_gaussian)
md_sugar_beet_weight$Treatment <- factor(as.character(md_sugar_beet_weight$Treatment), levels = order_levels)

fitted_tot_sugar_beet_weight <- predict(Sugar_beet_weight_gaussian, type = "response")                 
fixed_tot_sugar_beet_weight  <- predict(Sugar_beet_weight_gaussian, type = "response", re.form = NA)   
obs_tot_sugar_beet_weight    <- md_sugar_beet_weight$Sugar_beet_weight_per_m2

stopifnot(length(obs_tot_sugar_beet_weight) == length(fitted_tot_sugar_beet_weight),
          length(fixed_tot_sugar_beet_weight) == length(fitted_tot_sugar_beet_weight))

partial_resid_sugar_beet_weight <- obs_tot_sugar_beet_weight - fitted_tot_sugar_beet_weight + fixed_tot_sugar_beet_weight

pd_plot_sugar_beet_weight <- dplyr::mutate(md_sugar_beet_weight, partial_resid_sugar_beet_weight = partial_resid_sugar_beet_weight)



# 3) Plot
names(custom_colors) <- order_levels

sugar_beet_weight_plot <- ggplot() +
  geom_jitter(
    data = pd_plot_sugar_beet_weight,
    aes(x = Treatment, y = partial_resid_sugar_beet_weight, color = Treatment),
    width = 0.12, size = 1.6, shape = 21, fill = "white", alpha = 0.55
  ) +
  geom_errorbar(
    data = emm_sugar_beet_weight_df,
    aes(x = Treatment, ymin = asymp.LCL, ymax = asymp.UCL, color = Treatment),
    width = 0.12, linewidth = 0.9
  ) +
  geom_point(
    data = emm_sugar_beet_weight_df,
    aes(x = Treatment, y = response, fill = Treatment, color = Treatment),
    size = 4.2, shape = 21, stroke = 0.7
  ) +
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_color_manual(values = custom_colors, guide = "none") +
  scale_fill_manual(values = custom_colors, guide = "none") +
  labs(x = NULL, y = expression(Sugar~beet~biomass~(g~m^{-2}))) +
  coord_cartesian(clip = "off") +
  scale_y_continuous(expand = expansion(mult = c(0.02, 0.12))) +
  theme_classic(base_size = 12) +
  theme(
    axis.text.x = element_text(size = 13, margin = margin(t = 6)),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    panel.grid = element_blank(),
    legend.position = "none"
  )

sugar_beet_weight_plot



### C. Service crop aboveground biomass---- 


# 1) Model frame + extract numeric offset (then back to Survey_area_sugar_beet_sqm)
md_Service_crop_weight  <- model.frame(Service_crop_weight_gamma)      # data actually used by model
off_Service_crop_weight <- model.offset(md_Service_crop_weight)          # numeric vector of log(survey area)
stopifnot(!is.null(off_Service_crop_weight))
np_Service_crop_weight  <- exp(off_Service_crop_weight)                  # recover Survey_area_sugar_beet_sqm
stopifnot(length(np_Service_crop_weight) == nrow(md_Service_crop_weight), all(np_Service_crop_weight > 0, na.rm = TRUE))

# 2) Estimated marginal means (per plant, offset = 0)
rg_Service_crop_weight <- ref_grid(Service_crop_weight_gamma, offset = 0)
emm_Service_crop_weight_df <- as.data.frame(
  emmeans(rg_Service_crop_weight, ~ Treatment, type = "response", weights = "proportional")
)

# Standardize CI column names if emmeans returns asymp.LCL/UCL
if (!("lower.CL" %in% names(emm_Service_crop_weight_df)) && "asymp.LCL" %in% names(emm_Service_crop_weight_df)) {
  emm_Service_crop_weight_df <- emm_Service_crop_weight_df %>%
    rename(lower.CL = asymp.LCL, upper.CL = asymp.UCL)
}

# Standardize the predicted column name to `response` (handles rate/prob/emmean/etc.)
resp_candidates <- c("response", "rate", "prob", "p", "emmean")
resp_col <- intersect(resp_candidates, names(emm_Service_crop_weight_df))
if (length(resp_col) == 0) {
  stop("No response-like column found in emmeans output. Columns are: ",
       paste(names(emm_Service_crop_weight_df), collapse = ", "))
}
emm_Service_crop_weight_df <- emm_Service_crop_weight_df %>%
  rename(response = !!resp_col[1])


# enforce treatment order in all data frames
emm_Service_crop_weight_df$Treatment <- factor(as.character(emm_Service_crop_weight_df$Treatment),
                                       levels = order_levels)
md_Service_crop_weight$Treatment     <- factor(as.character(md_Service_crop_weight$Treatment),
                                       levels = order_levels)



# 3) Partial residuals per plant
# fitted with RE + offset on response scale (counts)
fitted_tot_Service_crop_weight <- predict(Service_crop_weight_gamma, type = "response")                  
# fixed effects only, includes offset, on response scale
fixed_tot_Service_crop_weight  <- exp(predict(Service_crop_weight_gamma, type = "link", re.form = NA))   
obs_tot_Service_crop_weight    <- md_Service_crop_weight$Service_crop_weight

stopifnot(all(fitted_tot_Service_crop_weight > 0, na.rm = TRUE),
          all(fixed_tot_Service_crop_weight  > 0, na.rm = TRUE),
          all(np_Service_crop_weight         > 0, na.rm = TRUE))

partial_resid_Service_crop_weight_pp <- (obs_tot_Service_crop_weight / fitted_tot_Service_crop_weight) * 
  (fixed_tot_Service_crop_weight / np_Service_crop_weight)

pd_plot_Service_crop_weight <- md_Service_crop_weight %>%
  mutate(
    Treatment = factor(Treatment, levels = order_levels),
    partial_resid_Service_crop_weight = partial_resid_Service_crop_weight_pp
  )

# 4) Plot (per-plant scale for everything)
Service_crop_weight_plot <- ggplot() +
  # partial residuals (per plant)
  geom_jitter(
    data = pd_plot_Service_crop_weight,
    aes(x = Treatment, y = partial_resid_Service_crop_weight, color = Treatment),
    width = 0.1, size = 2, shape = 21, fill = "white", alpha = 0.7
  ) +
  # EMMs + 95% CI (per plant)
  geom_errorbar(
    data = emm_Service_crop_weight_df,
    aes(x = Treatment, ymin = lower.CL, ymax = upper.CL, color = Treatment),
    width = 0.12, linewidth = 0.9
  ) +
  geom_point(
    data = emm_Service_crop_weight_df,
    aes(x = Treatment, y = response, color = Treatment),
    size = 4
  ) +
  # your order & display names for the x-axis
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_color_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  scale_fill_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  labs(x = NULL, y = expression(Service~crop~biomass~(g~m^{-2}))) +
  coord_cartesian(clip = "off") +
  theme_classic(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    legend.position = "none",
    plot.margin = margin(t = 14, r = 10, b = 10, l = 10)
  )

print(Service_crop_weight_plot)




### D. Weed aboveground biomass---- 

# -- 1) EMMs on the log scale, then back-transform for plotting
emm_log <- emmeans(
  Weed_weight_gaussian_loga,
  ~ Treatment | Year,             # means for each Treatment within Year
  type = "link",                  # IMPORTANT: stay on the model (log) scale
  weights = "proportional"
)

emm_Weed_df <- as.data.frame(emm_log) %>%
  # back-transform CIs and point estimates; optionally bias-correct the mean
  mutate(
    response  = exp(emmean),          # geometric mean on original scale
    lower.CL  = exp(asymp.LCL),
    upper.CL  = exp(asymp.UCL)
  )

# (Optional) bias-correct the *mean* for log-normal: E[Y] = exp(mu + sigma^2/2)
# If you want that instead of geometric mean, uncomment the next line:
# emm_Weed_df$response <- exp(emm_Weed_df$emmean + 0.5 * sigma(Weed_weight_gaussian_loga)^2)

# enforce factor order
emm_Weed_df$Treatment <- factor(as.character(emm_Weed_df$Treatment), levels = order_levels)
emm_Weed_df$Year <- factor(emm_Weed_df$Year, levels = c("2024", "2025"))

# -- 2) Partial residuals: work on the LOG scale, then exponentiate back
md_Weed_weight <- model.frame(Weed_weight_gaussian_loga)

# Get the response actually used by the model (this is already log(Weed_weight_per_m2))
log_obs   <- model.response(md_Weed_weight)

# Predictions are on the model (log) scale for this model
fitted_log <- predict(Weed_weight_gaussian_loga, type = "response")                # with RE, log scale
fixed_log  <- predict(Weed_weight_gaussian_loga, type = "response", re.form = NA)  # fixed only, log scale

# Additive on log-scale -> multiplicative on original scale
partial_log        <- (log_obs - fitted_log) + fixed_log
partial_resid_Weed <- exp(partial_log)   # back to g m^-2

pd_plot_Weed <- md_Weed_weight %>%
  dplyr::mutate(
    Treatment = factor(as.character(Treatment), levels = order_levels),
    Year      = factor(Year, levels = c("2024","2025")),
    partial_resid_Weed_weight = partial_resid_Weed
  )


# -- 3) Plot (Treatment on x; Year dodged; EMMs filled; residuals white-filled)
dodge_w  <- 0.60
pos_d    <- position_dodge(width = dodge_w)
shape_vals <- c("2024" = 22, "2025" = 24)  # filled square & triangle

Weed_weight_plot_by_year_gaussian_loga <- ggplot() +
  # Partial residuals (multiplicative), jittered & dodged by Year
  geom_jitter(
    data = pd_plot_Weed,
    aes(x = Treatment,
        y = partial_resid_Weed_weight,
        color = Treatment,
        shape = Year,
        group = Year),
    position = position_jitterdodge(jitter.width = 0.10, jitter.height = 0,
                                    dodge.width = dodge_w),
    size = 2, stroke = 0.3, fill = "white", alpha = 0.7
  ) +
  # EMM 95% CI (on original scale), dodged by Year
  geom_errorbar(
    data = dplyr::filter(emm_Weed_df, is.finite(lower.CL), is.finite(upper.CL)),
    aes(x = Treatment, ymin = lower.CL, ymax = upper.CL,
        color = Treatment, group = Year),
    position = pos_d, width = 0.12, linewidth = 0.9
  ) +
  # EMM points (FILLED; on original scale), dodged by Year
  geom_point(
    data = dplyr::filter(emm_Weed_df, is.finite(response)),
    aes(x = Treatment, y = response,
        color = Treatment, fill = Treatment, shape = Year, group = Year),
    position = pos_d, size = 3.8, stroke = 0.6
  ) +
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_color_manual(values = custom_colors, breaks = order_levels, labels = custom_labels, guide = "none") +
  scale_fill_manual(values = custom_colors,  breaks = order_levels, labels = custom_labels, guide = "none") +
  scale_shape_manual(values = shape_vals) +
  scale_y_log10(
    labels = scales::label_number(accuracy = 0.1),
    breaks = scales::breaks_log(n = 6),
    limits = c(NA, NA)  # or set e.g. c(0.01, 100) if you want fixed range
  ) +
  labs(x = NULL, y = expression(Weed~biomass~(g~m^{-2}))) +
  coord_cartesian(clip = "off") +
  theme_classic(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.y = element_text(margin = margin(r = 8), size = 13),
    legend.position = "none",
    plot.margin = margin(t = 14, r = 10, b = 10, l = 10)
  )


print(Weed_weight_plot_by_year_gaussian_loga)



### Combined figure---- 
Plant_assessment_glm <- 
  (sugar_beet_nb_plot + sugar_beet_weight_plot + Service_crop_weight_plot +  Weed_weight_plot_by_year_gaussian_loga)+
  plot_layout(ncol = 2, nrow= 2) +
  plot_annotation(tag_levels = "A")

Plant_assessment_glm


#Save figure
ggsave("Plant_assessment.png", Plant_assessment_glm, width = 10, height = 8, dpi = 300, bg = "white")










## 6.3 Service crop plants by year and site ---- 
View(pa)



##### Define labels, colors and order before plotting---- 
custom_colors <- c(
  "1" = "#007439",
  "2" = "purple4",
  "3" = "#E2D200FF"
)

custom_labels <- c(
  "1" = "Barley",
  "2" = "Clover",
  "3" = "Mustard"
)

order_levels <- c("1", "2", "3")


# plots by Site
######Lage---- 
Service_crop_Lage <- pa %>%
  filter(Site == "Lage")



#Actual plot
Service_crop_Lage_plot <- ggplot(Service_crop_Lage, aes(x = factor(Treatment, levels = c("1", "2", "3")), y = Service_crop_nb, fill = Treatment)) +
  # main boxplot
  geom_boxplot(
    width  = 0.6,
    outlier.shape = NA,
    color  = "black",
    alpha  = 0.85
  ) +
  # add jitter points
  geom_jitter(
    width  = 0.15,
    shape  = 21,
    size   = 2.5,
    stroke = 0.3,
    color  = "black",     # border color
    alpha  = 0.8,
    aes(fill = Treatment)      # fill from group color
  ) +
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_fill_manual(
    values = custom_colors,
    breaks = order_levels,
    labels = custom_labels,
    name   = "Treatment"
  ) +
  geom_hline(yintercept = 20, linetype = "dashed", color = "grey40", linewidth = 0.7) + # 20 plants threshold
  scale_y_continuous(limits = c(0, 300)) + # fixed  scale
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line        = element_line(colour = "black", linewidth = 0.5),
    axis.text        = element_text(size = 13),
    axis.title       = element_text(size = 14, face = "bold"),
    plot.margin      = margin(5, 15, 5, 5)          # room for y-axis label
  ) +
  labs(
    x   = "Treatment",
    y   = "Service crop plants",
    fill = NULL
  ) +
  guides(fill = "none")

Service_crop_Lage_plot




######Vade---- 
Service_crop_Vade <- pa %>%
  filter(Site == "Vade")


#Actual plot
Service_crop_Vade_plot <- ggplot(Service_crop_Vade, aes(x = factor(Treatment, levels = c("1", "2", "3")), y = Service_crop_nb, fill = Treatment)) +
  # main boxplot
  geom_boxplot(
    width  = 0.6,
    outlier.shape = NA,
    color  = "black",
    alpha  = 0.85
  ) +
  # add jitter points
  geom_jitter(
    width  = 0.15,
    shape  = 21,
    size   = 2.5,
    stroke = 0.3,
    color  = "black",     # border color
    alpha  = 0.8,
    aes(fill = Treatment)      # fill from group color
  ) +
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_fill_manual(
    values = custom_colors,
    breaks = order_levels,
    labels = custom_labels,
    name   = "Treatment"
  ) +
  geom_hline(yintercept = 20, linetype = "dashed", color = "grey40", linewidth = 0.7) + # 20 plants threshold
  scale_y_continuous(limits = c(0, 300)) + # fixed  scale
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line        = element_line(colour = "black", linewidth = 0.5),
    axis.text        = element_text(size = 13),
    axis.title       = element_text(size = 14, face = "bold"),
    plot.margin      = margin(5, 15, 5, 5)          # room for y-axis label
  ) +
  labs(
    x   = "Treatment",
    y   = "Service crop plants",
    fill = NULL
  ) +
  guides(fill = "none")

Service_crop_Vade_plot



###### Adelholm---- 
Service_crop_Adelholm <- pa %>%
  filter(Site == "Adelholm")

#Actual plot
Service_crop_Adelholm_plot <- ggplot(Service_crop_Adelholm, aes(x = factor(Treatment, levels = c("1", "2", "3")), y = Service_crop_nb, fill = Treatment)) +
  # main boxplot
  geom_boxplot(
    width  = 0.6,
    outlier.shape = NA,
    color  = "black",
    alpha  = 0.85
  ) +
  # add jitter points
  geom_jitter(
    width  = 0.15,
    shape  = 21,
    size   = 2.5,
    stroke = 0.3,
    color  = "black",     # border color
    alpha  = 0.8,
    aes(fill = Treatment)      # fill from group color
  ) +
  geom_hline(yintercept = 20, linetype = "dashed", color = "grey40", linewidth = 0.7) + # 20 plants threshold
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_fill_manual(
    values = custom_colors,
    breaks = order_levels,
    labels = custom_labels,
    name   = "Treatment"
  ) +
  scale_y_continuous(limits = c(0, 300)) + # fixed  scale
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line        = element_line(colour = "black", linewidth = 0.5),
    axis.text        = element_text(size = 13),
    axis.title       = element_text(size = 14, face = "bold"),
    plot.margin      = margin(5, 15, 5, 5)          # room for y-axis label
  ) +
  labs(
    x   = "Treatment",
    y   = "Service crop plants",
    fill = NULL
  ) +
  guides(fill = "none")

Service_crop_Adelholm_plot



###### Glumslov---- 
Service_crop_Glumslov <- pa %>%
  filter(Site == "Glumslov")



#Actual plot
Service_crop_Glumslov_plot <- ggplot(Service_crop_Glumslov, aes(x = factor(Treatment, levels = c("1", "2", "3")), y = Service_crop_nb, fill = Treatment)) +
  # main boxplot
  geom_boxplot(
    width  = 0.6,
    outlier.shape = NA,
    color  = "black",
    alpha  = 0.85
  ) +
  # add jitter points
  geom_jitter(
    width  = 0.15,
    shape  = 21,
    size   = 2.5,
    stroke = 0.3,
    color  = "black",     # border color
    alpha  = 0.8,
    aes(fill = Treatment)      # fill from group color
  ) +
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_fill_manual(
    values = custom_colors,
    breaks = order_levels,
    labels = custom_labels,
    name   = "Treatment"
  ) +
  geom_hline(yintercept = 20, linetype = "dashed", color = "grey40", linewidth = 0.7) + # 20 plants threshold
  scale_y_continuous(limits = c(0, 300)) + # fixed  scale
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line        = element_line(colour = "black", linewidth = 0.5),
    axis.text        = element_text(size = 13),
    axis.title       = element_text(size = 14, face = "bold"),
    plot.margin      = margin(5, 15, 5, 5)          # room for y-axis label
  ) +
  labs(
    x   = "Treatment",
    y   = "Service crop plants",
    fill = NULL
  ) +
  guides(fill = "none")

Service_crop_Glumslov_plot



###### Lagestorp---- 
Service_crop_Lagestorp <- pa %>%
  filter(Site == "Lagestorp")





#Actual plot
Service_crop_Lagestorp_plot <- ggplot(Service_crop_Lagestorp, aes(x = factor(Treatment, levels = c("1", "2", "3")), y = Service_crop_nb, fill = Treatment)) +
  # main boxplot
  geom_boxplot(
    width  = 0.6,
    outlier.shape = NA,
    color  = "black",
    alpha  = 0.85
  ) +
  # add jitter points
  geom_jitter(
    width  = 0.15,
    shape  = 21,
    size   = 2.5,
    stroke = 0.3,
    color  = "black",     # border color
    alpha  = 0.8,
    aes(fill = Treatment)      # fill from group color
  ) +
  scale_x_discrete(limits = order_levels, labels = custom_labels) +
  scale_fill_manual(
    values = custom_colors,
    breaks = order_levels,
    labels = custom_labels,
    name   = "Treatment"
  ) +
  geom_hline(yintercept = 20, linetype = "dashed", color = "grey40", linewidth = 0.7) + # 20 plants threshold
  scale_y_continuous(limits = c(0, 300)) + # fixed  scale
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line        = element_line(colour = "black", linewidth = 0.5),
    axis.text        = element_text(size = 13),
    axis.title       = element_text(size = 14, face = "bold"),
    plot.margin      = margin(5, 15, 5, 5)          # room for y-axis label
  ) +
  labs(
    x   = "Treatment",
    y   = "Service crop plants",
    fill = NULL
  ) +
  guides(fill = "none")

Service_crop_Lagestorp_plot






###### Combined figure---- 
Service_crop_by_field <- 
  (Service_crop_Lage_plot + Service_crop_Vade_plot + Service_crop_Adelholm_plot + Service_crop_Glumslov_plot +  Service_crop_Lagestorp_plot)+
  plot_layout(ncol = 3, nrow= 2) +
  plot_annotation(tag_levels = "A")

Service_crop_by_field


#Save figure
ggsave("Service_crop_by_field.png", Service_crop_by_field, width = 10, height = 8, dpi = 300, bg = "white")