## This R Script details the data analysis process described in the associated paper
## This script can be followed to repeat the analysis and recreate the main figures

# To protect privacy of the participating farmers we are unable to publicly share any geodata associated with the project.
# However, this information can be requested by contacting the corresponding author.

# Required packages
library(tidyverse)
library(lme4)
library(glmmTMB)
library(performance)
library(car)
library(vegan)
library(ggplot2)
library(ggpubr)
library(MuMIn)
library(DHARMa)
library(piecewiseSEM)
library("GGally")

############################## DATA ANALYSIS ##############################################

############ LOAD IN PEST AND NATURAL ENEMY OBSERVATION DATASETS #####

PP <- read.csv("Pest_Pressure.csv") # Pest pressure observations
AD <- read.csv("Adult_Damage .csv") # Adult feeding damage observations
LL <- read.csv("Larva_Damage.csv") # Larval load observations
TD <- read.csv("Insect_Traps.csv") # Natural enemy observations (trap data)

# Note there are some NAs in the larval load dataset  where the null observations (0s) were not input
LL[7:10][is.na(LL[7:10])] <- 0 # Replaces the NAs with the null observations

# Update column types and specify which should be considered factors
PP$Field <- as.factor(PP$Field)
PP$Transect <- as.factor(PP$Transect)
PP$Week <- as.factor(PP$Week)

AD$Field <- as.factor(AD$Field)
AD$Transect <- as.factor(AD$Transect)
AD$Assessment_Area <- as.factor(AD$Assessment_Area)
AD$Assessment_Round <- as.factor(AD$Assessment_Round)
AD$Percent_Damage_Logit <- (logit(AD$Percent_Damage)) # logit transoformation

LL$Field <- as.factor(LL$Field)
LL$Transect <- as.factor(LL$Transect)
LL$Assessment_Area <- as.factor(LL$Assessment_Area)
LL$Assessment_Round <- 2 # recode assessment round so data will merge with the cumulative pest pressure data
LL$Assessment_Round <- as.factor(LL$Assessment_Round)

TD$Field <- as.factor(TD$Field)
TD$Transect <- as.factor(TD$Transect)
TD$Assessment_Round <- as.factor(TD$Assessment_Round)
TD$Trap_Type <- as.factor(TD$Trap_Type)

# Load in agronomic/adjacent OSR dataset
Adjacent_OSR <- read.csv("Rapeseed_Site_Info.csv") # This data includes the agronomic information for each site alongside Y/N for adjacent rapeseed and distance from previous season rapeseed field
Adjacent_OSR$Field <- as.factor(Adjacent_OSR$Field)

# Merge the pest pressure, damage, and larvae data with the agronomic datafiles and check for NAs
PP.A <- full_join(Adjacent_OSR, PP)
which(is.na(PP.A), arr.ind=TRUE) # Note - this code checks for NA values. This is a good commonsense check for successful merging of datasets, as unsuccessful merging can create additional rows and NA values
AD.A <- full_join(Adjacent_OSR, AD)
which(is.na(AD.A), arr.ind=TRUE) # 
LL.A <- full_join(Adjacent_OSR, LL)
which(is.na(LL.A), arr.ind=TRUE) # 

############ POOL data to reduce zero-inflation. Data pooling procedure is described in the manuscript #####
##### Pest pressure
PP.Pool <- PP.A %>%
  group_by(Field, Transect, Week) %>%
  summarise(Pressure_Cummulative = sum(Pressure_Cummulative), Prox_OSR_2021 = max(Distance_OSR_2021_m),
            Adj_OSR_2022 = max(Adj_OSR_2022),
            .groups = "keep") # This pooled data is the cumulative beetle counts across all weeks / visits split at field-transect level

which(is.na(PP.Pool), arr.ind=TRUE) #

##### Adult damage
AD.Pool <- AD.A %>%
  group_by(Field, Transect, Assessment_Round, Assessment_Area) %>%
  summarise(Percent_Damage_Logit = mean(Percent_Damage_Logit), Percent_Damage = mean(Percent_Damage),
            Prox_OSR_2021 = max(Distance_OSR_2021_m),
            Adj_OSR_2022 = max(Adj_OSR_2022),
            .groups = "keep")

which(is.na(AD.Pool), arr.ind=TRUE) # 

##### Larval load 
LL.Pool <- LL.A %>%
  group_by(Field, Transect, Assessment_Area) %>%
  summarise(Larvae_Total_Mean = mean(Larvae_Total),
            Larvae_Total_Sum = sum(Larvae_Total),
            Plant_Diameter = mean(Plant_Diameter),Prox_OSR_2021 = max(Distance_OSR_2021_m),
            Adj_OSR_2022 = max(Adj_OSR_2022),
            .groups = "keep")

which(is.na(LL.Pool), arr.ind=TRUE) # 

###### Natural enemy
# Trap data has one sample point per transect, but data should be summed across trap types
# Trap data pooling/summing
TD.Pool <- TD %>%
  group_by(Field, Transect, Assessment_Round) %>%
  summarise_at(vars(Coleoptera:Hym_Apidea), sum, na.rm = TRUE)
which(is.na(TD.Pool), arr.ind=TRUE) # 

# Calculate the abundance of potential natural enemies by summing - Carabids, Spiders, Braconid wasps, Ichnumonid wasps
TD.Pool$PotentialNatEn <- rowSums(TD.Pool[, c(18,26,27,28,29,30,31,32,33)]) # CHECK THE COLUMN NUMBERS CORRESPOND TO THE DESIRED ARTHROPOD GROUPS
# Calculate Shannon's diversity for potential natural enemy diversity
TD.Pool$PotentialNatEnDiv <- diversity(TD.Pool[, c(18,26,27,28,29,30,31,32,33)], index = "shannon", base = exp(1)) # CHECK THE COLUMN NUMBERS CORRESPOND TO THE DESIRED ARTHROPOD GROUPS

# Check for correlation between Natural enemy diversity and abundance
cor.test(TD.Pool$PotentialNatEn, TD.Pool$PotentialNatEnDiv) # No correlation

###################################### DATA ANALYSIS: Pest pressure, adult damage, larvae load, natural enemy ######
## Merge data for damage correlations. This creates the datasets for the initial analysis
PP.Pool.Sub <- subset(PP.Pool, Week %in% c ("39", "43")) # Subset pest pressure data to initial pressure and final cummulative pressure (week 39 and 43)
PP.Pool.Sub$Assessment_Round <- PP.Pool.Sub$Week # Recode the week number to the corresponding assessment round
PP.Pool.Sub$Assessment_Round <- recode(PP.Pool.Sub$Assessment_Round, "'39' = '1';
                                       '43' = '2'")
#
PP.AD.Pool <- full_join(x = PP.Pool.Sub, y = AD.Pool) # merge pest pressure and adult damage data for assessment of the influence pest pressure has on adult damage

# For Larvae Load only final cummulative pest pressure data used
PP.Sub.LL <- subset(PP.Pool.Sub, Assessment_Round %in% c ("2"))
PP.LL.Pool <- full_join(x = PP.Sub.LL, y = LL.Pool) # merge pest pressure and larvae load for assessment of the influence pest pressure has on larvae load

Dam1 <- subset(PP.AD.Pool, Assessment_Round %in% c ("1")) # Assessment round 1 - First damage assessment
which(is.na(Dam1), arr.ind=TRUE) # 
Dam2 <- subset(PP.AD.Pool, Assessment_Round %in% c ("2")) # Assessment round 2 - Second damage assessment
which(is.na(Dam2), arr.ind=TRUE) #
Dam3 <- subset(PP.LL.Pool, Assessment_Round %in% c ("2")) # Assessment round 3 - Larvae load assessment
which(is.na(Dam3), arr.ind=TRUE) # 

##### Data analysis - relationships between beetle abundance and damage/larvae

# Damage assessment 1
Damage.Pressure.a <- lmer(Percent_Damage_Logit ~ Pressure_Cummulative+Transect+
                            (1|Field),
                          data = Dam1)
summary(Damage.Pressure.a) # summary of the statistical test
plot(Damage.Pressure.a) # Plot model and check distributions
Anova(Damage.Pressure.a) # Analysis of deviance test for significance
vif(Damage.Pressure.a) # Check the variance inflation factors for all variables
AICc(Damage.Pressure.a) # Extract the AICc values

# Damage assessment 2
Damage.Pressure.b <- lmer(Percent_Damage_Logit ~ Pressure_Cummulative+Transect+
                            (1|Field),
                   data = Dam2)
summary(Damage.Pressure.b)
plot(Damage.Pressure.b)
Anova(Damage.Pressure.b)
vif(Damage.Pressure.b)
AICc(Damage.Pressure.b)

# Larvae load
Larvae.Pressure.c.2 <- glmmTMB(Larvae_Total_Sum ~ Pressure_Cummulative+Transect+
                                 +(1|Field),  
              family="nbinom2",
              data=PP.LL.Pool)
summary(Larvae.Pressure.c.2)
Anova(Larvae.Pressure.c.2)
check_collinearity(Larvae.Pressure.c.2) # This code calculated VIF in negative binomial models
AICc(Larvae.Pressure.c.2)
plot(DHARMa::simulateResiduals(Larvae.Pressure.c.2)) # DHARMa package is used to plot residuals for negative binomial models

##### Data analysis - impact of OSR proximity and adjacency on pest pressure, adult damage, and larvae load
# Pest pressure
PP.OSR.a.2 <- glmmTMB(Pressure_Cummulative ~ Adj_OSR_2022+sqrt(Prox_OSR_2021)+Transect+Assessment_Round+(1|Field),  
                               family="nbinom2",
                               data=PP.Pool.Sub)
summary(PP.OSR.a.2)
Anova(PP.OSR.a.2)
check_collinearity(PP.OSR.a.2)
AICc(PP.OSR.a.2)
plot(DHARMa::simulateResiduals(PP.OSR.a.2))

# Crop damage
AD.OSR.a <- lmer(Percent_Damage_Logit ~ Adj_OSR_2022+sqrt(Prox_OSR_2021)+Transect+Assessment_Round+(1|Field),
                 data = AD.Pool)
summary(AD.OSR.a)
plot(AD.OSR.a)
Anova(AD.OSR.a)
vif(AD.OSR.a)
AICc(AD.OSR.a)

# Larvae load
LL.OSR.a.2 <- glmmTMB(Larvae_Total_Sum ~ Adj_OSR_2022+sqrt(Prox_OSR_2021)+Transect+(1|Field),  
                      family="nbinom2",
                      data=LL.Pool)
summary(LL.OSR.a.2)
Anova(LL.OSR.a.2)
check_collinearity(LL.OSR.a.2)
AICc(LL.OSR.a.2)
plot(DHARMa::simulateResiduals(LL.OSR.a.2))

##### Data analysis: Impact of natural enemies on pest pressure, adult damage, and larvae load
NatEnPressDam <- subset(TD.Pool, Assessment_Round %in% c ("1", "2")) # Subset to the different assessment rounds. Round 1 and 2 for analysis with pest pressure and damage data
NatEmLar <- subset(TD.Pool, Assessment_Round %in% c ("3")) # Trap data assessment round 3 for analysis with larvae data

# merge trap data with pressure, damage, and larvae datasets - for this analysis, need to pool across transects
PP.Nat.Pool <- full_join(x = PP.Pool.Sub, y = NatEnPressDam) # already transect-level
which(is.na(TD.Pool), arr.ind=TRUE) # 

#
AD.Pool.Trans <- AD %>%
  group_by(Field, Transect, Assessment_Round) %>%
  summarise(Percent_Damage_Logit = mean(Percent_Damage_Logit), Percent_Damage = mean(Percent_Damage),
            .groups = "keep")
AD.Nat.Pool <- full_join(x = AD.Pool.Trans, y = NatEnPressDam)
which(is.na(AD.Nat.Pool), arr.ind=TRUE) # 

#
LL.Pool.Trans <- LL %>%
  group_by(Field, Transect) %>%
  summarise(Larvae_Total_Sum = sum(Larvae_Total),
            Plant_Diameter = mean(Plant_Diameter),
            .groups = "keep")
LL.Pool.Trans$Assessment_Round <- 3
LL.Pool.Trans$Assessment_Round <- as.factor(LL.Pool.Trans$Assessment_Round)
LL.Nat.Pool <- full_join(x = LL.Pool.Trans, y = NatEmLar)
which(is.na(LL.Nat.Pool), arr.ind=TRUE) #

# Run the analysis to examine how pest pressure, adult damage, and larvae load are affected by the natural enemies

# Pest pressure
PP.Nat.a.2 <- glmmTMB(Pressure_Cummulative ~ PotentialNatEn+PotentialNatEnDiv+Transect+Assessment_Round+
                        (1|Field), 
                      family="nbinom2",
                      data=PP.Nat.Pool)
summary(PP.Nat.a.2)
Anova(PP.Nat.a.2)
check_collinearity(PP.Nat.a.2)
AICc(PP.Nat.a.2)
plot(DHARMa::simulateResiduals(PP.Nat.a.2))

# Crop damage
AD.Nat.a <- lmer(Percent_Damage_Logit ~ PotentialNatEn+PotentialNatEnDiv+Transect+
                   Assessment_Round+(1|Field), 
                 data = AD.Nat.Pool)
summary(AD.Nat.a)
plot(AD.Nat.a)
Anova(AD.Nat.a)
vif(AD.Nat.a)
AICc(AD.Nat.a)

# Larvae load
LL.Nat.a.2 <- glmmTMB(Larvae_Total_Sum ~ PotentialNatEn+PotentialNatEnDiv+Transect+
                        (1|Field),  
                      family="nbinom2",
                      data=LL.Nat.Pool)
summary(LL.Nat.a.2)
Anova(LL.Nat.a.2)
check_collinearity(LL.Nat.a.2)
AICc(LL.Nat.a.2)
plot(DHARMa::simulateResiduals(LL.Nat.a.2))


############################## LANDSCAPE ANALYSIS ##########################################
##### LOAD IN AND CREATE LANDSCAPE DATASETS #####

# Load in the landscape datasets. There are three datasets in total, one for each spatial scale.
# Datasets have been extracted from processed QGIS data
# To protect privacy of the participating farmers, the associated geodata are unable to be publicly shared.

# 500 m radii
Landscape1 <- read.csv("Landscape_500m_Radii.csv")
Landscape1$Field_ID
Landscape1$Field <- recode(Landscape1$Field_ID, "'Feld 16' = '16'; 'Feld 17' = '17'; 'Feld 18' = '18';'Feld 2' = '2';
                             'Feld 20' = '20';'Feld 21' = '21';'Feld 25' = '25';'Feld 27' = '27';
                             'Feld 3' = '3';'Feld 30' = '30';'Feld 31' = '31';'Feld 32' = '32';
                             'Feld 5' = '5';'Feld 6' = '6'") # Recode the field codes so they match across the datasets
Landscape1$Field <- as.factor(Landscape1$Field)

Landscape500m <- Landscape1 # Rename the dataset

# 750 m radii
Landscape2 <- read.csv("Landscape_750m_Radii.csv")
Landscape2$Field_ID
Landscape2$Field <- recode(Landscape2$Field_ID, "'Feld 16' = '16'; 'Feld 17' = '17'; 'Feld 18' = '18';'Feld 2' = '2';
                             'Feld 20' = '20';'Feld 21' = '21';'Feld 25' = '25';'Feld 27' = '27';
                             'Feld 3' = '3';'Feld 30' = '30';'Feld 31' = '31';'Feld 32' = '32';
                             'Feld 5' = '5';'Feld 6' = '6'") # Recode the field codes so they match across the datasets
Landscape2$Field <- as.factor(Landscape2$Field)
Landscape2$Field

Landscape750 <- Landscape2 # Rename the dataset

# 1000 m radii
Landscape3 <- read.csv("Landscape_1000m_Radii.csv")
Landscape3$Field_ID
Landscape3$Field <- recode(Landscape3$Field_ID, "'Feld 16' = '16'; 'Feld 17' = '17'; 'Feld 18' = '18';'Feld 2' = '2';
                             'Feld 20' = '20';'Feld 21' = '21';'Feld 25' = '25';'Feld 27' = '27';
                             'Feld 3' = '3';'Feld 30' = '30';'Feld 31' = '31';'Feld 32' = '32';
                             'Feld 5' = '5';'Feld 6' = '6'")# Recode the field codes so they match across the datasets
Landscape3$Field <- as.factor(Landscape3$Field)

Landscape1000 <- Landscape3 # Rename the dataset

###### Merge landscape data with observation data #####
# Combine the field data and the landscape data into several further datasets: One for each observation-landscape combination

Land.PP.500 <- full_join(PP.Pool.Sub, Landscape500m) # landscape - 500 m; pest pressure
Land.PP.500$Transect <- as.factor(Land.PP.500$Transect)
Land.PP.500$Adj_OSR_2022 <- as.factor(Land.PP.500$Adj_OSR_2022)
which(is.na(Land.PP.500), arr.ind=TRUE) # This code checks for NA values. It is also a good common sense check when merging data, if NAs appear it suggests the two datasets might not have merged correctly

Land.AD.500 <- full_join(Landscape500m, AD.Pool) # Landscape - 500 m; adult damage
Land.AD.500$Transect <- as.factor(Land.AD.500$Transect)
Land.AD.500$Adj_OSR_2022 <- as.factor(Land.AD.500$Adj_OSR_2022)
Land.AD.500$Assessment_Round <- as.factor(Land.AD.500$Assessment_Round)
which(is.na(Land.AD.500), arr.ind=TRUE) # 

Land.LL.500 <- full_join(Landscape500m, LL.Pool) # Landscape - 500 m; larvae abundance
Land.LL.500$Transect <- as.factor(Land.LL.500$Transect)
Land.LL.500$Adj_OSR_2022 <- as.factor(Land.LL.500$Adj_OSR_2022)
which(is.na(Land.LL.500), arr.ind=TRUE) # 
#

Land.PP.750 <- full_join(PP.Pool.Sub, Landscape750) # Landscape - 750 m; pest pressure
Land.PP.750$Transect <- as.factor(Land.PP.750$Transect)
Land.PP.750$Adj_OSR_2022 <- as.factor(Land.PP.750$Adj_OSR_2022)
which(is.na(Land.PP.750), arr.ind=TRUE) # 

Land.AD.750 <- full_join(Landscape750, AD.Pool) # Landscape - 750 m; adult damage
Land.AD.750$Transect <- as.factor(Land.AD.750$Transect)
Land.AD.750$Adj_OSR_2022 <- as.factor(Land.AD.750$Adj_OSR_2022)
Land.AD.750$Assessment_Round <- as.factor(Land.AD.750$Assessment_Round)
which(is.na(Land.AD.750), arr.ind=TRUE) #

Land.LL.750 <- full_join(Landscape750, LL.Pool) # Landscape - 750 m; larvae abundance
Land.LL.750$Transect <- as.factor(Land.LL.750$Transect)
Land.LL.750$Adj_OSR_2022 <- as.factor(Land.LL.750$Adj_OSR_2022)
which(is.na(Land.LL.750), arr.ind=TRUE) # 
#

Land.PP.1000 <- full_join(Landscape1000, PP.Pool.Sub) # Landscape - 1000 m; pest pressure
Land.PP.1000$Transect <- as.factor(Land.PP.1000$Transect)
Land.PP.1000$Adj_OSR_2022 <- as.factor(Land.PP.1000$Adj_OSR_2022)
which(is.na(Land.PP.1000), arr.ind=TRUE) # 

Land.AD.1000 <- full_join(Landscape1000, AD.Pool) # Landscape - 1000 m; adult damage
Land.AD.1000$Transect <- as.factor(Land.AD.1000$Transect)
Land.AD.1000$Adj_OSR_2022 <- as.factor(Land.AD.1000$Adj_OSR_2022)
Land.AD.1000$Assessment_Round <- as.factor(Land.AD.1000$Assessment_Round)
which(is.na(Land.AD.1000), arr.ind=TRUE) # 

Land.LL.1000 <- full_join(Landscape1000, LL.Pool) # Landscape - 1000 m; larvae abundance
Land.LL.1000$Transect <- as.factor(Land.LL.1000$Transect)
Land.LL.1000$Adj_OSR_2022 <- as.factor(Land.LL.1000$Adj_OSR_2022)
which(is.na(Land.LL.1000), arr.ind=TRUE) # 

### Natural enemy data
NatEnPressDam <- subset(TD.Pool, Assessment_Round %in% c ("1", "2")) # Subset to the different assessment rounds. Round 1 and 2 for analysis with pest pressure and damage data
NatEmLar <- subset(TD.Pool, Assessment_Round %in% c ("3")) # Trap data assessment round 3 for analysis with larvae data



###### Data analysis - Pest pressure #####
# 500 m model

PP.Land.500.a <- glmmTMB(Pressure_Cummulative ~ Prop_500m_Y2021_WinterRape+
                Prop_500m_Y2022_WinterRape+
                CropDiv21+CropDiv22+SNH_Prop_500m+Transect+Assessment_Round+(1|Field),  
              family="nbinom2",
              data=Land.PP.500)
summary(PP.Land.500.a)
car::Anova(PP.Land.500.a)
check_collinearity(PP.Land.500.a)
AICc(PP.Land.500.a)
plot(DHARMa::simulateResiduals(PP.Land.500.a))

# 750 m model

PP.Land.750.b <- glmmTMB(Pressure_Cummulative ~ Prop_750m_Y2021_WinterRape+
                           Prop_750m_Y2022_WinterRape+
                CropDiv21+CropDiv22+Prop_SNH_750m+Transect+Assessment_Round+(1|Field),  
              family="nbinom2",
              data=Land.PP.750)
summary(PP.Land.750.b)
car::Anova(PP.Land.750.b)
check_collinearity(PP.Land.750.b)
AICc(PP.Land.750.b)
plot(DHARMa::simulateResiduals(PP.Land.750.b))

# 1000 m model

PP.Land.1000.a.2 <- glmmTMB(Pressure_Cummulative ~ Prop_1000m_Y2021_WinterRape+
                              Prop_1000m_Y2022_WinterRape+
                           CropDiv21+CropDiv22+Prop_SNH_1000m+Transect+Assessment_Round+(1|Field),  
                         family="nbinom2",
                         data=Land.PP.1000)
summary(PP.Land.1000.a.2)
car::Anova(PP.Land.1000.a.2)
check_collinearity(PP.Land.1000.a.2) # High co-linearity/VIF value for 2021 crop heterogeneity, remove variable from analysis
AICc(PP.Land.1000.a.2)
#
PP.Land.1000.a.3 <- glmmTMB(Pressure_Cummulative ~ Prop_1000m_Y2021_WinterRape+
                              Prop_1000m_Y2022_WinterRape+
                              CropDiv22+Prop_SNH_1000m+Transect+Assessment_Round+(1|Field),  
                            family="nbinom2",
                            data=Land.PP.1000)
summary(PP.Land.1000.a.3)
car::Anova(PP.Land.1000.a.3)
check_collinearity(PP.Land.1000.a.3)
AICc(PP.Land.1000.a.3)
plot(DHARMa::simulateResiduals(PP.Land.1000.a.3))

###### Data analysis - adult damage ######

# 500 m radius model

AD.Land.500.a <- lmer(Percent_Damage_Logit ~ Prop_500m_Y2021_WinterRape+
                        Prop_500m_Y2022_WinterRape+
                        CropDiv21+CropDiv22+SNH_Prop_500m+Transect+Assessment_Round+(1|Field),
                  data = Land.AD.500)
summary(AD.Land.500.a)
plot(AD.Land.500.a)
car::Anova(AD.Land.500.a)
vif(AD.Land.500.a)
AICc(AD.Land.500.a)

# 750 m radius model

AD.Land.750.a <- lmer(Percent_Damage_Logit ~ Prop_750m_Y2021_WinterRape+
                        Prop_750m_Y2022_WinterRape+
                        CropDiv21+CropDiv22+Prop_SNH_750m+Transect+Assessment_Round+(1|Field),
                      data = Land.AD.750)
summary(AD.Land.750.a)
plot(AD.Land.750.a)
car::Anova(AD.Land.750.a)
vif(AD.Land.750.a)
AICc(AD.Land.750.a)

# 1000 m radius model

AD.Land.1000.a <- lmer(Percent_Damage_Logit ~ Prop_1000m_Y2021_WinterRape+
                         Prop_1000m_Y2022_WinterRape+
                        CropDiv21+CropDiv22+Prop_SNH_1000m+Assessment_Round+(1|Field),
                      data = Land.AD.1000)
summary(AD.Land.1000.a)
plot(AD.Land.1000.a)
car::Anova(AD.Land.1000.a)
vif(AD.Land.1000.a) # High VIF for 2021 crop heterogeneity; remove variable
AICc(AD.Land.1000.a)

#
AD.Land.1000.b <- lmer(Percent_Damage_Logit ~ Prop_1000m_Y2021_WinterRape+
                         Prop_1000m_Y2022_WinterRape+
                         CropDiv22+Prop_SNH_1000m+Transect+Assessment_Round+(1|Field),
                       data = Land.AD.1000)
summary(AD.Land.1000.b)
plot(AD.Land.1000.b)
car::Anova(AD.Land.1000.b)
vif(AD.Land.1000.b)
AICc(AD.Land.1000.b)
###### Data analysis - larvae load #####

# 500 m radius model

LL.Land.500.a.2 <- glmmTMB(Larvae_Total_Sum ~ Prop_500m_Y2021_WinterRape+
                             Prop_500m_Y2022_WinterRape+
                              CropDiv21+CropDiv22+SNH_Prop_500m+Transect+(1|Field),  
                            family="nbinom2",
                            data=Land.LL.500)
summary(LL.Land.500.a.2)
Anova(LL.Land.500.a.2)
check_collinearity(LL.Land.500.a.2)
AICc(LL.Land.500.a.2)
plot(DHARMa::simulateResiduals(LL.Land.500.a.2))

# 750 m radius model

LL.Land.500.a.3 <- glmmTMB(Larvae_Total_Sum ~ Prop_750m_Y2021_WinterRape+
                             Prop_750m_Y2022_WinterRape+
                             CropDiv21+CropDiv22+Prop_SNH_750m+Transect+(1|Field),  
                           family="nbinom2",
                           data=Land.LL.750)
summary(LL.Land.500.a.3)
Anova(LL.Land.500.a.3)
check_collinearity(LL.Land.500.a.3)
AICc(LL.Land.500.a.3)
plot(DHARMa::simulateResiduals(LL.Land.500.a.3))

# 1000 m radius model

LL.Land.1000.a.2 <- glmmTMB(Larvae_Total_Sum ~ Prop_1000m_Y2021_WinterRape+
                              Prop_1000m_Y2022_WinterRape+
                             CropDiv21+CropDiv22+Prop_SNH_1000m+Transect+(1|Field),  
                           family="nbinom2",
                           data=Land.LL.1000)
summary(LL.Land.1000.a.2)
car::Anova(LL.Land.1000.a.2)
check_collinearity(LL.Land.1000.a.2) # High VIF/co-linearity for 2021 crop heterogeneity, remove and repeat analysis
AICc(LL.Land.1000.a.2)

# 
LL.Land.1000.a.3 <- glmmTMB(Larvae_Total_Sum ~ Prop_1000m_Y2021_WinterRape+
                              Prop_1000m_Y2022_WinterRape+
                              CropDiv22+Prop_SNH_1000m+Transect+(1|Field),  
                            family="nbinom2",
                            data=Land.LL.1000)
summary(LL.Land.1000.a.3)
Anova(LL.Land.1000.a.3)
check_collinearity(LL.Land.1000.a.3)
AICc(LL.Land.1000.a.3)
plot(DHARMa::simulateResiduals(LL.Land.1000.a.3))

##### Data analysis - natural enemy #####
# Create the natural enemy landscape datasets, similar to the processes above for the pest pressure, damage, and larvae data

Land.NE.500 <- full_join(TD.Pool, Landscape500m) # Landscape - 500 m
which(is.na(Land.NE.500), arr.ind=TRUE) # 

Land.NE.750 <- full_join(TD.Pool, Landscape750) # Landscape - 500 m; larvae abundance
which(is.na(Land.NE.750), arr.ind=TRUE) #

Land.NE.1000 <- full_join(Landscape1000, TD.Pool) # Landscape - 500 m; larvae abundance
which(is.na(Land.NE.1000), arr.ind=TRUE) #

## Natural enemy abudnance models
# 500 m spatial models
NE.Land.500.a.2 <- glmmTMB(PotentialNatEn ~ Prop_500m_Y2021_WinterRape+
                             Prop_500m_Y2022_WinterRape+
                             CropDiv21+CropDiv22+SNH_Prop_500m+Transect+Assessment_Round+(1|Field), 
                      family="nbinom2",
                      data=Land.NE.500)
summary(NE.Land.500.a.2)
Anova(NE.Land.500.a.2)
check_collinearity(NE.Land.500.a.2)
AICc(NE.Land.500.a.2)
plot(DHARMa::simulateResiduals(NE.Land.500.a.2))

#

NE.Land.750.a.2 <- glmmTMB(PotentialNatEn ~ Prop_750m_Y2021_WinterRape+
                             Prop_750m_Y2022_WinterRape+
                             CropDiv21+CropDiv22+Prop_SNH_750m+Transect+Assessment_Round+(1|Field), 
                           family="nbinom2",
                           data=Land.NE.750)
summary(NE.Land.750.a.2)
Anova(NE.Land.750.a.2)
check_collinearity(NE.Land.750.a.2)
AICc(NE.Land.750.a.2)
plot(DHARMa::simulateResiduals(NE.Land.750.a.2))

#
NE.Land.1000.a.2 <- glmmTMB(PotentialNatEn ~ Prop_1000m_Y2021_WinterRape+
                             Prop_1000m_Y2022_WinterRape+
                             CropDiv21+CropDiv22+Prop_SNH_1000m+Transect+Assessment_Round+(1|Field), 
                           family="nbinom2",
                           data=Land.NE.1000)
summary(NE.Land.1000.a.2)
Anova(NE.Land.1000.a.2)
check_collinearity(NE.Land.1000.a.2) # High VIF for 2021 crop heterogeneity
AICc(NE.Land.1000.a.2)
plot(DHARMa::simulateResiduals(NE.Land.1000.a.2))

#
NE.Land.1000.a.3 <- glmmTMB(PotentialNatEn ~ Prop_1000m_Y2021_WinterRape +
                              Prop_1000m_Y2022_WinterRape+
                              CropDiv22+Prop_SNH_1000m+Transect+Assessment_Round+(1|Field), 
                            family="nbinom2",
                            data=Land.NE.1000)
summary(NE.Land.1000.a.3)
Anova(NE.Land.1000.a.3)
check_collinearity(NE.Land.1000.a.3)
AICc(NE.Land.1000.a.3)
plot(DHARMa::simulateResiduals(NE.Land.1000.a.3))

#### Natural enemy diversity analysis

NE.Land.500.b <- lmer(PotentialNatEnDiv ~ Prop_500m_Y2021_WinterRape+
                        Prop_500m_Y2022_WinterRape+
                        CropDiv21+CropDiv22+SNH_Prop_500m+Transect+Assessment_Round+(1|Field),
                      data = Land.NE.500)
summary(NE.Land.500.b)
plot(NE.Land.500.b)
car::Anova(NE.Land.500.b)
vif(NE.Land.500.b)
AICc(NE.Land.500.b)


#
NE.Land.750.b <- lmer(PotentialNatEnDiv ~ Prop_750m_Y2021_WinterRape+
                        Prop_750m_Y2022_WinterRape+
                        CropDiv21+CropDiv22+Prop_SNH_750m+Transect+Assessment_Round+(1|Field),
                      data = Land.NE.750)
summary(NE.Land.750.b)
plot(NE.Land.750.b)
car::Anova(NE.Land.750.b)
vif(NE.Land.750.b)
AICc(NE.Land.750.b)

#
NE.Land.1000.c <- lmer(PotentialNatEnDiv ~ Prop_1000m_Y2021_WinterRape+
                         Prop_1000m_Y2022_WinterRape+
                        CropDiv21+CropDiv22+Prop_SNH_1000m+Transect+Assessment_Round+(1|Field),
                      data = Land.NE.1000)
summary(NE.Land.1000.c)
plot(NE.Land.1000.c)
car::Anova(NE.Land.1000.c)
vif(NE.Land.1000.c) # High VIF for 2021 heterogeneity
AIC(NE.Land.1000.c)

NE.Land.1000.d <- lmer(PotentialNatEnDiv ~ Prop_1000m_Y2021_WinterRape+
                         Prop_1000m_Y2022_WinterRape+
                         CropDiv22+Prop_SNH_1000m+Transect+Assessment_Round+(1|Field),
                       data = Land.NE.1000)
summary(NE.Land.1000.d)
plot(NE.Land.1000.d)
car::Anova(NE.Land.1000.d)
vif(NE.Land.1000.d)
AIC(NE.Land.1000.d)


########################################### END OF STANDARD DATA ANALYSIS ####################################

########################################## CODE FOR PLOTS ###########################################
# Plot themese, axis names etc.,.
PlotTheme <- theme_classic(base_size = 9)
Jitter_Theme <- geom_jitter(size = 0.5, alpha = 0.15)
Boxplot_Theme <- geom_boxplot(outlier.alpha = 0.00)
GLM_Line <- geom_smooth(method = "glm", formula = y~x,
                        method.args = list(family = negative.binomial(theta = 1, link = log)),
                        size = 0.5)
LM_Line <- geom_smooth(method = "lm", formula = y~x, size = 0.5)
PressureLab39 <- "Beetle pressure (week 39)"
PressureLabCumm <- "Cumulative pressure"
DamageLab <- "Leaf damage (%)"
LarvaeLab <- "Larvae per quadrat"

SameYearRapsTitle <- "Rapeseed proportion: Cropping year"
PreviousYearHet <- "Crop heterogeneity: Previous year"
SamesYearHet <- "Crop heterogeneity: Cropping year"
PropSNHTitle <- "Seminatural habitat proprtion"

Assessement_Label <- c("Round 1", "Round 2")
Assessement_LabelNE <- c("Round 1", "Round 2", "Round 3")
names(Assessement_Label) <- c("1", "2")
names(Assessement_LabelNE) <- c("1", "2", "3")
AssessmentFacet <- facet_wrap(~Assessment_Round, labeller = labeller(Assessment_Round = Assessement_Label))
NatAbLabel <- c("Total abundance of natural enemy groups")
AssessmentFacetNE <- facet_wrap(~Assessment_Round, labeller = labeller(Assessment_Round = Assessement_LabelNE))

##### FIG 3 #####
# Panel A-C: Impact of beetle pressure on adult damage and larvae load
PressDam <- ggplot(Dam1, aes(y = Percent_Damage, x = Pressure_Cummulative))+
  Jitter_Theme+
  PlotTheme+ylab(DamageLab)+
  xlab(PressureLab39)+LM_Line+ggtitle("***")
PressDam

PressDam2 <- ggplot(Dam2, aes(y = Percent_Damage, x = Pressure_Cummulative))+
  Jitter_Theme+
  PlotTheme+ylab(DamageLab)+
  xlab(PressureLabCumm)+LM_Line+ggtitle("**")
PressDam2

PressDam3 <- ggplot(Dam3, aes(y = Larvae_Total_Sum, x = Pressure_Cummulative))+
  Jitter_Theme+
  PlotTheme+ylab(LarvaeLab)+
  xlab(PressureLabCumm)+GLM_Line+ggtitle("***")
PressDam3

# Panel D-F: Impact of an adjacent OSR field on pressure, damage, and larvae load
Beetle.OSR.A <- ggplot(Land.PP.500, aes(y = Pressure_Cummulative, x = Adj_OSR_2022))+
  Boxplot_Theme+Jitter_Theme+PlotTheme+ylab(PressureLabCumm)+
  xlab(element_blank())+ggtitle("***")
Beetle.OSR.A

Beetle.OSR.B <- ggplot(Land.AD.500, aes(y = Percent_Damage, x = Adj_OSR_2022))+
  Boxplot_Theme+Jitter_Theme+
  PlotTheme+ylab(DamageLab)+
  xlab("Adjacent to rapeseed crop")+ggtitle("**")
Beetle.OSR.B

Beetle.OSR.C <- ggplot(Land.LL.500, aes(y = Larvae_Total_Sum, x = Adj_OSR_2022))+
  Boxplot_Theme+Jitter_Theme+PlotTheme+ylab(LarvaeLab)+
  xlab(element_blank())+ggtitle("*")
Beetle.OSR.C

Fig.3 <- ggarrange(PressDam, PressDam2, PressDam3, Beetle.OSR.A, Beetle.OSR.B, Beetle.OSR.C,
                   labels = c("A", "B", "C", "D", "E", "F"), 
                            font.label = list(size = 8), nrow = 2, ncol = 3)
Fig.3
#
ggsave(
  "Fig_3.tiff",
  plot = Fig.3,
  device = "tiff",
  scale = 1,
  width = 135,
  height = 100,
  units = c("mm"),
  dpi = 300,
  limitsize = TRUE)


##### FIG 4 ######
# Landscape plots: Effect of rapeseed proportion on beetle pressure and adult damage

Beetle.A <- ggplot(Land.PP.500, aes(y = Pressure_Cummulative, x = Prop_500m_Y2022_WinterRape))+
  Jitter_Theme+PlotTheme+ylab(PressureLabCumm)+
  theme(legend.position = "none")+ylim(NA,125)+
  xlab(SameYearRapsTitle)+AssessmentFacet+
  GLM_Line+ggtitle("*")
Beetle.A

# adult damage
Beetle.B <- ggplot(Land.AD.500, aes(y = Percent_Damage, x = Prop_500m_Y2022_WinterRape))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  ylab(DamageLab)+
  xlab(SameYearRapsTitle)+LM_Line+ggtitle("*")
Beetle.B

Fig.4 <- ggarrange(Beetle.A, Beetle.B,
                   labels = c("AUTO"), common.legend = TRUE, legend = "bottom",
                   font.label = list(size = 8), nrow = 2, ncol = 1)
Fig.4

#
#
ggsave(
  "Fig_4.tiff",
  plot = Fig.4,
  device = "tiff",
  scale = 1,
  width = 80,
  height = 140,
  units = c("mm"),
  dpi = 300,
  limitsize = TRUE)

##### Fig 5 #####
# Landscape plots: Effect of rapeseed proportion and crop heterogeneity on larvae load

Larvae.A <- ggplot(Land.LL.500, aes(y = Larvae_Total_Sum, x = Prop_500m_Y2022_WinterRape))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  xlab(SameYearRapsTitle)+
  ylab(element_blank())+GLM_Line+ggtitle("**")
Larvae.A

Larvae.A.2 <- ggplot(Land.LL.500, aes(y = Larvae_Total_Sum, x = CropDiv21))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  xlab(PreviousYearHet)+
  ylab(element_blank())+GLM_Line+ggtitle("**")
Larvae.A.2

Fig.5.a <- ggarrange(Larvae.A, Larvae.A.2,
                   labels = c("A", "B"), common.legend = TRUE, legend = "bottom",
                   font.label = list(size = 7), nrow = 1, ncol = 2, hjust = -0.9)

Fig.5 <- annotate_figure(Fig.5.a,
                         left = text_grob(LarvaeLab,size = 8, rot = 90))
Fig.5 

#
ggsave(
  "Fig_5.tiff",
  plot = Fig.5,
  device = "tiff",
  scale = 1,
  width = 140,
  height = 80,
  units = c("mm"),
  dpi = 300,
  limitsize = TRUE)
##### SUPPLEMENTARY FIGS: S1 - S3 ######
Landscape1 <- read.csv("Landscape_500m_Radii_Correlations.csv")
Landscape1$Field_ID

Corr500m <- ggpairs(Landscape1[2:7], title = "500 m spatial scale",
                    upper = list(continuous = wrap("cor", size = 3)),
                    lower = list(continuous = wrap("points", size = 0.5)),
                    diag = list(continuous = wrap("densityDiag", size = 0.5)),
                    columnLabels = c("RP 2021", "CH 2021", "RP 2022", "CH 2022", "Change RP", "SNH"))+
  theme(axis.text = element_text(size = 7))

Corr500m

## Correlations with some, update datasheet and repeat
ggsave(
  "Fig_S1.tiff",
  plot = Corr500m,
  device = "tiff",
  scale = 1,
  width = 160,
  height = 130,
  units = c("mm"),
  dpi = 300,
  limitsize = TRUE)

#
Landscape2 <- read.csv("Landscape_750m_Radii_Correlations.csv") 
Landscape2$Field_ID


Corr750m <- ggpairs(Landscape2[2:7], title = "750 m spatial scale",
                    upper = list(continuous = wrap("cor", size = 3)),
                    lower = list(continuous = wrap("points", size = 0.5)),
                    diag = list(continuous = wrap("densityDiag", size = 0.5)),
                    columnLabels = c("RP 2021", "CH 2021", "RP 2022", "CH 2022", "Change RP", "SNH"))+
  theme(axis.text = element_text(size = 7))
Corr750m

## Correlations with some, update datasheet and repeat
ggsave(
  "Fig_S2.tiff",
  plot = Corr750m,
  device = "tiff",
  scale = 1,
  width = 160,
  height = 130,
  units = c("mm"),
  dpi = 300,
  limitsize = TRUE)
#
Landscape3 <- read.csv("Landscape_1000m_Radii_Correlations.csv")
Landscape3$Field_ID

Corr1000m <- ggpairs(Landscape3[2:7], title = "1000 m spatial scale",
                     upper = list(continuous = wrap("cor", size = 3)),
                     lower = list(continuous = wrap("points", size = 0.5)),
                     diag = list(continuous = wrap("densityDiag", size = 0.5)),
                     columnLabels = c("RP 2021", "CH 2021", "RP 2022", "CH 2022", "Change RP", "SNH"))+
  theme(axis.text = element_text(size = 7))
Corr1000m

## Correlations with some, update datasheet and repeat
ggsave(
  "Fig_S3.tiff",
  plot = Corr1000m,
  device = "tiff",
  scale = 1,
  width = 160,
  height = 130,
  units = c("mm"),
  dpi = 300,
  limitsize = TRUE)


##### SUPPLEMENTARY FIGS: S7 & S8 #####
## NatEn
NatEnPlotA1 <- ggplot(Land.NE.500, aes(y = PotentialNatEn, x = CropDiv21))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  ylab(element_blank())+
  xlab(PreviousYearHet)+GLM_Line+ggtitle("500 m")+AssessmentFacetNE
NatEnPlotA1

NatEnPlotA2 <- ggplot(Land.NE.500, aes(y = PotentialNatEn, x = CropDiv22))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  ylab(element_blank())+
  xlab(SamesYearHet)+GLM_Line+ggtitle("500 m")+AssessmentFacetNE
NatEnPlotA2

NatEnPlotA3 <- ggplot(Land.NE.750, aes(y = PotentialNatEn, x = CropDiv22))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  ylab(element_blank())+
  xlab(SamesYearHet)+GLM_Line+GLM_Line+ggtitle("750 m")+AssessmentFacetNE
NatEnPlotA3

NatEnPlotA4 <- ggplot(Land.NE.750, aes(y = PotentialNatEn, x = Prop_SNH_750m))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  ylab(element_blank())+
  xlab(PropSNHTitle)+GLM_Line+ggtitle("750 m")+AssessmentFacetNE
NatEnPlotA4

NatEnPlotA5 <- ggplot(Land.NE.1000, aes(y = PotentialNatEn, x = CropDiv22))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  ylab(element_blank())+
  xlab(SamesYearHet)+GLM_Line+GLM_Line+ggtitle("1000 m")+AssessmentFacetNE
NatEnPlotA5

NatEnPlotA6 <- ggplot(Land.NE.1000, aes(y = PotentialNatEn, x = Prop_SNH_1000m))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  ylab(element_blank())+
  xlab(PropSNHTitle)+GLM_Line+ggtitle("1000 m")+AssessmentFacetNE
NatEnPlotA6

Fig.s7 <- ggarrange(NatEnPlotA1, NatEnPlotA2,
                    NatEnPlotA3, NatEnPlotA4,
                    NatEnPlotA5, NatEnPlotA6,
                     labels = c("A", "B", "C", "D", "E", "F"), common.legend = TRUE, legend = "bottom",
                     font.label = list(size = 7), nrow = 3, ncol = 2, hjust = -0.9)
Fig.s7

Fig.s7b <- annotate_figure(Fig.s7,
                         left = text_grob(NatAbLabel,size = 8, rot = 90))
Fig.s7b 

ggsave(
  "Fig_S7.tiff",
  plot = Fig.s7b,
  device = "tiff",
  scale = 1,
  width = 150,
  height = 200,
  units = c("mm"),
  dpi = 300,
  limitsize = TRUE)

##
NatEnPlotD1 <- ggplot(Land.NE.1000, aes(y = PotentialNatEnDiv, x = Prop_1000m_Y2021_WinterRape))+
  Jitter_Theme+PlotTheme+
  theme(legend.position = "none")+
  ylab("Natural enemy diversity")+
  xlab("Rapeseed proportion: Previous year")+LM_Line+AssessmentFacetNE
NatEnPlotD1

ggsave(
  "Fig_S8.tiff",
  plot = NatEnPlotD1,
  device = "tiff",
  scale = 1,
  width = 80,
  height = 80,
  units = c("mm"),
  dpi = 300,
  limitsize = TRUE)


########################################## STRUCTURAL EQUATION MODEL #########################################
## SEM analysis

#
Land.PP.500$Adj_OSR_2022
Land.PP.500$Adj_OSR_2022 <- plyr::revalue(Land.PP.500$Adj_OSR_2022, c("Y"=1)) # Recode to numeric and pass to SEM
Land.PP.500$Adj_OSR_2022 <- plyr::revalue(Land.PP.500$Adj_OSR_2022, c("N"=0))
Land.PP.500$Adj_OSR_2022 <-  as.numeric(Land.PP.500$Adj_OSR_2022)
Land.PP.500$Assessment_Round <-  as.numeric(Land.PP.500$Assessment_Round)

#
Land.NE.500$Assessment_Round <-  as.numeric(Land.NE.500$Assessment_Round)

#
Land.AD.500$Adj_OSR_202
Land.AD.500$Adj_OSR_2022 <- plyr::revalue(Land.AD.500$Adj_OSR_2022, c("Y"=1)) # Recode to numeric and pass to SEM
Land.AD.500$Adj_OSR_2022 <- plyr::revalue(Land.AD.500$Adj_OSR_2022, c("N"=0))
Land.AD.500$Adj_OSR_2022 <-  as.numeric(Land.AD.500$Adj_OSR_2022)
Land.AD.500$Assessment_Round <-  as.numeric(Land.AD.500$Assessment_Round)

Land.LL.500$Adj_OSR_2022 <- plyr::revalue(Land.LL.500$Adj_OSR_2022, c("Y"=1)) # Recode to numeric and pass to SEM
Land.LL.500$Adj_OSR_2022 <- plyr::revalue(Land.LL.500$Adj_OSR_2022, c("N"=0))
Land.LL.500$Adj_OSR_2022 <-  as.numeric(Land.LL.500$Adj_OSR_2022)


Land.LL.750$Adj_OSR_2022 <- plyr::revalue(Land.LL.750$Adj_OSR_2022, c("Y"=1)) # Recode to numeric and pass to SEM
Land.LL.750$Adj_OSR_2022 <- plyr::revalue(Land.LL.750$Adj_OSR_2022, c("N"=0))
Land.LL.750$Adj_OSR_2022 <-  as.numeric(Land.LL.750$Adj_OSR_2022)

Land.LL.1000$Adj_OSR_2022 <- plyr::revalue(Land.LL.1000$Adj_OSR_2022, c("Y"=1)) # Recode to numeric and pass to SEM
Land.LL.1000$Adj_OSR_2022 <- plyr::revalue(Land.LL.1000$Adj_OSR_2022, c("N"=0))
Land.LL.1000$Adj_OSR_2022 <-  as.numeric(Land.LL.1000$Adj_OSR_2022)

### NOTE QUIRK WITH psem - sometimes the model needs to be built "from scratch" by adding in the individual models in a step-wise manner
## Otherwise an "NA in dataset" error occurs

# Join the datasets into SEM datasets
SEM.Data.500.Damage <- full_join(x = Land.PP.500, y = Land.AD.500)
Land.PP.Sub <- subset(Land.PP.500, Assessment_Round %in% c ("2"))
SEM.Data.500.Larvae <- full_join(x = Land.PP.Sub, y = Land.LL.500)

SEM.Data.500.NE.Sub <- subset(Land.NE.500, Assessment_Round %in% c ("1", "2"))
SEM.Data.500.NE <- full_join(x = Land.PP.500, y = SEM.Data.500.NE.Sub)

which(is.na(SEM.Data.500.Damage), arr.ind=TRUE)
which(is.na(SEM.Data.500.Larvae), arr.ind=TRUE)
which(is.na(Land.PP.500), arr.ind=TRUE)
which(is.na(SEM.Data.500.NE), arr.ind=TRUE)

# SEM code
SEM.Land.NE <- psem(glmmTMB(Pressure_Cummulative ~ 
                              Prop_500m_Y2022_WinterRape+
                              CropDiv21+CropDiv22+
                              Adj_OSR_2022+Assessment_Round+(1|Field),
                            data=Land.PP.500),
                    
                    lmer(Percent_Damage_Logit ~ Pressure_Cummulative+
                           Prop_500m_Y2022_WinterRape+
                           CropDiv21+CropDiv22+
                           Adj_OSR_2022+Assessment_Round+(1|Field),
                         data = SEM.Data.500.Damage,
                         na.action = na.omit), 
                    
                    glmmTMB(Larvae_Total_Sum ~ Pressure_Cummulative+
                              Prop_500m_Y2022_WinterRape+
                              CropDiv21+CropDiv22+
                              Adj_OSR_2022+(1|Field), 
                            data = SEM.Data.500.Larvae),
                    
                    glmmTMB(PotentialNatEn ~ 
                              Prop_500m_Y2022_WinterRape+
                              CropDiv21+CropDiv22+Assessment_Round+(1|Field),
                            data=SEM.Data.500.NE)
                    
) # Note - model might bring up "NA in dataset" error, but datasets have been checked and no NA are present in the data

summary(SEM.Land.NE, conserve = TRUE)
coefs(SEM.Land.NE)
fisherC(SEM.Land.NE, conserve = TRUE)
AIC(SEM.Land.NE)
#