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

#This R script contains all statistical analyses presented in the manuscript "Perennial flower strips increase pollinator and natural enemy abundance but have limited effects on pest control in adjacent crops"
# 
# by Neus Rodriguez-Gasol, Maria Viketoft, Elodie Chapurlat, Johan A. Stenberg, Mattias Jonsson and 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.




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



# 2. Load packages ----

# packages  used
library(dplyr)       
library(ggplot2)     
library(glmmTMB)     
library(DHARMa)      
library(performance) 
library(car)         
library(emmeans)     
library(parameters)  
library(broom.mixed) 
library(purrr)       
library(tidyr)
library(openxlsx)
library(ggpubr)
library(lme4)
library(boot)




# 3. Load datasets & data handling ----


##    3.1. Vegetation and floral area ---- 

##Load dataset
Vegetation <- read.csv("SZ_Vegetation.csv", header=T, sep=";", na.strings="NA") 
Vegetation$Strip_type       <- as.factor(Vegetation$Strip_type)
Vegetation$Field <- as.factor(Vegetation$Field)
Vegetation$Survey <- as.factor(Vegetation$Survey)
Vegetation$Replicate <- as.factor(Vegetation$Replicate)
str(Vegetation)


# Reorder factor levels
Vegetation$Strip_type <- factor(Vegetation$Strip_type, levels = c("Flower_strip", "Control"))
View(Vegetation)


#Summarize data
Vegetation <- Vegetation %>%
  group_by(Field, Strip_type) %>%
  summarize(across(where(is.numeric), ~ sum(., na.rm = TRUE)), .groups = "drop") # Summarize data
View(Vegetation)


#Transformations
#Sum_Floralarea --> shows the sum of floral area (numbers of floral units x average floral area of each species) in the 16 quadrats of 0.6m x0.6m, so cm² per 5.76 m²
Vegetation$Floral_area_cm2_per_m2 <- Vegetation$Sum_Floralarea / 5.76 ## Convert Sum_Floralarea from cm² per 5.76m² → cm² per m²
Vegetation$Floral_area_rounded = round(Vegetation$Floral_area_cm2_per_m2) # Round floral area to avoid decimals to fit poisson

#Perc_cover --> File shows the sum, but we want the average percentage of vegetation cover in 5.76 (16 quadrats of 0.6mx0.6m)
Vegetation$Perc_cover <- Vegetation$Perc_cover / 16  # Convert Perc_cover from sum % → average % per quadrat (divide by 16)
Vegetation$Perc_cover <- Vegetation$Perc_cover / 100  # Convert percentage into proportion (0-1 scale)
View(Vegetation)



##    3.2. Pollinators ---- 

## Load dataset
Pollinators <- read.csv("SZ_Pollinators.csv", header=T, sep=";", na.strings="NA") 
Pollinators$Field    <- as.factor(Pollinators$Field)
Pollinators$Crop    <- as.factor(Pollinators$Crop)
Pollinators$Strip_type    <- as.factor(Pollinators$Strip_type)
Pollinators$Survey    <- as.factor(Pollinators$Survey)
str(Pollinators)

# Reorder factor levels
Pollinators$Strip_type <- factor(Pollinators$Strip_type, levels = c("Flower_strip", "Control"))


## Summarize data
Pollinators <- Pollinators %>%
  group_by(Field, Crop, Strip_type) %>%
  summarize(across(where(is.numeric), ~ sum(., na.rm = TRUE)), .groups = "drop")
View(Pollinators)



##    3.3. Yellow sticky traps ----  

##Load dataset
YST <- read.csv("SZ_YST.csv", header=T, sep=";", na.strings="NA") 
YST$Field    <- as.factor(YST$Field)
YST$Crop    <- as.factor(YST$Crop)
YST$Strip_type    <- as.factor(YST$Strip_type)
YST$Location   <- as.factor(YST$Location)
YST$Survey   <- as.factor(YST$Survey)
YST$Station   <- as.factor(YST$Station)
str(YST)

# Reorder factor levels
YST$Strip_type <- factor(YST$Strip_type, levels = c("Flower_strip", "Control"))
YST$Location <- factor(YST$Location, levels = c("Margin", "Crop"))


##Summarize data
YST <- YST %>%
  group_by(Field, Crop, Location, Strip_type) %>%
  summarize(across(where(is.numeric), ~ sum(., na.rm = TRUE)), .groups = "drop")
View(YST)
print(YST)


##    3.4. Pitfall traps ---- 

##Load dataset
Pitfall <- read.csv("SZ_Pitfall.csv", header=T, sep=";", na.strings="NA") 
Pitfall$Field    <- as.factor(Pitfall$Field)
Pitfall$Crop    <- as.factor(Pitfall$Crop)
Pitfall$Strip_type    <- as.factor(Pitfall$Strip_type)
Pitfall$Location   <- as.factor(Pitfall$Location)
Pitfall$Station   <- as.factor(Pitfall$Station)
Pitfall$Survey   <- as.factor(Pitfall$Survey)
str(Pitfall)

# Reorder factor levels
Pitfall$Strip_type <- factor(Pitfall$Strip_type, levels = c("Flower_strip", "Control"))
Pitfall$Location <- factor(Pitfall$Location, levels = c("Margin", "Crop"))


##Summarize data
Pitfall <- Pitfall %>%
  group_by(Field, Crop, Location, Strip_type) %>%
  summarize(across(where(is.numeric), ~ sum(., na.rm = TRUE)), .groups = "drop")
View(Pitfall)




##    3.5. Bait lamina strips ---- 

##Load dataset
BLS <- read.csv("SZ_Baitlaminastrips.csv", header=T, sep=";", na.strings="NA") 
BLS$Field <- as.factor(BLS$Field)
BLS$Crop <- as.factor(BLS$Crop)
BLS$Location <- as.factor(BLS$Location)
BLS$Strip_type <- as.factor(BLS$Strip_type)
BLS$Station <- as.factor(BLS$Station)
BLS$Lamina_code <- as.factor(BLS$Lamina_code)
str(BLS)

## Reorder factor levels
BLS$Strip_type <- factor(BLS$Strip_type, levels = c("Flower_strip", "Control"))
BLS$Location <- factor(BLS$Location, levels = c("Margin", "Crop"))

##Summarize data
BLS <- BLS %>%
  group_by(Field, Crop, Location, Strip_type) %>%
  summarize(across(where(is.numeric), ~ sum(., na.rm = TRUE)), .groups = "drop")
View(BLS)

## Create new variable (decomposition_rate)
BLS$Decomposition_rate = BLS$Number_empty_holes/BLS$Total_number_holes
View(BLS)




##    3.6. Tiller counts ---- 
##Load dataset
TC <- read.csv("SZ_Tillercounts.csv", header=T, sep=";", na.strings="NA") 
TC$Field    <- as.factor(TC$Field)
TC$Crop    <- as.factor(TC$Crop)
TC$Strip_type    <- as.factor(TC$Strip_type)
TC$Station   <- as.factor(TC$Station)
TC$Survey   <- as.factor(TC$Survey)
str(TC)

# Reorder factor levels
TC$Strip_type <- factor(TC$Strip_type, levels = c("Flower_strip", "Control"))

##Summarize data
TC <- TC %>%
  group_by(Field, Crop, Strip_type) %>%
  summarize(across(where(is.numeric), ~ sum(., na.rm = TRUE)), .groups = "drop")
View(TC)








##    3.7. Aphid cards ----
##Load dataset
AC <- read.csv("SZ_Aphidcards.csv", header=T, sep=";", na.strings="NA") 
AC$Card_type    <- as.factor(AC$Card_type)
AC$Field    <- as.factor(AC$Field)
AC$Crop    <- as.factor(AC$Crop)
AC$Strip_type    <- as.factor(AC$Strip_type)
AC$Station   <- as.factor(AC$Station)
AC$Replicate   <- as.factor(AC$Replicate)
str(AC)
View(AC)

# Reorder factor levels
AC$Strip_type <- factor(AC$Strip_type, levels = c("Flower_strip", "Control"))

##Summarize data
AC <- AC %>%
  group_by(Card_type, Field, Crop, Strip_type) %>%
  summarize(across(where(is.numeric), ~ sum(., na.rm = TRUE)), .groups = "drop")
View(AC)


####### Create new variable (predation_rate) and shrink it in order to be able to fit a beta_binomial distribution
AC$Predation_rate = (((AC$Aphids_eaten /AC$Aphids_set_up)*(AC$Aphids_set_up - 1) +0.5)/AC$Aphids_set_up)
View(AC)


##Create separate data frames for Leaf & Ground cards
#Leaf cards
AC_leaf <- AC %>% filter(Card_type == "Leaf")
AC_leaf <- AC_leaf %>%
  rename(Predation_rate_leaf = Predation_rate)
View(AC_leaf)

#Ground cards
AC_ground <- AC %>% filter(Card_type == "Ground")
AC_ground <- AC_ground %>%
  rename(Predation_rate_ground = Predation_rate)
View(AC_ground)






# 4. Models - Margin habitats ----

##    4.1. Vegetation and floral area ---- 

###    (a) Vegetation cover ---- 

####### Modify variable in order to be able to fit a beta_binomial distribution
# Substitute 1 with a slightly smaller value and divide by 100
Vegetation$Perc_cover_adj <- ifelse(Vegetation$Perc_cover == 1.000000, 0.999999, Vegetation$Perc_cover)
View(Vegetation)

####### Beta_family - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Perc_cover_beta <- glmmTMB(Perc_cover_adj ~ Strip_type + (1|Field), 
                           family = beta_family(link = "logit"), 
                           data = Vegetation)

hist(residuals(Perc_cover_beta))
plot(simulateResiduals(Perc_cover_beta)) 
testDispersion(Perc_cover_beta) 

####### Model significance
summary(Perc_cover_beta) 
car::Anova(Perc_cover_beta, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)




###    (b) Floral area ---- 

####### negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Floral_area_nbinom <- glmmTMB (Floral_area_rounded ~ Strip_type + (1|Field), 
                                family= nbinom2(),  
                                data = Vegetation)

hist(residuals(Floral_area_nbinom))
plot(simulateResiduals(Floral_area_nbinom)) 
testDispersion(Floral_area_nbinom) 


####### Gamma - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Floral_area_gamma <- glmmTMB (Sum_Floralarea ~ Strip_type + (1|Field), 
                              family= Gamma(link = "log"),  
                              data = Vegetation)

hist(residuals(Floral_area_gamma))
plot(simulateResiduals(Floral_area_gamma)) 
testDispersion(Floral_area_gamma) 


####### Lognormal - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Floral_area_lognormal <- glmmTMB (Sum_Floralarea ~ Strip_type + (1|Field), 
                                  family = lognormal,  
                                  data = Vegetation)

hist(residuals(Floral_area_lognormal))
plot(simulateResiduals(Floral_area_lognormal)) 
testDispersion(Floral_area_lognormal)


####### Model significance
summary(Floral_area_gamma) 
car::Anova(Floral_area_gamma, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)

summary(Floral_area_nbinom) 
car::Anova(Floral_area_nbinom, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)

summary(Floral_area_lognormal) 
car::Anova(Floral_area_lognormal, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)





##    4.2. Pollinators ---- 

### (a) Hoverflies ---- 
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Hoverfly_poisson <- glmmTMB(Hoverfly ~ Strip_type + (1|Field),
                            family = poisson(),
                            data = Pollinators)

hist(residuals(Hoverfly_poisson))
plot(simulateResiduals(Hoverfly_poisson)) 
testDispersion(Hoverfly_poisson) 


####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Hoverfly_nbinom <- glmmTMB(Hoverfly ~ Strip_type + (1|Field),
                           family = nbinom2(),
                           data = Pollinators)

hist(residuals(Hoverfly_nbinom))
plot(simulateResiduals(Hoverfly_nbinom)) 
testDispersion(Hoverfly_nbinom) 


####### Model significance
summary(Hoverfly_nbinom) 
car::Anova(Hoverfly_nbinom, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)





### (b) Honey bees ---- 

####### Poisson, Negative binomial, Binomial or Betabinomial don't work 

####### Transform response variable into binary format
Pollinators$Honey_bee_binary <- ifelse(Pollinators$Honey_bee > 0, 1, 0)

####### Firth Logistic Regression
Honey_bee_wilcox <- wilcox.test(Honey_bee_binary ~ Strip_type, data = Pollinators, exact = FALSE)




### (c) Bumblebees ---- 
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Bumblebee_poisson <- glmmTMB(Bumblebee ~ Strip_type + (1|Field),
                             family = poisson(),
                             data = Pollinators)

hist(residuals(Bumblebee_poisson))
plot(simulateResiduals(Bumblebee_poisson)) 
testDispersion(Bumblebee_poisson) 


####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Bumblebee_nbinom <- glmmTMB(Bumblebee ~ Strip_type + (1|Field),
                            family = nbinom2(),
                            data = Pollinators)

hist(residuals(Bumblebee_nbinom))
plot(simulateResiduals(Bumblebee_nbinom)) 
testDispersion(Bumblebee_nbinom) 


####### Model significance
summary(Bumblebee_nbinom) 
car::Anova(Bumblebee_nbinom, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)





### (d) Solitary bees ---- 

####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Solitary_bee_poisson <- glmmTMB(Solitary_bee ~ Strip_type + (1|Field),
                                family = poisson(),
                                data = Pollinators)

hist(residuals(Solitary_bee_poisson))
plot(simulateResiduals(Solitary_bee_poisson)) 
testDispersion(Solitary_bee_poisson) 


####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Solitary_bee_nbinom <- glmmTMB(Solitary_bee ~ Strip_type + (1|Field),
                               family = nbinom2(),
                               data = Pollinators)

hist(residuals(Solitary_bee_nbinom))
plot(simulateResiduals(Solitary_bee_nbinom)) 
testDispersion(Solitary_bee_nbinom) 


####### Model significance
summary(Solitary_bee_nbinom) 
car::Anova(Solitary_bee_nbinom, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)






# 5. Models - Margin & adjacent field habitats ----

##    5.1. Yellow sticky traps ---- 

###    5.1.1. Natural enemies ---- 

####    (a) Leaf-dwelling spiders ---- 
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Spiders_poisson <- glmmTMB(Spiders ~ Strip_type * Crop * Location + (1|Field),
                                family = poisson(), 
                                data = YST)

hist(residuals(YST_Spiders_poisson))
plot(simulateResiduals(YST_Spiders_poisson)) 
testDispersion(YST_Spiders_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Spiders_negbinom <- glmmTMB(Spiders ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = YST)

hist(residuals(YST_Spiders_negbinom))
plot(simulateResiduals(YST_Spiders_negbinom)) 
testDispersion(YST_Spiders_negbinom) 

####### Model significance
summary(YST_Spiders_negbinom) 
car::Anova(YST_Spiders_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



####    (b) Leaf-dwelling predatory beetles ---- 
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Predatory_beetles_poisson <- glmmTMB(Predatory_beetles ~ Strip_type * Crop * Location + (1|Field),
                               family = poisson(), 
                               data = YST)

hist(residuals(YST_Predatory_beetles_poisson))
plot(simulateResiduals(YST_Predatory_beetles_poisson)) 
testDispersion(YST_Predatory_beetles_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Predatory_beetles_negbinom <- glmmTMB(Predatory_beetles ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = YST)

hist(residuals(YST_Predatory_beetles_negbinom))
plot(simulateResiduals(YST_Predatory_beetles_negbinom)) 
testDispersion(YST_Predatory_beetles_negbinom) 

####### Model significance
summary(YST_Predatory_beetles_negbinom) 
car::Anova(YST_Predatory_beetles_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



####    (c) Leaf-dwelling predatory heteroptera ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Predatory_heteroptera_poisson <- glmmTMB(Predatory_heteroptera ~ Strip_type * Crop * Location + (1|Field),
                               family = poisson(), 
                               data = YST)

hist(residuals(YST_Predatory_heteroptera_poisson))
plot(simulateResiduals(YST_Predatory_heteroptera_poisson)) 
testDispersion(YST_Predatory_heteroptera_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Predatory_heteroptera_negbinom <- glmmTMB(Predatory_heteroptera ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = YST)

hist(residuals(YST_Predatory_heteroptera_negbinom))
plot(simulateResiduals(YST_Predatory_heteroptera_negbinom)) 
testDispersion(YST_Predatory_heteroptera_negbinom) 

####### Model significance
summary(YST_Predatory_heteroptera_negbinom) 
car::Anova(YST_Predatory_heteroptera_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



####    (d) Neuroptera ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Neuroptera_poisson <- glmmTMB(Neuroptera ~ Strip_type * Crop * Location + (1|Field),
                               family = poisson(), 
                               data = YST)

hist(residuals(YST_Neuroptera_poisson))
plot(simulateResiduals(YST_Neuroptera_poisson)) 
testDispersion(YST_Neuroptera_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Neuroptera_negbinom <- glmmTMB(Neuroptera ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = YST)

hist(residuals(YST_Neuroptera_negbinom))
plot(simulateResiduals(YST_Neuroptera_negbinom)) 
testDispersion(YST_Neuroptera_negbinom) 

####### Model significance
summary(YST_Neuroptera_negbinom) 
car::Anova(YST_Neuroptera_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



####    (e) Predatory thrips ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Predatory_thrips_poisson <- glmmTMB(Predatory_thrips ~ Strip_type * Crop * Location + (1|Field),
                               family = poisson(), 
                               data = YST)

hist(residuals(YST_Predatory_thrips_poisson))
plot(simulateResiduals(YST_Predatory_thrips_poisson)) 
testDispersion(YST_Predatory_thrips_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Predatory_thrips_negbinom <- glmmTMB(Predatory_thrips ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = YST)

hist(residuals(YST_Predatory_thrips_negbinom))
plot(simulateResiduals(YST_Predatory_thrips_negbinom)) 
testDispersion(YST_Predatory_thrips_negbinom) 

####### Model significance
summary(YST_Predatory_thrips_negbinom) 
car::Anova(YST_Predatory_thrips_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



####    (f) Syrphidae ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Syrphidae_poisson <- glmmTMB(Syrphidae ~ Strip_type * Crop * Location + (1|Field),
                               family = poisson(), 
                               data = YST)

hist(residuals(YST_Syrphidae_poisson))
plot(simulateResiduals(YST_Syrphidae_poisson)) 
testDispersion(YST_Syrphidae_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Syrphidae_negbinom <- glmmTMB(Syrphidae ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = YST)

hist(residuals(YST_Syrphidae_negbinom))
plot(simulateResiduals(YST_Syrphidae_negbinom)) 
testDispersion(YST_Syrphidae_negbinom) 

####### Model significance
summary(YST_Syrphidae_negbinom) 
car::Anova(YST_Syrphidae_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



####    (g) Parasitic wasps ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Parasitic_wasps_poisson <- glmmTMB(Parasitic_wasps ~ Strip_type * Crop * Location + (1|Field),
                               family = poisson(), 
                               data = YST)

hist(residuals(YST_Parasitic_wasps_poisson))
plot(simulateResiduals(YST_Parasitic_wasps_poisson)) 
testDispersion(YST_Parasitic_wasps_poisson) 


####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Parasitic_wasps_negbinom <- glmmTMB(Parasitic_wasps ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = YST)

hist(residuals(YST_Parasitic_wasps_negbinom))
plot(simulateResiduals(YST_Parasitic_wasps_negbinom)) 
testDispersion(YST_Parasitic_wasps_negbinom) 

####### Model significance
summary(YST_Parasitic_wasps_negbinom) 
car::Anova(YST_Parasitic_wasps_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)




###    5.1.2. Herbivores ---- 

####    (i) Aphids ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Aphids_poisson <- glmmTMB(Aphids ~ Strip_type * Crop * Location + (1|Field),
                              family = poisson(), 
                              data = YST)

hist(residuals(YST_Aphids_poisson))
plot(simulateResiduals(YST_Aphids_poisson)) 
testDispersion(YST_Aphids_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Aphids_negbinom <- glmmTMB(Aphids ~ Strip_type * Crop * Location + (1|Field),
                               family = nbinom2(), 
                               data = YST)

hist(residuals(YST_Aphids_negbinom))
plot(simulateResiduals(YST_Aphids_negbinom)) 
testDispersion(YST_Aphids_negbinom) 

####### Model significance
summary(YST_Aphids_negbinom) 
car::Anova(YST_Aphids_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)




####    (j) Cicadellidae ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Cicadellidae_poisson <- glmmTMB(Cicadellidae ~ Strip_type * Crop * Location + (1|Field),
                              family = poisson(), 
                              data = YST)

hist(residuals(YST_Cicadellidae_poisson))
plot(simulateResiduals(YST_Cicadellidae_poisson)) 
testDispersion(YST_Cicadellidae_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Cicadellidae_negbinom <- glmmTMB(Cicadellidae ~ Strip_type * Crop * Location + (1|Field),
                               family = nbinom2(), 
                               data = YST)

hist(residuals(YST_Cicadellidae_negbinom))
plot(simulateResiduals(YST_Cicadellidae_negbinom)) 
testDispersion(YST_Cicadellidae_negbinom) 

####### Model significance
summary(YST_Cicadellidae_negbinom) 
car::Anova(YST_Cicadellidae_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)





####    (k) Lygaeidae ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Lygaeidae_poisson <- glmmTMB(Lygaeidae ~ Strip_type * Crop * Location + (1|Field),
                              family = poisson(), 
                              data = YST)

hist(residuals(YST_Lygaeidae_poisson))
plot(simulateResiduals(YST_Lygaeidae_poisson)) 
testDispersion(YST_Lygaeidae_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Lygaeidae_negbinom <- glmmTMB(Lygaeidae ~ Strip_type * Crop * Location + (1|Field),
                               family = nbinom2(), 
                               data = YST)

hist(residuals(YST_Lygaeidae_negbinom))
plot(simulateResiduals(YST_Lygaeidae_negbinom)) 
testDispersion(YST_Lygaeidae_negbinom) 

####### Model significance
summary(YST_Lygaeidae_negbinom) 
car::Anova(YST_Lygaeidae_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)


####    (l) Sap beetles ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Nitidulidae_poisson <- glmmTMB(Nitidulidae ~ Strip_type * Crop * Location + (1|Field),
                              family = poisson(), 
                              data = YST)

hist(residuals(YST_Nitidulidae_poisson))
plot(simulateResiduals(YST_Nitidulidae_poisson)) 
testDispersion(YST_Nitidulidae_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Nitidulidae_negbinom <- glmmTMB(Nitidulidae ~ Strip_type * Crop * Location + (1|Field),
                               family = nbinom2(), 
                               data = YST)

hist(residuals(YST_Nitidulidae_negbinom))
plot(simulateResiduals(YST_Nitidulidae_negbinom)) 
testDispersion(YST_Nitidulidae_negbinom) 

####### Model significance
summary(YST_Nitidulidae_negbinom) 
car::Anova(YST_Nitidulidae_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)




####    (m) Psylloidea ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Psylloidea_poisson <- glmmTMB(Psylloidea ~ Strip_type * Crop * Location + (1|Field),
                              family = poisson(), 
                              data = YST)

hist(residuals(YST_Psylloidea_poisson))
plot(simulateResiduals(YST_Psylloidea_poisson)) 
testDispersion(YST_Psylloidea_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Psylloidea_negbinom <- glmmTMB(Psylloidea ~ Strip_type * Crop * Location + (1|Field),
                               family = nbinom2(), 
                               data = YST)

hist(residuals(YST_Psylloidea_negbinom))
plot(simulateResiduals(YST_Psylloidea_negbinom)) 
testDispersion(YST_Psylloidea_negbinom) 

####### Model significance
summary(YST_Psylloidea_negbinom) 
car::Anova(YST_Psylloidea_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



####    (n) Thysanoptera ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Thysanoptera_poisson <- glmmTMB(Thysanoptera ~ Strip_type * Crop * Location + (1|Field),
                               family = poisson(), 
                               data = YST)

hist(residuals(YST_Thysanoptera_poisson))
plot(simulateResiduals(YST_Thysanoptera_poisson)) 
testDispersion(YST_Thysanoptera_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
YST_Thysanoptera_negbinom <- glmmTMB(Thysanoptera ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = YST)

hist(residuals(YST_Thysanoptera_negbinom))
plot(simulateResiduals(YST_Thysanoptera_negbinom)) 
testDispersion(YST_Thysanoptera_negbinom) 

####### Model significance
summary(YST_Thysanoptera_negbinom) 
car::Anova(YST_Thysanoptera_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)






##    5.2. Pitfall traps ---- 

###    (a) Ground-dwelling spiders ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Pitfall_Spiders_poisson <- glmmTMB(Ground_Spiders ~ Strip_type * Crop * Location + (1|Field),
                               family = poisson(), 
                               data = Pitfall)

hist(residuals(Pitfall_Spiders_poisson))
plot(simulateResiduals(Pitfall_Spiders_poisson)) 
testDispersion(Pitfall_Spiders_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Pitfall_Spiders_negbinom <- glmmTMB(Ground_Spiders ~ Strip_type * Crop * Location + (1|Field),
                                family = nbinom2(), 
                                data = Pitfall)

hist(residuals(Pitfall_Spiders_negbinom))
plot(simulateResiduals(Pitfall_Spiders_negbinom)) 
testDispersion(Pitfall_Spiders_negbinom) 

####### Model significance
summary(Pitfall_Spiders_negbinom) 
car::Anova(Pitfall_Spiders_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)




###    (b) Carabid beeltes ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Pitfall_Carabidae_poisson <- glmmTMB(Carabidae ~ Strip_type * Crop * Location + (1|Field),
                                   family = poisson(), 
                                   data = Pitfall)

hist(residuals(Pitfall_Carabidae_poisson))
plot(simulateResiduals(Pitfall_Carabidae_poisson)) 
testDispersion(Pitfall_Carabidae_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Pitfall_Carabidae_negbinom <- glmmTMB(Carabidae ~ Strip_type * Crop * Location + (1|Field),
                                    family = nbinom2(), 
                                    data = Pitfall)

hist(residuals(Pitfall_Carabidae_negbinom))
plot(simulateResiduals(Pitfall_Carabidae_negbinom)) 
testDispersion(Pitfall_Carabidae_negbinom) 

####### Model significance
summary(Pitfall_Carabidae_negbinom) 
car::Anova(Pitfall_Carabidae_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)




###    (c) Rove beetles ----
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Pitfall_Staphylinidae_poisson <- glmmTMB(Staphylinidae ~ Strip_type * Crop * Location + (1|Field),
                                   family = poisson(), 
                                   data = Pitfall)

hist(residuals(Pitfall_Staphylinidae_poisson))
plot(simulateResiduals(Pitfall_Staphylinidae_poisson)) 
testDispersion(Pitfall_Staphylinidae_poisson) 

####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
Pitfall_Staphylinidae_negbinom <- glmmTMB(Staphylinidae ~ Strip_type * Crop * Location + (1|Field),
                                    family = nbinom2(), 
                                    data = Pitfall)

hist(residuals(Pitfall_Staphylinidae_negbinom))
plot(simulateResiduals(Pitfall_Staphylinidae_negbinom)) 
testDispersion(Pitfall_Staphylinidae_negbinom) 

####### Model significance
summary(Pitfall_Staphylinidae_negbinom) 
car::Anova(Pitfall_Staphylinidae_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)





##    5.3. Bait lamina strips ---- 

###    (a) Decomposition rate ----
# Check distribution for this one and all the beta or betabionomial so it is the same in all
# Do figure but put it in a separate figure
# Extract data for table also


####### Beta family - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
BL_beta <- glmmTMB(Decomposition_rate ~ Strip_type * Crop * Location + (1|Field), 
                   family = beta_family(), 
                   data = BLS)

hist(residuals(BL_beta))
plot(simulateResiduals(BL_beta)) 
testDispersion(BL_beta) 

####### binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
BL_binomial <- glmmTMB (Decomposition_rate ~ Strip_type * Crop * Location + (1|Field), 
                        family=betabinomial(link="logit"), 
                        weights= Total_number_holes,
                        data = BLS)


hist(residuals(BL_binomial))
plot(simulateResiduals(BL_binomial)) 
testDispersion(BL_binomial) 

####### Model significance
summary(BL_binomial) 
car::Anova(BL_binomial, type = 3) # to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)





# 6. Models - Adjacent field habitats ----

##    6.1. Tiller counts ---- 

###    (a) Natural enemies ---- 
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
TC_NE_poisson <- glmer (NE ~ Strip_type * Crop + (1|Field), 
                        family=poisson(),  
                        data = TC)

hist(residuals(TC_NE_poisson))
plot(simulateResiduals(TC_NE_poisson)) 
testDispersion(TC_NE_poisson) 


####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
TC_NE_negbinom <- glmmTMB (NE ~ Strip_type * Crop + (1|Field), 
                           family=nbinom2(),  
                           data = TC)

hist(residuals(TC_NE_negbinom))
plot(simulateResiduals(TC_NE_negbinom)) 
testDispersion(TC_NE_negbinom) 

####### Model significance
summary(TC_NE_negbinom) 
car::Anova(TC_NE_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



###    (b) Thrips ---- 
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
TC_Thrips_poisson <- glmer (Thrips_total ~ Strip_type * Crop + (1|Field), 
                            family=poisson(),  
                            data = TC)

hist(residuals(TC_Thrips_poisson))
plot(simulateResiduals(TC_Thrips_poisson)) 
testDispersion(TC_Thrips_poisson) 


####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
TC_Thrips_negbinom <- glmmTMB (Thrips_total ~ Strip_type * Crop + (1|Field), 
                               family=nbinom2(),  
                               data = TC)

hist(residuals(TC_Thrips_negbinom))
plot(simulateResiduals(TC_Thrips_negbinom)) 
testDispersion(TC_Thrips_negbinom) 

####### Model significance
summary(TC_Thrips_negbinom) 
car::Anova(TC_Thrips_negbinom, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)



###    (c) Aphids ---- 
####### Poisson - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
TC_Aphids_poisson <- glmer (Aphids ~ Strip_type * Crop + (1|Field), 
                            family=poisson(),  
                            data = TC)

hist(residuals(TC_Aphids_poisson))
plot(simulateResiduals(TC_Aphids_poisson)) 
testDispersion(TC_Aphids_poisson) 


####### Negative binomial - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
TC_Aphids_negbinom <- glmmTMB (Aphids ~ Strip_type * Crop + (1|Field), 
                               family=nbinom2(),  
                               data = TC)

hist(residuals(TC_Aphids_negbinom))
plot(simulateResiduals(TC_Aphids_negbinom)) 
testDispersion(TC_Aphids_negbinom) 

####### Model significance
summary(TC_Aphids_poisson) 
car::Anova(TC_Aphids_poisson, type=3) #ANOVA to get significance of the main factor, add type 3 if we want to investigate the interaction, otherwise type 2)





##    6.2. Aphid cards ---- 

###    (a) Predation rates - leaf ---- 
####### Beta family - Build model & validate
# Used shringking formula form Smithson & Verkuilen 2006
options(contrasts = c("contr.sum", "contr.poly"))
AC_leaf_beta <- glmmTMB(Predation_rate_leaf ~ Strip_type * Crop + (1|Field), 
                        family = beta_family(link = "logit"), 
                        data = AC_leaf)

hist(residuals(AC_leaf_beta))
plot(simulateResiduals(AC_leaf_beta)) 
testDispersion(AC_leaf_beta) 


####### Betabinomial - Build model & validate
# Used shringking formula form Smithson & Verkuilen 2006
options(contrasts = c("contr.sum", "contr.poly"))
AC_leaf_betabinomial <- glmmTMB(Predation_rate_leaf ~ Strip_type * Crop + (1|Field), 
                                family=betabinomial(link="logit"), 
                                weights= Aphids_set_up,
                                data = AC_leaf)

hist(residuals(AC_leaf_betabinomial))
plot(simulateResiduals(AC_leaf_betabinomial)) 
testDispersion(AC_leaf_betabinomial) 



####### Model significance
summary(AC_leaf_betabinomial) 
car::Anova(AC_leaf_betabinomial, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)



###    (b) Predation rates - ground ---- 
####### Beta family - Build model & validate
options(contrasts = c("contr.sum", "contr.poly"))
AC_ground_beta <- glmmTMB(Predation_rate_ground ~ Strip_type * Crop + (1|Field), 
                          family = beta_family(link = "logit"), 
                          data = AC_ground)

hist(residuals(AC_ground_beta))
plot(simulateResiduals(AC_ground_beta)) 
testDispersion(AC_ground_beta) 



####### Betabinomial - Build model & validate
# Used shringking formula form Smithson & Verkuilen 2006
options(contrasts = c("contr.sum", "contr.poly"))
AC_ground_betabinomial <- glmmTMB(Predation_rate_ground ~ Strip_type * Crop+ (1|Field), 
                                  family=betabinomial(link="logit"), 
                                  weights= Aphids_set_up,
                                  data = AC_ground)

hist(residuals(AC_ground_betabinomial))
plot(simulateResiduals(AC_ground_betabinomial)) 
testDispersion(AC_ground_betabinomial) 



####### Model significance
summary(AC_ground_betabinomial) 
car::Anova(AC_ground_betabinomial, type=3) #ANOVA, add type 3 to investigate interaction, otherwise type 2)











# 7. Post-hoc tests ----

##    7.1. Margin habitats ---- 

####### Model list
#Honey bees not included because it is a wilcox test result
Margin_modelslist <- list(
  Perc_cover_beta,
  Floral_area_nbinom,
  Hoverfly_nbinom,
  Bumblebee_nbinom,
  Solitary_bee_nbinom
)

####### Model names
Margin_modelsnames <- c(
  "Proportion vegetation cover",
  "Total floral area (cm2/m2)", 
  "Hoverflies", 
  "Bumblebees", 
  "Solitary bees"
)


####### Function to get emmeans and CI
# Create function to extract and transform EMMs
extract_emm_and_pairs_Margin <- function(models, names, factor = ~ Strip_type) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  # Combine and return
  return(list(
    Margin_EMMs = bind_rows(all_emms),
    Margin_Pairs = bind_rows(all_pairs)
  ))
}
  


# Run function 
margin_results <- extract_emm_and_pairs_Margin(Margin_modelslist, Margin_modelsnames)

# View results
print(margin_results$Margin_EMMs)
print(margin_results$Margin_Pairs)

# Save tables
write.xlsx(margin_results$Margin_EMMs, "Table_Margin_EMMs.xlsx")
write.xlsx(margin_results$Margin_Pairs, "Table_Margin_Pairs.xlsx")








##    7.2. Margin & adjacent field habitats ---- 

####### Model list
MarginField_modelslist <- list(
  YST_Spiders_negbinom,
  YST_Predatory_beetles_negbinom,
  YST_Predatory_heteroptera_negbinom,
  YST_Neuroptera_negbinom,
  YST_Predatory_thrips_negbinom,
  YST_Syrphidae_negbinom,
  YST_Parasitic_wasps_negbinom,
  YST_Aphids_negbinom,
  YST_Cicadellidae_negbinom,
  YST_Lygaeidae_negbinom,
  YST_Nitidulidae_negbinom,
  YST_Psylloidea_negbinom,
  YST_Thysanoptera_negbinom,
  Pitfall_Spiders_negbinom,
  Pitfall_Carabidae_negbinom,
  Pitfall_Staphylinidae_negbinom,
  BL_binomial
)

####### Model names
MarginField_modelsnames <- c(
  "Leaf-dwelling spiders",
  "Leaf-dwelling predatory beetles", 
  "Predatory bugs", 
  "Lacewings", 
  "Predatory thrips", 
  "Hoverflies",
  "Parasitic wasps",
  "Aphids",
  "Leafhoppers",
  "Seed bugs",
  "Sap beetles", 
  "Psyllids",
  "Thrips",
  "Ground-dwelling spiders",
  "Carabid beetles",
  "Rove beetles",
  "Decompostion rate"
)





### 7.2.1.Posthoc for treatment ####### 

# Create function to extract and transform EMMs
extract_emm_and_pairs_FieldMargin_treatment <- function(models, names, factor = ~ Strip_type) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  # Combine and return
  return(list(
    FieldMargin_treatment_EMMs = bind_rows(all_emms),
    FieldMargin_treatment_Pairs = bind_rows(all_pairs)
  ))
}



# Run function 
FieldMargin_treatment_results <- extract_emm_and_pairs_FieldMargin_treatment(MarginField_modelslist, MarginField_modelsnames)

# View results
print(FieldMargin_treatment_results$FieldMargin_treatment_EMMs)
print(FieldMargin_treatment_results$FieldMargin_treatment_Pairs)


# Save to Excel
write.xlsx(FieldMargin_treatment_results$FieldMargin_treatment_EMMs, "Table_FieldMargin_treatment_EMMs.xlsx")
write.xlsx(FieldMargin_treatment_results$FieldMargin_treatment_Pairs, "Table_FieldMargin_treatment_Pairs.xlsx")





### 7.2.2.Posthoc for location ####### 
# Create function to extract and transform EMMs
extract_emm_and_pairs_FieldMargin_location <- function(models, names, factor = ~ Location) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  # Combine and return
  return(list(
    FieldMargin_location_EMMs = bind_rows(all_emms),
    FieldMargin_location_Pairs = bind_rows(all_pairs)
  ))
}



# Run function 
FieldMargin_location_results <- extract_emm_and_pairs_FieldMargin_location(MarginField_modelslist, MarginField_modelsnames)

# View results
print(FieldMargin_location_results$FieldMargin_location_EMMs)
print(FieldMargin_location_results$FieldMargin_location_Pairs)


# Save to Excel
write.xlsx(FieldMargin_location_results$FieldMargin_location_EMMs, "Table_FieldMargin_location_EMMs.xlsx")
write.xlsx(FieldMargin_location_results$FieldMargin_location_Pairs, "Table_FieldMargin_location_Pairs.xlsx")









### 7.2.3.Posthoc for the interaction between treatment and location ####### 
####### Function to extract and transform EMMs and pairwise comparisons
extract_emm_and_pairs_FieldMargin_interaction <- function(models, names, factor = ~ Strip_type * Location) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  return(list(
    FieldMargin_interaction_EMMs = bind_rows(all_emms),
    FieldMargin_interaction_Pairs = bind_rows(all_pairs)
  ))
}

# Run function 
fieldmargin_results_interaction <- extract_emm_and_pairs_FieldMargin_interaction(MarginField_modelslist, MarginField_modelsnames)

# View results
print(fieldmargin_results_interaction$FieldMargin_interaction_EMMs)
print(fieldmargin_results_interaction$FieldMargin_interaction_Pairs)

# Save to Excel
write.xlsx(fieldmargin_results_interaction$FieldMargin_interaction_EMMs, "Table_FieldMargin_interaction_EMMs.xlsx")
write.xlsx(fieldmargin_results_interaction$FieldMargin_interaction_Pairs, "Table_FieldMargin_interaction_Pairs.xlsx")










##    7.3. Field habitats ---- 
####### Model list
Field_modelslist <- list(
  TC_NE_negbinom,
  TC_Thrips_negbinom,
  TC_Aphids_poisson,
  AC_leaf_betabinomial,
  AC_ground_betabinomial
)

####### Model names
Field_modelsnames <- c(
  "Natural enemies",
  "Thrips", 
  "Aphids", 
  "Aphid predation rates by leaf-dwellers", 
  "Aphid predation rates by ground-dwellers"
)

####### Function to get emmeans and CI
extract_emm_and_pairs_Field <- function(models, names, factor = ~ Strip_type) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  # Combine and return
  return(list(
    Field_EMMs = bind_rows(all_emms),
    Field_Pairs = bind_rows(all_pairs)
  ))
}

# Run function 
field_results <- extract_emm_and_pairs_Field(Field_modelslist, Field_modelsnames)

# View results
print(field_results$Field_EMMs)
print(field_results$Field_Pairs)

# Save tables
write.xlsx(field_results$Field_EMMs, "Table_Field_EMMs.xlsx")
write.xlsx(field_results$Field_Pairs, "Table_Field_Pairs.xlsx")





# 8. Tables ----

##    8.1. Margin habitats ---- 

### Table 1. Accumulated counts and model results. Margin characterization & pollinator observations ##########

# Precompute Honey bee counts and model estimates
##### As honey bees are done with wilcox text, this preliminary step is necessary, later will be added into the function with the other models #
honeybee_modelsummary <- data.frame(
  Model = "Honey bees",
  Test_statistics = Honey_bee_wilcox$statistic,
  p = Honey_bee_wilcox$p.value,
  Test = "Wilcoxon",
  Significance = case_when(
    Honey_bee_wilcox$p.value < 0.001 ~ "***",
    Honey_bee_wilcox$p.value < 0.01  ~ "**",
    Honey_bee_wilcox$p.value < 0.05  ~ "*",
    Honey_bee_wilcox$p.value < 0.1   ~ "(*)",
    TRUE                             ~ ""
  ),
  Sum_Flower_strip = sum(Pollinators$Honey_bee[Pollinators$Strip_type == "Flower_strip"], na.rm = TRUE),
  Sum_Control = sum(Pollinators$Honey_bee[Pollinators$Strip_type == "Control"], na.rm = TRUE)
)

print(honeybee_modelsummary)




####### Model list
Margin_modelslist <- list(
  Perc_cover_beta,
  Floral_area_nbinom,
  Hoverfly_nbinom,
  Honey_bee_wilcox,
  Bumblebee_nbinom,
  Solitary_bee_nbinom
)

####### Model names
Margin_modelsnames <- c(
  "Proportion vegetation cover",
  "Total floral area (cm2/m2)", 
  "Hoverflies", 
  "Honey bees", 
  "Bumblebees", 
  "Solitary bees"
)

####### Response variables
Margin_response_variables <- c(
  "Perc_cover_adj",         
  "Floral_area_rounded",        
  "Hoverfly",         
  "Honey_bee",         
  "Bumblebee",         
  "Solitary_bee"       
)

####### Loop through models to perform sums and ANOVAs
Margin_modelsummary <- lapply(seq_along(Margin_modelslist), function(i) {
  model <- Margin_modelslist[[i]]
  model_name <- Margin_modelsnames[i]
  response_var <- Margin_response_variables[i]
  
  if (model_name == "Honey bees") {
    # Inject precomputed Honey bees parameters
    return(honeybee_modelsummary)
  } else {
    # Analysis for the other models (e.g., glmmTMB)
    dataset <- model.frame(model)
    
    # Calculate sums
    strip_type_sums <- dataset %>%
      group_by(Strip_type) %>%
      summarise(
        Value = if (response_var == "Perc_cover_adj") {
          mean(.data[[response_var]], na.rm = TRUE)  # Use mean for vegetation cover
        } else {
          round(sum(.data[[response_var]], na.rm = TRUE))  # Sum for others
        },
        .groups = "drop"
      ) %>%
      pivot_wider(names_from = Strip_type, values_from = Value, names_prefix = "Value_")
    
    # Perform ANOVA (Type III)
    anova_result <- Anova(model, type = 3)
    
    # Extract results
    model_results <- data.frame(
      Model = model_name,
      Term = rownames(anova_result),
      Test_statistics = as.numeric(anova_result$`Chisq`),  # Use Chi-squared value as test statistic
      p = as.numeric(anova_result$`Pr(>Chisq)`),
      Test = case_when(
        model_name == "Proportion vegetation cover" ~ "Beta",
        TRUE ~ "Negative binomial"
      ),# Specify the model method
      Significance = case_when(
        as.numeric(anova_result$`Pr(>Chisq)`) < 0.001 ~ "***",
        as.numeric(anova_result$`Pr(>Chisq)`) < 0.01  ~ "**",
        as.numeric(anova_result$`Pr(>Chisq)`) < 0.05  ~ "*",
        as.numeric(anova_result$`Pr(>Chisq)`) < 0.1   ~ "(*)",
        TRUE                                         ~ ""
      )
    )
    
    # Combine results
    combined_results <- model_results %>%
      filter(Term == "Strip_type") %>%
      select(-Term) %>%
      bind_cols(strip_type_sums)
    
    return(combined_results)
  }
}) %>% bind_rows()

print(Margin_modelsummary)


####### Reorder and format the final table
Margin_modelsummary <- Margin_modelsummary %>%
  mutate(
    `Value_Flower_strip` = as.numeric(`Value_Flower_strip`),  # Ensure numeric
    `Value_Control` = as.numeric(`Value_Control`),  # Ensure numeric
    `Value_Flower_strip` = ifelse(
      Model == "Proportion vegetation cover",
      formatC(`Value_Flower_strip`, format = "f", digits = 6),  # 6 decimals for vegetation cover
      round(`Value_Flower_strip`)  # No decimals for other rows
    ),
    `Value_Control` = ifelse(
      Model == "Proportion vegetation cover",
      formatC(`Value_Control`, format = "f", digits = 6),  # 6 decimals for vegetation cover
      round(`Value_Control`)  # No decimals for other rows
    ),
    `Test_statistics` = round(as.numeric(`Test_statistics`), 2),  # Round to 2 decimals
    p = formatC(as.numeric(p), format = "f", digits = 3)  # Round p-value to 3 decimals
  )%>%
  transmute(
    Model = Model,
    `Flower strip (N)` = Value_Flower_strip,  # Rename column
    `Control (N)` = Value_Control,  # Rename column
    Test = Test,
    `Test statistics (W-Statistic or Chi-squared)` = Test_statistics,  # Rename column
    p = p,
    Significance = Significance
  )

print(Margin_modelsummary)

#######  Save table
write.xlsx(Margin_modelsummary, 'Table1.xlsx')









##    8.2. Margin & adjacent field habitats ---- 

### Table 2. Accumulated counts and model results. YST, pitfall traps & decomposition ##########

####### Model list
MarginField_modelslist <- list(
  YST_Spiders_negbinom,
  YST_Predatory_beetles_negbinom,
  YST_Predatory_heteroptera_negbinom,
  YST_Neuroptera_negbinom,
  YST_Predatory_thrips_negbinom,
  YST_Syrphidae_negbinom,
  YST_Parasitic_wasps_negbinom,
  YST_Aphids_negbinom,
  YST_Cicadellidae_negbinom,
  YST_Lygaeidae_negbinom,
  YST_Nitidulidae_negbinom,
  YST_Psylloidea_negbinom,
  YST_Thysanoptera_negbinom,
  Pitfall_Spiders_negbinom,
  Pitfall_Carabidae_negbinom,
  Pitfall_Staphylinidae_negbinom,
  BL_binomial
)

####### Model names
MarginField_modelsnames <- c(
  "Leaf-dwelling spiders",
  "Leaf-dwelling predatory beetles", 
  "Predatory bugs", 
  "Lacewings", 
  "Predatory thrips", 
  "Hoverflies",
  "Parasitic wasps",
  "Aphids",
  "Leafhoppers",
  "Seed bugs",
  "Sap beetles", 
  "Psyllids",
  "Thrips",
  "Ground-dwelling spiders",
  "Carabid beetles",
  "Rove beetles",
  "Decompostion rate"
)


####### Loop through models and generate the model output table
MarginField_modelsummary <- lapply(seq_along(MarginField_modelslist), function(i) {
  model <- MarginField_modelslist[[i]]
  model_name <- MarginField_modelsnames[i]
  
  # Dynamically assign residual distribution
  residual_distribution <- if (model_name == "Decompostion rate") {
    "Binomial"
  } else {
    "Negative binomial"
  }
  
  # Extract dataset and response variable
  dataset <- model.frame(model)
  response_var <- as.character(attr(terms(model), "variables"))[2]  # Extract response variable name
  
  # Calculate sums nested by Location and Strip_type
  nested_sums <- dataset %>%
    group_by(Location, Strip_type) %>%
    summarise(Sum = sum(.data[[response_var]], na.rm = TRUE), .groups = "drop") %>%
    unite("Location_StripType", Location, Strip_type) %>%  # Combine Location and Strip_type into a single column
    pivot_wider(names_from = Location_StripType, values_from = Sum, names_prefix = "Sum_")
  
  # Perform ANOVA (Type III)
  anova_result <- Anova(model, type = 3)
  
  # Format ANOVA results
  anova_summary <- data.frame(
    Term = rownames(anova_result),
    Chisq = round(as.numeric(anova_result$`Chisq`), 2),
    p_value = round(as.numeric(anova_result$`Pr(>Chisq)`), 3),
    Significance = case_when(
      as.numeric(anova_result$`Pr(>Chisq)`) < 0.001 ~ "***",
      as.numeric(anova_result$`Pr(>Chisq)`) < 0.01  ~ "**",
      as.numeric(anova_result$`Pr(>Chisq)`) < 0.05  ~ "*",
      as.numeric(anova_result$`Pr(>Chisq)`) < 0.1   ~ "(*)",
      TRUE                                         ~ ""
    )
  ) %>%
    pivot_wider(
      names_from = Term,
      values_from = c(Chisq, p_value, Significance),
      names_sep = "_",
      values_fill = list(Chisq = NA, p_value = NA, Significance = NA)
    )
  
  # Combine all results into a single row
  model_row <- data.frame(
    Model = model_name, 
    Residual_Distribution = residual_distribution
  ) %>%
    bind_cols(nested_sums, anova_summary)  # Combine all results
  
  return(model_row)
}) %>% bind_rows()


####### Define the desired column order
MarginField_desired_column_order <- c(
  "Model",
  "Sum_Margin_Flower_strip",
  "Sum_Margin_Control",
  "Sum_Crop_Flower_strip",
  "Sum_Crop_Control",
  "Residual_Distribution",
  "Chisq_Strip_type",   
  "p_value_Strip_type",  
  "Significance_Strip_type",
  "Chisq_Crop",   
  "p_value_Crop",  
  "Significance_Crop",
  "Chisq_Location",   
  "p_value_Location",  
  "Significance_Location",
  "Chisq_Strip_type:Crop",   
  "p_value_Strip_type:Crop",  
  "Significance_Strip_type:Crop",
  "Chisq_Strip_type:Location",   
  "p_value_Strip_type:Location",  
  "Significance_Strip_type:Location",
  "Chisq_Crop:Location",   
  "p_value_Crop:Location",  
  "Significance_Crop:Location",  
  "Chisq_Strip_type:Crop:Location",   
  "p_value_Strip_type:Crop:Location",  
  "Significance_Strip_type:Crop:Location"
)


#######  Reorder columns in the data frame
MarginField_modelsummary <- MarginField_modelsummary %>%
  select(all_of(MarginField_desired_column_order))

#######  Check the reordered results
View(MarginField_modelsummary)
print(MarginField_modelsummary)




#######  Save table
write.xlsx(MarginField_modelsummary, 'Table2.xlsx')



##In this table decomposition rate should be substituted by average (A) now the table is showing accumulated count (N) for all the models
Average_decomposition_rate <- BLS %>%
  group_by(Strip_type, Location) %>%   # Group by treatment and location
  summarize(Average_decomposition_rate = mean(Decomposition_rate, na.rm = TRUE))  # Calculate the average
print(Average_decomposition_rate)


##    8.3. Field habitats ---- 

### Table 3. Accumulated counts and model results. Tiller counts & aphid cards ##########

####### Model list
Field_modelslist <- list(
  TC_NE_negbinom ,
  TC_Thrips_negbinom,
  TC_Aphids_poisson,
  AC_leaf_betabinomial,
  AC_ground_betabinomial
)

####### Model names
Field_modelsnames <- c(
  "Natural enemies",
  "Thrips", 
  "Aphids", 
  "Aphid predation rates by leaf-dwellers", 
  "Aphid predation rates by ground-dwellers"
)

####### Response variables
Field_response_variables <- c(
  "NE",        
  "Thrips_total",         
  "Aphids",         
  "Predation_rate_leaf",         
  "Predation_rate_ground"       
)



####### Loop through models and generate the model output table
Field_modelsummary <- lapply(seq_along(Field_modelslist), function(i) {
  model <- Field_modelslist[[i]]
  model_name <- Field_modelsnames[i]
  
  # Dynamically assign residual distribution
  residual_distribution <- if (model_name %in% c("Aphid predation rates by leaf-dwellers", 
                                                 "Aphid predation rates by ground-dwellers")) {
    "Betabinomial"
  } else {
    "Negative binomial"
  }
  
  # Extract dataset and response variable
  dataset <- model.frame(model)
  response_var <- as.character(attr(terms(model), "variables"))[2]  # Extract response variable name
  
  # Calculate sums nested by Strip_type
  strip_type_sums <- dataset %>%
    group_by(Strip_type) %>%
    summarise(Sum = sum(.data[[response_var]], na.rm = TRUE), .groups = "drop") %>%
    pivot_wider(names_from = Strip_type, values_from = Sum, names_prefix = "Sum_")
  
  # Perform ANOVA (Type III)
  anova_result <- Anova(model, type = 3)
  
  # Format ANOVA results
  anova_summary <- data.frame(
    Term = c("Strip_type", "Crop", "Strip_type:Crop"),  # Include all relevant terms
    Chisq = round(as.numeric(anova_result$`Chisq`[match(c("Strip_type", "Crop", "Strip_type:Crop"), rownames(anova_result))]), 2),
    p_value = round(as.numeric(anova_result$`Pr(>Chisq)`[match(c("Strip_type", "Crop", "Strip_type:Crop"), rownames(anova_result))]), 3),
    Significance = case_when(
      as.numeric(anova_result$`Pr(>Chisq)`[match(c("Strip_type", "Crop", "Strip_type:Crop"), rownames(anova_result))]) < 0.001 ~ "***",
      as.numeric(anova_result$`Pr(>Chisq)`[match(c("Strip_type", "Crop", "Strip_type:Crop"), rownames(anova_result))]) < 0.01  ~ "**",
      as.numeric(anova_result$`Pr(>Chisq)`[match(c("Strip_type", "Crop", "Strip_type:Crop"), rownames(anova_result))]) < 0.05  ~ "*",
      as.numeric(anova_result$`Pr(>Chisq)`[match(c("Strip_type", "Crop", "Strip_type:Crop"), rownames(anova_result))]) < 0.1   ~ "(*)",
      TRUE                                                                                                               ~ ""
    )
  ) %>%
    pivot_wider(
      names_from = Term,
      values_from = c(Chisq, p_value, Significance),
      names_sep = "_",
      values_fill = list(Chisq = NA, p_value = NA, Significance = NA)
    )
  
  # Combine all results into a single row
  model_row <- data.frame(
    Model = model_name, 
    Residual_Distribution = residual_distribution
  ) %>%
    bind_cols(strip_type_sums, anova_summary)  # Combine all results
  
  return(model_row)
}) %>% bind_rows()

####### Define the desired column order
Field_desired_column_order <- c(
  "Model",
  "Sum_Flower_strip",
  "Sum_Control",
  "Residual_Distribution",
  "Chisq_Strip_type",   
  "p_value_Strip_type",  
  "Significance_Strip_type",
  "Chisq_Crop",   
  "p_value_Crop",  
  "Significance_Crop",
  "Chisq_Strip_type:Crop",   
  "p_value_Strip_type:Crop",  
  "Significance_Strip_type:Crop"
)

####### Reorder columns in the data frame
Field_modelsummary <- Field_modelsummary %>%
  select(all_of(Field_desired_column_order))

####### Check the reordered results
View(Field_modelsummary)
print(Field_modelsummary)

####### Save the table
write.xlsx(Field_modelsummary, "Table3.xlsx")



##In this table aphid predation rates should be substituted by average (A) now the table is showing accumulated count (N) for all the models
Aphid_predation_rate_leaf <- AC_leaf %>%
  group_by(Strip_type) %>%   # Group by treatment
  summarize(Aphid_predation_rate_leaf = mean(Predation_rate_leaf, na.rm = TRUE))  # Calculate the average
print(Aphid_predation_rate_leaf)

Aphid_predation_rate_ground <- AC_ground %>%
  group_by(Strip_type) %>%   # Group by treatment
  summarize(Aphid_predation_rate_ground = mean(Predation_rate_ground, na.rm = TRUE))  # Calculate the average
print(Aphid_predation_rate_ground)





# 9. Figures ----

##    9.1. Margin habitats ---- 


### Table SX. Preliminary table that gives the model coefficients necessary for Figure 1 ##########

### Precompute Honey bees median difference and confidence intervals (Wilcoxon doesn't directly give them)
# Calculate the observed median difference
median_diff <- with(Pollinators, 
                    median(Honey_bee[Strip_type == "Flower_strip"], na.rm = TRUE) - 
                      median(Honey_bee[Strip_type == "Control"], na.rm = TRUE))

# Define the bootstrap function
bootstrap_median_diff <- function(data, indices) {
  d <- data[indices, ]  # Resample the data
  flower_strip_median <- median(d$Honey_bee[d$Strip_type == "Flower_strip"], na.rm = TRUE)
  control_median <- median(d$Honey_bee[d$Strip_type == "Control"], na.rm = TRUE)
  return(flower_strip_median - control_median)
}

# Perform the bootstrap
set.seed(123)
boot_results <- boot(data = Pollinators, statistic = bootstrap_median_diff, R = 1000)

# Extract the confidence intervals
ci <- boot.ci(boot_results, type = "perc")$percent[4:5]  # 4th and 5th elements are the lower and upper bounds

# Print the results
print(median_diff)  # Median difference
print(ci)           # 95% confidence interval


# Calculate the standard error of the median difference
se_median_diff <- sd(boot_results$t)  # Bootstrap standard deviation of medians

# Calculate the parametric CI (normal approximation)
ci_parametric <- c(median_diff - 1.96 * se_median_diff, median_diff + 1.96 * se_median_diff)

# Log-transform the median difference and its parametric CI
log_median_diff <- log(median_diff)
log_ci_low <- log(ci_parametric[1])
log_ci_high <- log(ci_parametric[2])


# Update honeybee_parameters with the log-transformed values
honeybee_parameters <- data.frame(
  response = "Honey bees",
  Coefficient = log_median_diff,
  CI_low = log_ci_low,
  CI_high = log_ci_high,
  p = Honey_bee_wilcox$p.value,
  significance = case_when(
    Honey_bee_wilcox$p.value < 0.001 ~ "***",
    Honey_bee_wilcox$p.value < 0.01  ~ "**",
    Honey_bee_wilcox$p.value < 0.05  ~ "*",
    Honey_bee_wilcox$p.value < 0.1   ~ "(*)",
    TRUE                             ~ ""
  )
) %>%
  mutate(Parameter = "Strip_type1")

#Visualize
print(honeybee_parameters)







# Define a function to extract model parameters and add significance
extract_model_parameters <- function(model, response_name) {
  as.data.frame(model_parameters(model, type = 'response')) %>%
    mutate(response = response_name,
           significance = case_when(
             p < 0.001 ~ "***",
             p < 0.01  ~ "**",
             p < 0.05  ~ "*",
             p < 0.1   ~ "(*)",
             TRUE      ~ ""
           ))
}



# Apply the function to each model and combine the results for all models
Margin_parameters_combined <- bind_rows(
  extract_model_parameters(Perc_cover_beta, 'Proportion vegetation cover'),
  extract_model_parameters(Floral_area_nbinom, 'Total floral area (cm2/m2)'),
  extract_model_parameters(Hoverfly_nbinom, 'Hoverflies'),
  honeybee_parameters,  # Add Honey bees manually at the end
  extract_model_parameters(Bumblebee_nbinom, 'Bumblebees'),
  extract_model_parameters(Solitary_bee_nbinom, 'Solitary bees')
)

# Filter for relevant terms and format the data
Margin_parameters_combined <- Margin_parameters_combined %>%
  filter(Parameter == "Strip_type1") %>%
  select(response, Coefficient, CI_low, CI_high, p, significance)

# Check the combined data frame
print(Margin_parameters_combined)

####### Save table
write.xlsx(Margin_parameters_combined, 'TableSX.xlsx')




### Figure 1: Model coefficients. Margin characterization & pollinator observations ##########

# Create a new column to indicate groups and set factor levels
Margin_parameters_combined <- Margin_parameters_combined %>%
  mutate(
    group = case_when( # new column with groupings to define colors and labels later
      response %in% c("Proportion vegetation cover", "Total floral area (cm2/m2)") ~ "Margin characterization",
      response %in% c("Hoverflies", "Honey bees", "Bumblebees", "Solitary bees") ~ "Pollinators",
      TRUE ~ "Other"  # Fallback category
    ),
    # Desired order for the groups and variables within groups
    group = factor(group, levels = c("Margin characterization", "Pollinators")),
    response = factor(response, levels = c(
      "Proportion vegetation cover", 
      "Total floral area (cm2/m2)",
      "Hoverflies", 
      "Honey bees", 
      "Bumblebees", 
      "Solitary bees"
    )),
    y_position = CI_high + 0.2  # Offset for significance labels
  )



# Create the plot with defined upper y limit and fixed significance label height
Modelcoefficients_V_PO <- ggplot(Margin_parameters_combined, aes(x = response, y = Coefficient)) +
  geom_point(aes(color = group, shape = group, fill = group), size = 5) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high, color = group), width = 0, linewidth = 1.5, show.legend = FALSE) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black", linewidth = 0.5) +
  geom_text(aes(y = y_position, label = significance), vjust = -1.5, color = "black", size = 4) +
  scale_color_manual(values = c(
    "Margin characterization" = "#7A4D9C", 
    "Pollinators" = "#FF9A00"
  )) +
  scale_fill_manual(values = c(
    "Margin characterization" = "#7A4D9C", 
    "Pollinators" = "#FF9A00"
  )) +
  scale_shape_manual(values = c(
    "Margin characterization" = 21, 
    "Pollinators" = 24
  )) +
  labs(x = NULL, y = "Model coefficients ± 95% CI\n (flower strip vs. control)") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 0, size = 15, color = "black"), 
        axis.text.y = element_text(angle = 0, size = 15, color = "black"),
        axis.title = element_text(size = 15),
        panel.border = element_rect(color = "black", fill = NA, linewidth = 1), 
        panel.grid = element_blank(),
        legend.text = element_text(size = 15),
        legend.position = "bottom") +
  scale_x_discrete(labels = c("Proportion vegetation cover" = "Vegetation cover",
                              "Total floral area (cm2/m2)" = "Floral area",
                              "Hoverflies" = "Hoverflies",
                              "Honey bees" = "Honey bees",
                              "Bumblebees" = "Bumblebees", 
                              "Solitary bees" = "Solitary bees")) +
  ylim(c(-1.5, 4.5)) +  # Set the specific upper Y limit
  guides(color = "none", shape = guide_legend(title = NULL), fill = guide_legend(title = NULL))

# Print the plot to visualize it
print(Modelcoefficients_V_PO)


# Save the plot as a PNG file
ggsave(filename = "Figure1.png", plot = Modelcoefficients_V_PO, width = 12, height = 6, units = "in", dpi = 300)








##    9.2. Margin & field habitats ---- 


### Function to run for Figure 2 and 3. Interaction plots for yellow sticky trap and pitfall models ##########
create_emm_plot <- function(model, response_variable, data, label_info, strip_colors = c("Control" = "#00AFBB", "Flower_strip" = "#E7B800")) {
  if (!(response_variable %in% colnames(data))) {
    stop(paste("Response variable", response_variable, "not found in data."))
  }
  
  emm <- emmeans(model, ~ Strip_type * Location)
  emm_df <- as.data.frame(emm)
  
  emm_df$emmean <- exp(emm_df$emmean)
  emm_df$asymp.LCL <- exp(emm_df$asymp.LCL)
  emm_df$asymp.UCL <- exp(emm_df$asymp.UCL)
  
  data$Location <- factor(data$Location, levels = c("Margin", "Crop"))
  emm_df$Location <- factor(emm_df$Location, levels = c("Margin", "Crop"))
  
  fitted_values <- predict(model, type = "response")
  fixed_effects <- exp(predict(model, type = "link", re.form = NA))
  observed_values <- data[[response_variable]]
  
  if (length(observed_values) != length(fitted_values)) {
    stop("Mismatch in lengths between observed and fitted values")
  }
  
  partial_residuals <- observed_values - fitted_values + fixed_effects
  data[[paste0("partial_residuals_", response_variable)]] <- partial_residuals
  
  dodge_width <- 0.2
  emm_df <- emm_df %>%
    mutate(Location_num = as.numeric(Location),
           dodged_x = Location_num + dodge_width * ifelse(Strip_type == "Control", -0.2, 0.2))
  
  ggplot() +
    geom_point(data = data, aes(x = Location, y = data[[paste0("partial_residuals_", response_variable)]], color = Strip_type),
               shape = 1, size = 5, position = position_jitter(width = 0.2)) +
    geom_point(data = emm_df, aes(x = dodged_x, y = emmean, color = Strip_type), shape = 16, size = 5) +
    geom_errorbar(data = emm_df, aes(x = dodged_x, ymin = asymp.LCL, ymax = asymp.UCL, color = Strip_type),
                  width = 0, linewidth = 2.5) +
    geom_line(data = emm_df, aes(x = dodged_x, y = emmean, group = Strip_type, color = Strip_type), linewidth = 1) +
    scale_color_manual(values = strip_colors, labels = c("Flower strip", "Control")) +
    labs(x = "Location", y = label_info$y_label, color = "Strip Type") +
    theme_classic() +
    theme(
      axis.text = element_text(size = 15),
      axis.title = element_text(size = 15),
      legend.text = element_text(size = 15),
      legend.title = element_text(size = 15),
      legend.position = "bottom",
      panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
    ) +
    annotate("text", x = -Inf, y = Inf, label = label_info$annotation, hjust = -0.1, vjust = 1.1, size = 5, fontface = "italic")
}



### Figure 2. Interaction plots for natural enemies in yellow sticky trap models  ##########

#######  Clarify terms to add all the models into the function

# Models list
YST_NE_modelslist <- list(
  YST_Spiders_negbinom,
  YST_Predatory_beetles_negbinom,
  YST_Predatory_heteroptera_negbinom,
  YST_Neuroptera_negbinom,
  YST_Predatory_thrips_negbinom, 
  YST_Syrphidae_negbinom,
  YST_Parasitic_wasps_negbinom
)

# Model names
YST_NE_modelsnames <- c(
  "Leaf-dwelling spiders",
  "Leaf-dwelling predatory beetles",
  "Predatory bugs",
  "Lacewings", 
  "Predatory thrips", 
  "Hoverflies",
  "Parasitic wasps"
)

# Variable names
YST_NE_response_mapping <- list(
  "Leaf-dwelling spiders" = "Spiders",
  "Leaf-dwelling predatory beetles" = "Predatory_beetles",
  "Predatory bugs" = "Predatory_heteroptera", 
  "Lacewings" = "Neuroptera",
  "Predatory thrips" = "Predatory_thrips", 
  "Hoverflies" = "Syrphidae",
  "Parasitic wasps" = "Parasitic_wasps"
)

# Labels and annotations
YST_NE_labelsinfo <- list(
  Spiders = list(y_label = "Leaf-dwelling \nspiders", annotation = "  treatment: n.s.\n  location: n.s.\ntreatment : location: n.s."),
  Predatory_beetles = list(y_label = "Leaf-dwelling \npredatory beetles", annotation = "  treatment: ***\n  location: (*)\ntreatment : location: (*)"),
  Predatory_heteroptera = list(y_label = "Predatory bugs", annotation = "  treatment: *\n  location: n.s.\ntreatment : location: n.s."),
  Neuroptera = list(y_label = "Lacewings", annotation = "  treatment: n.s.\n  location: n.s.\ntreatment : location: n.s."),
  Predatory_thrips = list(y_label = "Predatory thrips", annotation = "  treatment: n.s.\n  location: n.s.\ntreatment : location: n.s."),
  Syrphidae = list(y_label = "Hoverflies", annotation = "  treatment: n.s.\n  location: *\ntreatment : location: n.s."),
  Parasitic_wasps = list(y_label = "Parasitic wasps", annotation = "  treatment: ***\n  location: ***\ntreatment : location: **")
)


#Build dataframe for yellow sticky traps
YST_NE_dataframes <- rep(list(YST), length(YST_NE_modelslist))

#######  Generate plots
YST_NE_plots <- mapply(function(model, name, data) {
  response <- YST_NE_response_mapping[[name]]
  label_info <- YST_NE_labelsinfo[[response]]
  create_emm_plot(model, response, data, label_info)
}, model = YST_NE_modelslist, name = YST_NE_modelsnames, data = YST_NE_dataframes, SIMPLIFY = FALSE)


#######  Display plots
for (p in YST_NE_plots) {
  print(p)
}


#######  Combine all plots into a single figure
YST_NE_combined_plot <- ggarrange(
  plotlist = YST_NE_plots,
  ncol = 3, nrow = 3,
  labels = LETTERS[1:length(YST_NE_plots)],
  common.legend = TRUE,
  legend = "bottom"
)

# Display the combined plot
print(YST_NE_combined_plot)

# Save the figure with adjusted height for longer plots
ggsave("Figure2.png", YST_NE_combined_plot, width = 16, height = 18) 





### Figure 3. Interaction plots for herbivores in yellow sticky trap models  ##########

#######  Clarify terms to add all the models into the function

# Models list
YST_H_modelslist <- list(
  YST_Aphids_negbinom,
  YST_Cicadellidae_negbinom,
  YST_Lygaeidae_negbinom,
  YST_Nitidulidae_negbinom,
  YST_Psylloidea_negbinom,
  YST_Thysanoptera_negbinom
)

# Model names
YST_H_modelsnames <- c(
  "Aphids",
  "Leafhoppers",
  "Seed bugs",
  "Sap beetles", 
  "Psyllids",
  "Thrips"
  )

# Variable names
YST_H_response_mapping <- list(
  "Aphids" = "Aphids",
  "Leafhoppers" = "Cicadellidae",
  "Seed bugs" = "Lygaeidae",
  "Sap beetles" = "Nitidulidae",
  "Psyllids" = "Psylloidea",
  "Thrips" = "Thysanoptera"
)


# Labels and annotations
YST_H_labelsinfo <- list(
  Aphids = list(y_label = "Aphids", annotation = "  treatment: *\n  location: n.s.\ntreatment : location: n.s."),
  Cicadellidae = list(y_label = "Leafhoppers", annotation = "  treatment: ***\n  location: n.s.\ntreatment : location: n.s."),
  Lygaeidae= list(y_label = "Seed bugs", annotation = "  treatment: (*)\n  location: n.s.\ntreatment : location: n.s."),
  Nitidulidae= list(y_label = "Sap beetles", annotation = "  treatment: ***\n  location: \ntreatment : location: ***"),
  Psylloidea= list(y_label = "Psyllids", annotation = "  treatment: n.s.\n  location: **\ntreatment : location: n.s."),
  Thysanoptera = list(y_label = "Thrips", annotation = "  treatment: ***\n  location: n.s.\ntreatment : location: *")
  )


#Build dataframe for yellow sticky traps
YST_H_dataframes <- rep(list(YST), length(YST_H_modelslist))

#######  Generate plots
YST_H_plots <- mapply(function(model, name, data) {
  response <- YST_H_response_mapping[[name]]
  label_info <- YST_H_labelsinfo[[response]]
  create_emm_plot(model, response, data, label_info)
}, model = YST_H_modelslist, name = YST_H_modelsnames, data = YST_H_dataframes, SIMPLIFY = FALSE)


#######  Display plots
for (p in YST_H_plots) {
  print(p)
}


#######  Combine all plots into a single figure
YST_H_combined_plot <- ggarrange(
  plotlist = YST_H_plots,
  ncol = 3, nrow = 2,
  labels = LETTERS[1:length(YST_H_plots)],
  common.legend = TRUE,
  legend = "bottom"
)

# Display the combined plot
print(YST_H_combined_plot)

# Save the figure with adjusted height for longer plots
ggsave("Figure3.png", YST_H_combined_plot, width = 16, height = 12) 





###Figure 4. Interaction plots for pitfall trap models ##########


#######  Clarify terms to add all the models into the function

# Models list
Pitfall_modelslist <- list(
  Pitfall_Spiders_negbinom,
  Pitfall_Carabidae_negbinom,
  Pitfall_Staphylinidae_negbinom
)

# Display names for models
Pitfall_modelsnames <- c(
  "Ground-dwelling spiders",
  "Carabid beetles",
  "Rove beetles"
)

# Variable names
Pitfall_response_mapping <- list(
  "Ground-dwelling spiders" = "Ground_Spiders",
  "Carabid beetles" = "Carabidae",
  "Rove beetles" = "Staphylinidae"
)

# Labels and annotations
Pitfall_labelsinfo <- list(
  Ground_Spiders = list( y_label = "Ground-dwelling \nspiders", annotation = "  treatment: **\n  location: n.s.\ntreatment : location: (*)"),
  Carabidae = list(y_label = "Carabid beetles", annotation = "  treatment: **\n  location: n.s.\ntreatment : location: n.s."),
  Staphylinidae = list(y_label = "Rove beetles", annotation = "  treatment: ***\n  location: n.s.\ntreatment : location: n.s.")
)


# Create dataframe for pitfall
Pitfall_dataframes <- rep(list(Pitfall), length(Pitfall_modelslist))


#######  Generate plots 
Pitfall_plots <- mapply(function(model, name, data) {
  response <- Pitfall_response_mapping[[name]]
  label_info <- Pitfall_labelsinfo[[response]]
  create_emm_plot(model, response, data, label_info)
}, model = Pitfall_modelslist, name = Pitfall_modelsnames, data = Pitfall_dataframes, SIMPLIFY = FALSE)

# Display plots
for (p in Pitfall_plots) {
  print(p)
}



#######  Combine all plots into a single figure
Pitfall_combined_plot <- ggarrange(
  plotlist = Pitfall_plots,
  ncol = 3, nrow = 1,
  labels = LETTERS[1:length(Pitfall_plots)],
  common.legend = TRUE,
  legend = "bottom"
)

# Display the combined plot
print(Pitfall_combined_plot)

# Save the figure with adjusted height for longer plots
ggsave("Figure4.png", Pitfall_combined_plot, width = 16, height = 6) 






###Figure 6. Interaction plot for decomposition rate ##########
# Extract estimated marginal means
emm_bl <- emmeans(BL_binomial, ~ Strip_type * Location)
emm_df_bl <- as.data.frame(emm_bl) %>%
  mutate(
    emmean = exp(emmean) / (1 + exp(emmean)),  # Back-transform from logit to probability scale
    asymp.LCL = exp(asymp.LCL) / (1 + exp(asymp.LCL)),
    asymp.UCL = exp(asymp.UCL) / (1 + exp(asymp.UCL))
  )

# Calculate partial residuals on the link scale (logit scale) and back-transform
link_fitted_bl <- predict(BL_binomial, type = "link")  # includes random effects
link_fixed_bl  <- predict(BL_binomial, type = "link", re.form = NA)  # fixed only

# Logit transformation: log(p / (1 - p)), so residuals on the link scale need this adjustment
partial_resid_link_bl <- log(BLS$Decomposition_rate / (1 - BLS$Decomposition_rate + 1e-3)) - link_fitted_bl + link_fixed_bl

# Back-transform residuals from logit scale to probability scale
BLS$partial_resids_bl <- 1 / (1 + exp(-partial_resid_link_bl))


# Adjust EMMs for dodging in the plot
dodge_width <- 0.2
emm_df_bl <- emm_df_bl %>%
  mutate(
    Location_num = as.numeric(Location),
    dodged_x = Location_num + dodge_width * ifelse(Strip_type == "Flower_strip", -0.2, 0.2)
  )

# Plot
Decomposition_plot <- ggplot() +
  geom_point(data = BLS,
             aes(x = Location, y = partial_resids_bl, color = Strip_type),
             shape = 1, size = 5, position = position_jitter(width = 0.2)) +
  geom_point(data = emm_df_bl,
             aes(x = dodged_x, y = emmean, color = Strip_type),
             shape = 16, size = 5) +
  geom_errorbar(data = emm_df_bl,
                aes(x = dodged_x, ymin = asymp.LCL, ymax = asymp.UCL, color = Strip_type),
                width = 0, linewidth = 2.5) +
  geom_line(data = emm_df_bl,
            aes(x = dodged_x, y = emmean, group = Strip_type, color = Strip_type),
            linewidth = 1) +
  scale_color_manual(values = c("Control" = "#00AFBB", "Flower_strip" = "#E7B800"), 
                     labels = c("Flower strip", "Control")) + 
  labs(x = "Location", y = "Decomposition rate", color = "Strip type") +
  theme_classic() +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 15),
    legend.text = element_text(size = 15),
    legend.title = element_text(size = 15),
    legend.position = "bottom",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
  ) +
  annotate("text", x = -Inf, y = Inf,
           label = "  treatment: n.s.\n  location: **\ntreatment : location: (*)", 
           hjust = -0.1, vjust = 1.1, size = 5, fontface = "italic")

# Print the plot
print(Decomposition_plot)



# Save the figure with adjusted height for longer plots
ggsave("Figure6.png", Decomposition_plot, width = 6, height = 6) 










##    9.3. Field habitats ---- 

### Table SXX: Preliminary table that gives the model coefficients necessary for Figure 3 ##########
# Define a function to extract model parameters and add significance
extract_model_parameters <- function(model, response_name) {
  as.data.frame(model_parameters(model, type = 'response')) %>%
    mutate(response = response_name,
           significance = case_when(
             p < 0.001 ~ "***",
             p < 0.01  ~ "**",
             p < 0.05  ~ "*",
             p < 0.1   ~ "(*)",
             TRUE      ~ ""
           ))
}

# Apply the function to each model and combine the results
Field_parameters_combined <- bind_rows(
  extract_model_parameters(TC_NE_negbinom, 'Natural enemies'),
  extract_model_parameters(TC_Thrips_negbinom, 'Thrips'),
  extract_model_parameters(TC_Aphids_poisson, 'Aphids'),
  extract_model_parameters(AC_leaf_betabinomial, 'Aphid predation rates by leaf-dwellers'),
  extract_model_parameters(AC_ground_betabinomial, 'Aphid predation rates by ground-dwellers')
)

# Filter for relevant terms and format the data
Field_parameters_combined <- Field_parameters_combined %>%
  filter(Parameter == "Strip_type1") %>%
  select(response, Coefficient, CI_low, CI_high, p, significance)

# Check the combined data frame
View(Field_parameters_combined)
print(Field_parameters_combined)

####### Save table
write.xlsx(Field_parameters_combined, 'TableSXX.xlsx')






### Figure 5: Model coefficients. Tiller counts & aphid cards  ##########

# Create a new column to indicate groups and set factor levels
Field_parameters_combined <- Field_parameters_combined %>%
  mutate(
    group = case_when( # new column with groupings to define colors and labels later
      response %in% c("Natural enemies", "Thrips", "Aphids") ~ "Tiller counts",
      response %in% c("Aphid predation rates by leaf-dwellers", "Aphid predation rates by ground-dwellers") ~ "Aphid cards",
      TRUE ~ "Other"  # Fallback category
    ),
    # Desired order for the groups and variables within groups
    group = factor(group, levels = c("Tiller counts", "Aphid cards")),
    response = factor(response, levels = c(
      "Natural enemies",
      "Thrips", 
      "Aphids",
      "Aphid predation rates by leaf-dwellers",
      "Aphid predation rates by ground-dwellers"
    )),
    y_position = CI_high + 0.2  # Offset for significance labels
  )



# Create the plot with defined upper y limit and fixed significance label height
Model_coefficients_TC_AC <- ggplot(Field_parameters_combined, aes(x = response, y = Coefficient)) +
  geom_point(aes(color = group, shape = group, fill = group), size = 5) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high, color = group), width = 0, linewidth = 1.5, show.legend = FALSE) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black", linewidth = 0.5) +
  geom_text(aes(y = y_position, label = significance), vjust = -1.5, color = "black", size = 4) +
  scale_color_manual(values = c(
    "Tiller counts" = "#7A4D9C", 
    "Aphid cards" = "#FF9A00"
  )) +
  scale_fill_manual(values = c(
    "Tiller counts" = "#7A4D9C", 
    "Aphid cards" = "#FF9A00"
  )) +
  scale_shape_manual(values = c(
    "Tiller counts" = 21, 
    "Aphid cards" = 24
  )) +
  labs(x = NULL, y = "Model coefficients ± 95% CI\n (flower strip vs. control)") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 0, size = 15, color = "black"), 
        axis.text.y = element_text(angle = 0, size = 15, color = "black"),
        axis.title = element_text(size = 15),
        panel.border = element_rect(color = "black", fill = NA, linewidth = 1), 
        panel.grid = element_blank(),
        legend.text = element_text(size = 15),
        legend.position = "bottom") +
  scale_x_discrete(labels = c(" Natural enemies" = " Natural enemies",
                              "Thrips" = "Thrips",
                              "Aphids" = "Aphids",
                              "Aphid predation rates by leaf-dwellers" = "Aphid predation rates \nby leaf-dwellers",
                              "Aphid predation rates by ground-dwellers" = "Aphid predation rates \nby ground-dwellers")) +
  ylim(c(-1, 1)) +  # Set the specific upper Y limit
  guides(color = "none", shape = guide_legend(title = NULL), fill = guide_legend(title = NULL))

# Print the plot to visualize it
print(Model_coefficients_TC_AC)


# Save the plot as a PNG file
ggsave(filename = "Figure5.png", plot = Model_coefficients_TC_AC, width = 12, height = 6, units = "in", dpi = 300)







# 10. Supplementary material ----

### Table S5. Posthoc for crop ##########

#### S5.1. Margin & field habitats ##########
####### Model list
MarginField_modelslist_S51 <- list(
  YST_Parasitic_wasps_negbinom,
  YST_Aphids_negbinom,
  YST_Nitidulidae_negbinom,
  YST_Psylloidea_negbinom
)

####### Model names
MarginField_modelsnames_S51 <- c(
  "Parasitic wasps",
  "Aphids",
  "Sap beetles", 
  "Psyllids"
)



# Create function to extract and transform EMMs
extract_emm_and_pairs_FieldMargin_S51 <- function(models, names, factor = ~ Crop) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  # Combine and return
  return(list(
    FieldMargin_EMMs_S51 = bind_rows(all_emms),
    FieldMargin_Pairs_S51 = bind_rows(all_pairs)
  ))
}



# Run function 
FieldMargin_crop_results_S51 <- extract_emm_and_pairs_FieldMargin_S51(MarginField_modelslist_S51, MarginField_modelsnames_S51)

# View results
print(FieldMargin_crop_results_S51$FieldMargin_EMMs_S51)
print(FieldMargin_crop_results_S51$FieldMargin_Pairs_S51)




#### S5.2. Field habitats ##########
####### Model list
Field_modelslist_S52 <- list(
  TC_Aphids_poisson,
  AC_leaf_betabinomial,
  AC_ground_betabinomial
)

####### Model names
Field_modelsnames_S52 <- c(
  "Aphids", 
  "Aphid predation rates by leaf-dwellers", 
  "Aphid predation rates by ground-dwellers"
)

####### Function to get emmeans and CI
extract_emm_and_pairs_Field_S52 <- function(models, names, factor = ~ Crop) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  # Combine and return
  return(list(
    Field_EMMs_S52 = bind_rows(all_emms),
    Field_Pairs_S52 = bind_rows(all_pairs)
  ))
}

# Run function 
field_results_S52 <- extract_emm_and_pairs_Field_S52(Field_modelslist_S52, Field_modelsnames_S52)

# View results
print(field_results_S52$Field_EMMs_S52)
print(field_results_S52$Field_Pairs_S52)

























 
### Figure S1. Interaction (Treatment*Crop) ##########


#### S1.1.Psyllids in yellow sticky traps ##########

# Extract estimated marginal means
emm_psy <- emmeans(YST_Psylloidea_negbinom, ~ Strip_type * Crop)
emm_df_psy <- as.data.frame(emm_psy) %>%
  mutate(
    emmean = exp(emmean),
    asymp.LCL = exp(asymp.LCL),
    asymp.UCL = exp(asymp.UCL)
  )

# Calculate partial residuals on link scale and back-transform
link_fitted_psy <- predict(YST_Psylloidea_negbinom, type = "link")
link_fixed_psy  <- predict(YST_Psylloidea_negbinom, type = "link", re.form = NA)
partial_resid_link_psy <- log(YST$Psylloidea + 1e-3) - link_fitted_psy + link_fixed_psy
YST$partial_resids_psy <- exp(partial_resid_link_psy)


# Plot
Psyllids_treatment_crop_plot <- ggplot() +
  geom_point(data = YST,
             aes(x = Crop, y = partial_resids_psy, color = Strip_type),
             shape = 1, size = 5, position = position_jitter(width = 0.2)) +
  geom_point(data = emm_df_psy,
             aes(x = Crop, y = emmean, color = Strip_type),
             shape = 16, size = 5, position = position_dodge(width = 0.4)) +
  geom_errorbar(data = emm_df_psy,
                aes(x = Crop, ymin = asymp.LCL, ymax = asymp.UCL, color = Strip_type),
                width = 0, linewidth = 2.5, position = position_dodge(width = 0.4)) +
  geom_line(data = emm_df_psy,
            aes(x = Crop, y = emmean, group = Strip_type, color = Strip_type),
            linewidth = 1, position = position_dodge(width = 0.4)) +
  scale_color_manual(values = c("Control" = "#00AFBB", "Flower_strip" = "#E7B800"), 
                     labels = c("Flower strip", "Control")) +  # Rename legend labels
  scale_x_discrete(labels = c("Spring_barley" = "Spring barley", "Winter_wheat" = "Winter wheat")) +  # Re-label x-axis
  labs(x = "Crop", y = "Psyllids", color = "Strip type") +
  theme_classic() +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 15),
    legend.text = element_text(size = 15),
    legend.title = element_text(size = 15),
    legend.position = "bottom",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
  ) +
  annotate("text", x = -Inf, y = Inf,
           label = "  treatment: n.s.\n  crop: n.s.\n treatment : crop: **",
           hjust = -0.1, vjust = 1.1, size = 5, fontface = "italic")


# Print the plot to visualize it
print(Psyllids_treatment_crop_plot)



# Save the figure with adjusted height for longer plots
ggsave("FigureS1.1.png", Psyllids_treatment_crop_plot, width = 6, height = 6) 






#### S1.2.Ground-dwelling spiders in pitfall traps ##########

# Extract estimated marginal means
emm_sp <- emmeans(Pitfall_Spiders_negbinom, ~ Strip_type * Crop)
emm_df_sp <- as.data.frame(emm_sp) %>%
  mutate(
    emmean = exp(emmean),
    asymp.LCL = exp(asymp.LCL),
    asymp.UCL = exp(asymp.UCL)
  )

# Calculate partial residuals on link scale and back-transform
link_fitted_sp <- predict(Pitfall_Spiders_negbinom, type = "link")
link_fixed_sp  <- predict(Pitfall_Spiders_negbinom, type = "link", re.form = NA)
partial_resid_link_sp <- log(Pitfall$Ground_Spiders + 1e-3) - link_fitted_sp + link_fixed_sp
Pitfall$partial_resids_sp <- exp(partial_resid_link_sp)


# Plot
Ground_Spiders_treatment_crop_plot <- ggplot() +
  geom_point(data = Pitfall,
             aes(x = Crop, y = partial_resids_sp, color = Strip_type),
             shape = 1, size = 5, position = position_jitter(width = 0.2)) +
  geom_point(data = emm_df_sp,
             aes(x = Crop, y = emmean, color = Strip_type),
             shape = 16, size = 5, position = position_dodge(width = 0.4)) +
  geom_errorbar(data = emm_df_sp,
                aes(x = Crop, ymin = asymp.LCL, ymax = asymp.UCL, color = Strip_type),
                width = 0, linewidth = 2.5, position = position_dodge(width = 0.4)) +
  geom_line(data = emm_df_sp,
            aes(x = Crop, y = emmean, group = Strip_type, color = Strip_type),
            linewidth = 1, position = position_dodge(width = 0.4)) +
  scale_color_manual(values = c("Control" = "#00AFBB", "Flower_strip" = "#E7B800"), 
                     labels = c("Flower strip", "Control")) +  # Rename legend labels
  scale_x_discrete(labels = c("Spring_barley" = "Spring barley", "Winter_wheat" = "Winter wheat")) +  # Re-label x-axis
  labs(x = "Crop", y = "Ground-dwelling spiders", color = "Strip type") +
  theme_classic() +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 15),
    legend.text = element_text(size = 15),
    legend.title = element_text(size = 15),
    legend.position = "bottom",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
  ) +
  annotate("text", x = -Inf, y = Inf,
           label = "  treatment: n.s.\n  crop: n.s.\n treatment : crop: **",
           hjust = -0.1, vjust = 1.1, size = 5, fontface = "italic")


# Print the plot to visualize it
print(Ground_Spiders_treatment_crop_plot)



# Save the figure with adjusted height for longer plots
ggsave("FigureS1.2.png", Ground_Spiders_treatment_crop_plot, width = 6, height = 6) 








#######  Combine all plots into a single figure #######  
Treatment_Crop_combined_plot <- ggarrange(
  plotlist = list(Psyllids_treatment_crop_plot, Ground_Spiders_treatment_crop_plot),
  ncol = 2, nrow = 1,
  labels = c("A", "B"),
  common.legend = TRUE,
  legend = "bottom"
)

# Display the combined plot
print(Treatment_Crop_combined_plot)

# Save the figure with adjusted height for longer plots
ggsave("FigureS1.png", Treatment_Crop_combined_plot, width = 10.6, height = 6) 





#######  Post-hoc analisis for the interacitons  #######  
####### Model list
MarginField_modelslist <- list(
  YST_Spiders_negbinom,
  YST_Predatory_beetles_negbinom,
  YST_Predatory_heteroptera_negbinom,
  YST_Neuroptera_negbinom,
  YST_Predatory_thrips_negbinom,
  YST_Syrphidae_negbinom,
  YST_Parasitic_wasps_negbinom,
  YST_Aphids_negbinom,
  YST_Cicadellidae_negbinom,
  YST_Lygaeidae_negbinom,
  YST_Nitidulidae_negbinom,
  YST_Psylloidea_negbinom,
  YST_Thysanoptera_negbinom,
  Pitfall_Spiders_negbinom,
  Pitfall_Carabidae_negbinom,
  Pitfall_Staphylinidae_negbinom,
  BL_binomial
)

####### Model names
MarginField_modelsnames <- c(
  "Leaf-dwelling spiders",
  "Leaf-dwelling predatory beetles", 
  "Predatory bugs", 
  "Lacewings", 
  "Predatory thrips", 
  "Hoverflies",
  "Parasitic wasps",
  "Aphids",
  "Leafhoppers",
  "Seed bugs",
  "Sap beetles", 
  "Psyllids",
  "Thrips",
  "Ground-dwelling spiders",
  "Carabid beetles",
  "Rove beetles",
  "Decompostion rate"
)



####### Function to extract and transform EMMs and pairwise comparisons
extract_emm_and_pairs_FieldMargin_CropTreatment <- function(models, names, factor = ~ Crop * Strip_type) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  return(list(
    FieldMargin_CropTreatment_EMMs = bind_rows(all_emms),
    FieldMargin_CropTreatment_Pairs = bind_rows(all_pairs)
  ))
}

# Run function 
fieldmargin_results_CropTreatment <- extract_emm_and_pairs_FieldMargin_CropTreatment(MarginField_modelslist, MarginField_modelsnames)

# View results
print(fieldmargin_results_CropTreatment$FieldMargin_CropTreatment_EMMs)
print(fieldmargin_results_CropTreatment$FieldMargin_CropTreatment_Pairs)











### Figure S2. Interaction (Treatment*Crop)  ##########


#### S2. Aphids in tiller counts ##########
# Extract estimated marginal means
emm_atc <- emmeans(TC_Aphids_poisson, ~ Strip_type * Crop)
emm_df_atc <- as.data.frame(emm_atc) %>%
  mutate(
    emmean = exp(emmean),
    asymp.LCL = exp(asymp.LCL),
    asymp.UCL = exp(asymp.UCL)
  )

# Calculate partial residuals on link scale and back-transform
link_fitted_atc <- predict(TC_Aphids_poisson, type = "link")
link_fixed_atc  <- predict(TC_Aphids_poisson, type = "link", re.form = NA)
partial_resid_link_atc <- log(TC$Aphids + 1e-3) - link_fitted_atc + link_fixed_atc
TC$partial_resids_atc <- exp(partial_resid_link_atc)


# Plot
TC_aphids_treatment_crop_plot <- ggplot() +
  geom_point(data = TC,
             aes(x = Crop, y = partial_resids_atc, color = Strip_type),
             shape = 1, size = 5, position = position_jitter(width = 0.2)) +
  geom_point(data = emm_df_atc,
             aes(x = Crop, y = emmean, color = Strip_type),
             shape = 16, size = 5, position = position_dodge(width = 0.4)) +
  geom_errorbar(data = emm_df_atc,
                aes(x = Crop, ymin = asymp.LCL, ymax = asymp.UCL, color = Strip_type),
                width = 0, linewidth = 2.5, position = position_dodge(width = 0.4)) +
  geom_line(data = emm_df_atc,
            aes(x = Crop, y = emmean, group = Strip_type, color = Strip_type),
            linewidth = 1, position = position_dodge(width = 0.4)) +
  scale_color_manual(values = c("Control" = "#00AFBB", "Flower_strip" = "#E7B800"), 
                     labels = c("Flower strip", "Control")) +  # Rename legend labels
  scale_x_discrete(labels = c("Spring_barley" = "Spring barley", "Winter_wheat" = "Winter wheat")) +  # Re-label x-axis
  labs(x = "Crop", y = "Aphids", color = "Strip type") +
  theme_classic() +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 15),
    legend.text = element_text(size = 15),
    legend.title = element_text(size = 15),
    legend.position = "bottom",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
  ) +
  annotate("text", x = -Inf, y = Inf,
           label = "  treatment: n.s.\n   crop: *\n treatment : crop: ***",
           hjust = -0.1, vjust = 1.1, size = 5, fontface = "italic")


# Print the plot to visualize it
print(TC_aphids_treatment_crop_plot)



# Save the figure with adjusted height for longer plots
ggsave("FigureS2.png", TC_aphids_treatment_crop_plot, width = 6, height = 6) 





#### Post-hoc test ##########
####### Model list
Field_modelslist <- list(
  TC_NE_negbinom,
  TC_Thrips_negbinom,
  TC_Aphids_poisson
)

####### Model names
Field_modelsnames <- c(
  "Natural enemies",
  "Thrips", 
  "Aphids"
)


####### Function to get emmeans and CI
extract_emm_and_pairs_Field_interaction <- function(models, names, factor = ~ Strip_type * Crop) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  # Combine and return
  return(list(
    Field_interaction_EMMs = bind_rows(all_emms),
    Field_interaction_Pairs = bind_rows(all_pairs)
  ))
}

# Run function 
field_interaction_results <- extract_emm_and_pairs_Field_interaction(Field_modelslist, Field_modelsnames)

# View results
print(field_interaction_results$Field_interaction_EMMs)
print(field_interaction_results$Field_interaction_Pairs)









### Figure S3. Interaction (Crop*Location)   ##########


#### S3.1. Parasitic wasps in yellow sticky traps   ##########

# Extract estimated marginal means
emm_pw <- emmeans(YST_Parasitic_wasps_negbinom, ~ Location * Crop)
emm_df_pw <- as.data.frame(emm_pw) %>%
  mutate(
    emmean = exp(emmean),
    asymp.LCL = exp(asymp.LCL),
    asymp.UCL = exp(asymp.UCL)
  )

# Calculate partial residuals on link scale and back-transform
link_fitted_pw <- predict(YST_Parasitic_wasps_negbinom, type = "link")
link_fixed_pw  <- predict(YST_Parasitic_wasps_negbinom, type = "link", re.form = NA)
partial_resid_link_pw <- log(YST$Parasitic_wasps + 1e-3) - link_fitted_pw + link_fixed_pw
YST$partial_resids_pw <- exp(partial_resid_link_pw)

# Plot
Parasitic_wasps_crop_location_plot <- ggplot() +
  geom_point(data = YST,
             aes(x = Location, y = partial_resids_pw, color = Crop),
             shape = 1, size = 5, position = position_jitter(width = 0.2)) +
  geom_point(data = emm_df_pw,
             aes(x = Location, y = emmean, color = Crop),
             shape = 16, size = 5, position = position_dodge(width = 0.4)) +
  geom_errorbar(data = emm_df_pw,
                aes(x = Location, ymin = asymp.LCL, ymax = asymp.UCL, color = Crop),
                width = 0, linewidth = 2.5, position = position_dodge(width = 0.4)) +
  geom_line(data = emm_df_pw,
            aes(x = Location, y = emmean, group = Crop, color = Crop),
            linewidth = 1, position = position_dodge(width = 0.4)) +
  scale_color_manual(values = c("Spring_barley" = "#E66100", "Winter_wheat" = "#5D3A9B"),
                     labels = c("Spring barley", "Winter wheat")) +
  labs(x = "Location", y = "Parasitic wasps", color = "Crop") +
  theme_classic() +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 15),
    legend.text = element_text(size = 15),
    legend.title = element_text(size = 15),
    legend.position = "bottom",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
  ) +
  annotate("text", x = -Inf, y = Inf,
           label = "  location: ***\n  crop: *\n  location : crop: *",
           hjust = -0.1, vjust = 1.1, size = 5, fontface = "italic")


# Print the plot to visualize it
print(Parasitic_wasps_crop_location_plot)



# Save the figure with adjusted height for longer plots
ggsave("FigureS3.1.png", Parasitic_wasps_crop_location_plot, width = 6, height = 6) 




#### S3.2. Seed bugs in yellow sticky traps   ##########

# Extract estimated marginal means
emm_lg <- emmeans(YST_Lygaeidae_negbinom, ~ Location * Crop)
emm_df_lg <- as.data.frame(emm_lg) %>%
  mutate(
    emmean = exp(emmean),
    asymp.LCL = exp(asymp.LCL),
    asymp.UCL = exp(asymp.UCL)
  )

# Calculate partial residuals on link scale and back-transform
link_fitted_lg <- predict(YST_Lygaeidae_negbinom, type = "link")
link_fixed_lg  <- predict(YST_Lygaeidae_negbinom, type = "link", re.form = NA)
partial_resid_link_lg <- log(YST$Lygaeidae + 1e-3) - link_fitted_lg + link_fixed_lg
YST$partial_resids_lg <- exp(partial_resid_link_lg)

# Plot
Seed_bugs_crop_location_plot <- ggplot() +
  geom_point(data = YST,
             aes(x = Location, y = partial_resids_lg, color = Crop),
             shape = 1, size = 5, position = position_jitter(width = 0.2)) +
  geom_point(data = emm_df_lg,
             aes(x = Location, y = emmean, color = Crop),
             shape = 16, size = 5, position = position_dodge(width = 0.4)) +
  geom_errorbar(data = emm_df_lg,
                aes(x = Location, ymin = asymp.LCL, ymax = asymp.UCL, color = Crop),
                width = 0, linewidth = 2.5, position = position_dodge(width = 0.4)) +
  geom_line(data = emm_df_lg,
            aes(x = Location, y = emmean, group = Crop, color = Crop),
            linewidth = 1, position = position_dodge(width = 0.4)) +
  scale_color_manual(values = c("Spring_barley" = "#E66100", "Winter_wheat" = "#5D3A9B"),
                     labels = c("Spring barley", "Winter wheat")) +
  labs(x = "Location", y = "Seed bugs", color = "Crop") +
  theme_classic() +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 15),
    legend.text = element_text(size = 15),
    legend.title = element_text(size = 15),
    legend.position = "bottom",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
  ) +
  annotate("text", x = -Inf, y = Inf,
           label = "  location: n.s.\n  crop: n.s.\n  location : crop: *",
           hjust = -0.1, vjust = 1.1, size = 5, fontface = "italic")


# Print the plot to visualize it
print(Seed_bugs_crop_location_plot)



# Save the figure with adjusted height for longer plots
ggsave("FigureS3.2.png", Seed_bugs_crop_location_plot, width = 6, height = 6) 




#### S3.3. Rove beetles in pitfall traps   ##########

# Extract estimated marginal means
emm_rv <- emmeans(Pitfall_Staphylinidae_negbinom, ~ Location * Crop)
emm_df_rv <- as.data.frame(emm_rv) %>%
  mutate(
    emmean = exp(emmean),
    asymp.LCL = exp(asymp.LCL),
    asymp.UCL = exp(asymp.UCL)
  )

# Calculate partial residuals on link scale and back-transform
link_fitted_rv <- predict(Pitfall_Staphylinidae_negbinom, type = "link")
link_fixed_rv  <- predict(Pitfall_Staphylinidae_negbinom, type = "link", re.form = NA)
partial_resid_link_rv <- log(Pitfall$Staphylinidae + 1e-3) - link_fitted_rv + link_fixed_rv
Pitfall$partial_resids_rv <- exp(partial_resid_link_rv)


# Plot
Staphylinidae_crop_location_plot <- ggplot() +
  geom_point(data = Pitfall,
             aes(x = Location, y = partial_resids_rv, color = Crop),
             shape = 1, size = 5, position = position_jitter(width = 0.2)) +
  geom_point(data = emm_df_rv,
             aes(x = Location, y = emmean, color = Crop),
             shape = 16, size = 5, position = position_dodge(width = 0.4)) +
  geom_errorbar(data = emm_df_rv,
                aes(x = Location, ymin = asymp.LCL, ymax = asymp.UCL, color = Crop),
                width = 0, linewidth = 2.5, position = position_dodge(width = 0.4)) +
  geom_line(data = emm_df_rv,
            aes(x = Location, y = emmean, group = Crop, color = Crop),
            linewidth = 1, position = position_dodge(width = 0.4)) +
  scale_color_manual(values = c("Spring_barley" = "#E66100", "Winter_wheat" = "#5D3A9B"),
                     labels = c("Spring barley", "Winter wheat")) +
  labs(x = "Location", y = "Rove beetles", color = "Crop") +
  theme_classic() +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 15),
    legend.text = element_text(size = 15),
    legend.title = element_text(size = 15),
    legend.position = "bottom",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
  ) +
  annotate("text", x = -Inf, y = Inf,
           label = "  location: n.s.\n  crop: n.s.\n location : crop: **",
           hjust = -0.1, vjust = 1.1, size = 5, fontface = "italic")


# Print the plot to visualize it
print(Staphylinidae_crop_location_plot)



# Save the figure with adjusted height for longer plots
ggsave("FigureS3.3.png", Staphylinidae_crop_location_plot, width = 6, height = 6) 




#######  Combine all plots into a single figure #######  
Crop_Location_combined_plot <- ggarrange(
  plotlist = list(Parasitic_wasps_crop_location_plot, Seed_bugs_crop_location_plot, Staphylinidae_crop_location_plot),
  ncol = 3, nrow = 1,
  labels = c("A", "B", "C"),
  common.legend = TRUE,
  legend = "bottom"
)

# Display the combined plot
print(Crop_Location_combined_plot)

# Save the figure with adjusted height for longer plots
ggsave("FigureS3.png", Crop_Location_combined_plot, width = 16, height = 6) 




#######  Post-hoc analisis for the interactions  #######  
####### Model list
MarginField_modelslist_S3 <- list(
  YST_Parasitic_wasps_negbinom,
  YST_Lygaeidae_negbinom,
  Pitfall_Staphylinidae_negbinom
)

####### Model names
MarginField_modelsnames_S3 <- c(
  "Parasitic wasps",
  "Seed bugs",
  "Rove beetles"
)


####### Function to extract and transform EMMs and pairwise comparisons
extract_emm_and_pairs_FieldMargin_S3 <- function(models, names, factor = ~ Crop * Location) {
  all_emms <- list()
  all_pairs <- list()
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    name <- names[[i]]
    
    # Get the link function
    link_fun <- tryCatch({
      family(model)$link
    }, error = function(e) {
      model$family$link
    })
    
    # Get EMMs
    emm <- emmeans(model, factor)
    emm_df <- as.data.frame(emm)
    
    # Back-transform if needed
    if (link_fun == "log") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), exp))
    } else if (link_fun == "logit") {
      emm_df <- emm_df %>%
        mutate(across(any_of(c("emmean", "asymp.LCL", "asymp.UCL")), plogis))
    }
    emm_df$Model <- name
    
    # Pairwise comparisons
    contrast_df <- as.data.frame(pairs(emm, adjust = "tukey"))
    if ("estimate" %in% names(contrast_df)) {
      if (link_fun == "log") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), exp))
      } else if (link_fun == "logit") {
        contrast_df <- contrast_df %>%
          mutate(across(any_of(c("estimate", "lower.CL", "upper.CL")), plogis))
      }
    }
    contrast_df$Model <- name
    
    all_emms[[name]] <- emm_df
    all_pairs[[name]] <- contrast_df
  }
  
  return(list(
    FieldMargin_EMMs_S3 = bind_rows(all_emms),
    FieldMargin_Pairs_S3 = bind_rows(all_pairs)
  ))
}

# Run function 
fieldmargin_results_S3 <- extract_emm_and_pairs_FieldMargin_S3(MarginField_modelslist_S3, MarginField_modelsnames_S3)

# View results
print(fieldmargin_results_S3$FieldMargin_EMMs_S3)
print(fieldmargin_results_S3$FieldMargin_Pairs_S3)