---
title: "Hyjazie and Forrest, 2023: Data Organization and Code"
date: "April 20, 2023"
output:
  pdf_document: default
  word_document: default
---

```{r, echo = FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE, error = FALSE, tidy.opts=list(width.cutoff=80), tidy=TRUE)

```

# Set-Up

## Load Packages

```{r Packages, include = FALSE}

# Packages for Data Organization
library(readxl)
library(tidyverse)
library(tidyr)
library(dplyr)
library(plyr)
library(purrr)

# Packages for Models
library(lme4)
library(glmmTMB) 
library(MASS)
library(lmtest)
library(lmerTest)
library(scales)
library(RVAideMemoire)
library(DHARMa) 

# Packages for Plots
library(ggplot2)
library(gridExtra)
library(grid)
library(lattice)
library(car)
library(RColorBrewer)
library(cowplot)

```

# Main Manuscript

## Data Organization: Year 1

### 1. Flower-Visitor Abundance

#### 1.1 *Osmia* spp. Abundance

```{r Data Osmia Abundance Year 1, results = "hide"}

# Loading Excel file
Abundance_Osmia_Y1 <- read_excel("Osmia_Abundance_1.xlsx")

# Merging variable names (Orchard (the logistical orchard name according to orchard pairings) and Site ID) to create unique identifiers for every site in new column "Site"
Abundance_Osmia_1 <-  Abundance_Osmia_Y1 %>% unite('Site', Orchard, `Site ID`, remove = FALSE)

# Verifying that Site = 24 and Orchard = 16
Abundance_Osmia_1_Site_Count <- n_distinct(Abundance_Osmia_1$Site)
Abundance_Osmia_1_Site_Count

Abundance_Osmia_1_Orchard_Count <- n_distinct(Abundance_Osmia_1$`Orchard ID`)
Abundance_Osmia_1_Orchard_Count

# Summarizing *Osmia* counts per transect walk according to plot, Orchard ID (the proper orchard name), Round, and Site
Abundance_Osmia_1 <- Abundance_Osmia_1 %>% group_by(Plot, `Orchard ID`, Round, Site) %>%
  summarise_at(vars(Count), list(Count=sum))

```

#### 1.2 All Flower Visitors

All flower visitors during year 1. 

```{r Data Total Abundance, results = "hide"}

# Loading Excel file
Abundance_Total_Y1 <- read_excel("Total_Abundance_1.xlsx")

# Merging variable names (Orchard and Site ID) to create unique identifiers for every site in new column "Site"
Abundance_Total_1 <- Abundance_Total_Y1 %>% unite('Site', Orchard, `Site ID`, remove = FALSE)

# Summarizing counts per transect walk according to plot, Orchard ID (the proper orchard name), Round, and Site
Abundance_Total_1 <- Abundance_Total_1 %>% group_by(Plot, `Orchard ID`, Round, Site) %>%
  summarise_at(vars(Count), list(Count=sum))

```

#### 1.3 Small Flower Visitors

Year 1; excludes honeybees and bumblebees.

```{r Selective Data Abundance, results = "hide"}

# Loading Excel file
Abundance_Select_Y1 <- read_excel("Total_Abundance_1.xlsx")

# Merging variable names (Orchard and Site ID) to create unique identifiers for every site in new column "Site"
Abundance_Select_1 <- Abundance_Select_Y1 %>% unite('Site', Orchard, `Site ID`, remove = FALSE)

#Creating small insect subset
Abundance_Select_1 <- Abundance_Select_1 %>% filter(Taxa != "Apis", Taxa != "Bombus", Taxa != "Bombus Queen")

# Summarizing counts per transect walk according to plot, Orchard ID (the proper orchard name), Round, and Site
Abundance_Select_1 <- Abundance_Select_1 %>% group_by(Plot, `Orchard ID`, Round, Site) %>%
  summarise_at(vars(Count), list(Count=sum))

```


### 2. Apple Yield and Quality

```{r Data Apple, results = "hide"}

# Loading Excel File
Apple_1 <- read_excel("Apple_1.xlsx")

# Merging variable names (Orchard and Site ID) to create unique identifiers for every site in new column "Site"
Apple_1 <-  Apple_1 %>% unite('Site', Orchard, `Site ID`, remove = FALSE)

# Re-labeling columns
Apple_1$Mean_Circum <- as.numeric(Apple_1$`Mean Non-Zero Circumference (cm)`)

# Creating columns for total fruit and bud numbers per tree (N + S branches)
Apple_1$Buds <- Apple_1$`N Bud #` + Apple_1$`S Bud #`
Apple_1$Pre_Fruit <- Apple_1$`N Pre-Thinning Fruit #` + Apple_1$`S Pre-Thinning Fruit #`
Apple_1$Post_Fruit <- Apple_1$`N Post-Thinning Fruit #` + Apple_1$`S Post-Thinning Fruit #`

# Replacing NA values for post-fruit apple with 0
Apple_1 <- Apple_1 %>%
  mutate(Post_Fruit=replace_na(Post_Fruit, 0))

# Calculating fruit set pre- and post-thinning
Apple_1$Fruit_Set_Pre <- Apple_1$Pre_Fruit/Apple_1$Buds
Apple_1$Fruit_Set_Post <- Apple_1$Post_Fruit/Apple_1$Buds

# Verifying that Site = 24 and Orchard = 16
Apple_1_Site_Count <- n_distinct(Apple_1$Site)
Apple_1_Site_Count

Apple_1_Orchard_Count <- n_distinct(Apple_1$`Orchard ID`)
Apple_1_Orchard_Count

```

## Results: Year 1

### 1. Flower-Visitor Abundance

```{r  Abundance Mean (All Taxa)}

# Sum of counts by Taxa per Round
Total_Abundance_1_Sum <- Abundance_Total_Y1 %>%
  group_by(Taxa, Round) %>%
  dplyr::summarize_at(vars(Count), list(Count=sum))
Total_Abundance_1_Sum

```

#### 1.1 *Osmia* spp. Abundance

```{r Osmia Abundance Mean and se}

#Means and SEs by plot (treatment) across sites and rounds
Abundance_Osmia_Means_1 <- Abundance_Osmia_1 %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Count), list(Mean_Count = mean, SE_Count = se))
Abundance_Osmia_Means_1

# Overall Percent Difference between treatments for *Osmia* based on mean counts for each sampling round:
Abundance_Osmia_Means_1 %>%
  summarise(percent_diff = ((Mean_Count[2] - Mean_Count[1])/mean(Mean_Count)))

# *Osmia* summed across Rounds
Abundance_Osmia_Values_1 <- Abundance_Osmia_1 %>%
  group_by(Plot, Site) %>%
  dplyr::summarize_at(vars(Count), list(Total_Osmia_Count = sum))
Abundance_Osmia_Values_1

#Means by plot (treatment) of abundance values summed across rounds
Abundance_Osmia_Totals_1 <- Abundance_Osmia_Values_1 %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Total_Osmia_Count), list(Mean_Total_Osmia = mean, SE_Count = se))
Abundance_Osmia_Totals_1

# Overall Percent Difference between treatments for *Osmia*, based on site totals:
Abundance_Osmia_Totals_1 %>%
  summarise(percent_diff = ((Mean_Total_Osmia[2]-Mean_Total_Osmia[1])/mean(Mean_Total_Osmia)))

# Mean and SE of Percent Difference between Site pairs:
Abundance_Osmia_Percent_Diff_by_Site <- pivot_wider(Abundance_Osmia_Values_1, names_from = Plot, values_from = Total_Osmia_Count)
Abundance_Osmia_Percent_Diff_by_Site$Diff <- Abundance_Osmia_Percent_Diff_by_Site$Treatment - Abundance_Osmia_Percent_Diff_by_Site$Control
Abundance_Osmia_Percent_Diff_by_Site$Mean_Abd <- (Abundance_Osmia_Percent_Diff_by_Site$Treatment+Abundance_Osmia_Percent_Diff_by_Site$Control)/2 
Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff <- 100*Abundance_Osmia_Percent_Diff_by_Site$Diff/Abundance_Osmia_Percent_Diff_by_Site$Mean_Abd
Abundance_Osmia_Percent_Diff_by_Site
mean(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff); sd(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff)/sqrt(length(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff))

```

##### 1.1.1 Models

```{r Osmia Abundance Model (Year 1)}

# Is there a difference between control and treatment plots for Osmia abundance, and does the treatment effect ("Plot") vary depending on sampling round ("Round")?

#First look at data distribution:
hist(Abundance_Osmia_1$Count)

#Initial model (Poisson):
Osmia_Abund_GLMM_1 <- glmer(Count ~ Plot*Round + (1|`Orchard ID`/Site), 
                    family=poisson(link = "log"),
                    data=Abundance_Osmia_1)

# Model diagnostics:
overdisp.glmer(Osmia_Abund_GLMM_1)
testDispersion(Osmia_Abund_GLMM_1)
simulationOutput <- simulateResiduals(fittedModel = Osmia_Abund_GLMM_1, plot = F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
#Some problems here, so try zero-inflated model fit instead.

#Zero-inflated Poisson using glmmTMB:
Osmia_Abund_ZI_1 <- glmmTMB(Count ~ Plot*Round + (1|`Orchard ID`/Site), 
                    ziformula = ~1 , family = "poisson",
                    data=Abundance_Osmia_1)

#Diagnostics
testDispersion(Osmia_Abund_ZI_1)
simulationOutput <- simulateResiduals(fittedModel = Osmia_Abund_ZI_1, plot = F)
plot(simulationOutput)

#Results of selected model:
summary(Osmia_Abund_ZI_1)

# Without Round as an Interaction (for comparison in anova)
Osmia_Abund_ZI_1a <- glmmTMB(Count ~ Plot + Round + (1|`Orchard ID`/Site), 
                    ziformula = ~1 , family = "poisson",
                    data=Abundance_Osmia_1)
anova(Osmia_Abund_ZI_1, Osmia_Abund_ZI_1a)
#Interaction NS.

# Test significance of Round term:
Osmia_Abund_ZI_1b <- glmmTMB(Count ~ Plot + (1|`Orchard ID`/Site), 
                    ziformula = ~1 , family = "poisson",
                    data=Abundance_Osmia_1)
anova(Osmia_Abund_ZI_1b, Osmia_Abund_ZI_1a)

```


#### 1.2 All Flower Visitors

```{r Pollinator Abundance Mean and SE}

# Sums by plot (treatment) and site for all flower visitors
Total_Abundance_by_Site <- Abundance_Total_1 %>%
  group_by(Site, Plot) %>%
  dplyr::summarize_at(vars(Count), list(Total_Abd=sum))
Total_Abundance_by_Site

# Means by plot (treatment) for all flower visitors
Total_Abundance_by_Plot <- Total_Abundance_by_Site %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Total_Abd), list(Mean_Totals=mean, SE_Totals = se))
Total_Abundance_by_Plot

# Overall Percent Difference in flower visitor numbers
Total_Abundance_by_Plot %>%
  summarise(percent_diff = ((Mean_Totals[2]-Mean_Totals[1])/mean(Mean_Totals)))

# Mean Percent Difference across Sites
Abundance_Total_Values_Percent_Diff_by_Site <- pivot_wider(Total_Abundance_by_Site, names_from = Plot, values_from = Total_Abd)
Abundance_Total_Values_Percent_Diff_by_Site$Diff <- Abundance_Total_Values_Percent_Diff_by_Site$Treatment - Abundance_Total_Values_Percent_Diff_by_Site$Control
Abundance_Total_Values_Percent_Diff_by_Site$Mean_Abd <- (Abundance_Total_Values_Percent_Diff_by_Site$Treatment+Abundance_Total_Values_Percent_Diff_by_Site$Control)/2 
Abundance_Total_Values_Percent_Diff_by_Site$Percent_Diff <- 100*Abundance_Total_Values_Percent_Diff_by_Site$Diff/Abundance_Total_Values_Percent_Diff_by_Site$Mean_Abd
Abundance_Total_Values_Percent_Diff_by_Site
mean(Abundance_Total_Values_Percent_Diff_by_Site$Percent_Diff)

```

##### 1.2.1 Models

```{r Pollinator Abundance Model (Year 1)}

# Is there a different between control and treatment plots for total flower-visitor abundance, and does the treatment effect ("Plot") vary depending on sampling round ("Round")?

#First look at data distribution:
hist(Abundance_Total_1$Count)

#Initial model (Poisson):
Total_Abund_GLMM_1 <- glmer(Count ~ Plot * Round + (1|`Orchard ID`/Site), 
                    family=poisson(link = "log"),
                    data=Abundance_Total_1)

# Model diagnostics:
overdisp.glmer(Total_Abund_GLMM_1)
testDispersion(Total_Abund_GLMM_1)
simulationOutput <- simulateResiduals(fittedModel = Total_Abund_GLMM_1, plot = F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
#Serious overdispersion problem here, but not zero-inflation, so try negative binomial model fit instead.

#Negative binomial using glmmTMB:
Total_Abund_NB_1 <- glmmTMB(Count ~ Plot*Round + (1|`Orchard ID`/Site), 
                    family = nbinom1,
                    data=Abundance_Total_1)

#Diagnostics
testDispersion(Total_Abund_NB_1)
simulationOutput <- simulateResiduals(fittedModel = Total_Abund_NB_1, plot = F)
plot(simulationOutput)

#Results of selected model:
summary(Total_Abund_NB_1)

# Without Round as an Interaction (for comparison in anova)
Total_Abund_NB_1a <- glmmTMB(Count ~ Plot + Round + (1|`Orchard ID`/Site), 
                    family = nbinom1,
                    data=Abundance_Total_1)
#Test significance of interaction term:
anova(Total_Abund_NB_1, Total_Abund_NB_1a)

# Test significance of Round term:
Total_Abund_NB_1b <- glmmTMB(Count ~ Plot + (1|`Orchard ID`/Site), 
                    family = nbinom1,
                    data=Abundance_Total_1)
anova(Total_Abund_NB_1b, Total_Abund_NB_1a)

```

#### 1.3 Small Flower Visitors

```{r Small Flower Abundance Mean and SE}

# Sums by plot (treatment) and site for all flower visitors EXCLUDING APIS AND BOMBUS
Select_Abundance_by_Site <- Abundance_Select_1 %>%
  group_by(Site, Plot) %>%
  dplyr::summarize_at(vars(Count), list(Total_Abd=sum))
Select_Abundance_by_Site

# Means by plot (treatment) 
Select_Abundance_by_Plot <- Select_Abundance_by_Site %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Total_Abd), list(Mean_Totals=mean, SE_Totals = se))
Select_Abundance_by_Plot

# Overall Percent Difference 
Select_Abundance_by_Plot %>%
  summarise(percent_diff = ((Mean_Totals[2]-Mean_Totals[1])/mean(Mean_Totals)))

# Mean Percent Difference across Sites
Abundance_Select_Values_Percent_Diff_by_Site <- pivot_wider(Select_Abundance_by_Site, names_from = Plot, values_from = Total_Abd)
Abundance_Select_Values_Percent_Diff_by_Site$Diff <- Abundance_Select_Values_Percent_Diff_by_Site$Treatment - Abundance_Select_Values_Percent_Diff_by_Site$Control
Abundance_Select_Values_Percent_Diff_by_Site$Mean_Abd <- (Abundance_Select_Values_Percent_Diff_by_Site$Treatment+Abundance_Select_Values_Percent_Diff_by_Site$Control)/2 
Abundance_Select_Values_Percent_Diff_by_Site$Percent_Diff <- 100*Abundance_Select_Values_Percent_Diff_by_Site$Diff/Abundance_Select_Values_Percent_Diff_by_Site$Mean_Abd
Abundance_Select_Values_Percent_Diff_by_Site
mean(Abundance_Select_Values_Percent_Diff_by_Site$Percent_Diff)

```

##### 1.3.1 Models

```{r Small Flower-Visitor Abundance Model (Year 1)}

# Is there a different between control and treatment plots for small flower-visitor abundance, and does the treatment effect ("Plot") vary depending on sampling round ("Round")?

#First look at data distribution:
hist(Abundance_Select_1$Count)

#Initial model (Poisson):
Select_Abund_GLMM_1 <- glmer(Count ~ Plot*Round + (1|`Orchard ID`/Site), 
                    family=poisson(link = "log"),
                    data=Abundance_Select_1)

# Model diagnostics:
overdisp.glmer(Select_Abund_GLMM_1)
testDispersion(Select_Abund_GLMM_1)
simulationOutput <- simulateResiduals(fittedModel = Select_Abund_GLMM_1, plot = F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
#Some problems here, so try zero-inflated model fit instead.

#Zero-inflated Poisson using glmmTMB:
Select_Abund_ZI_1 <- glmmTMB(Count ~ Plot*Round + (1|`Orchard ID`/Site), 
                    ziformula = ~1 , family = "poisson",
                    data=Abundance_Select_1)

#Diagnostics
testDispersion(Select_Abund_ZI_1)
simulationOutput <- simulateResiduals(fittedModel = Select_Abund_ZI_1, plot = F)
plot(simulationOutput)
#Still some issues, but not as bad.

#Negative binomial using glmmTMB:
Select_Abund_NB_1 <- glmmTMB(Count ~ Plot*Round + (1|`Orchard ID`/Site), 
                    family = nbinom1,
                    data=Abundance_Select_1)

#Diagnostics
testDispersion(Select_Abund_NB_1)
simulationOutput <- simulateResiduals(fittedModel = Select_Abund_NB_1, plot = F)
plot(simulationOutput)
#Much better, though still some overdispersion.

#Results of selected model:
summary(Select_Abund_NB_1)

#Negative binomial without interaction, to test significance:
Select_Abund_NB_1a <- glmmTMB(Count ~ Plot + Round + (1|`Orchard ID`/Site), 
                    family = nbinom1,
                    data=Abundance_Select_1)
anova(Select_Abund_NB_1a, Select_Abund_NB_1)
#Interaction NS.

# Test significance of Round term:
Select_Abund_NB_1b <- glmmTMB(Count ~ Plot + (1|`Orchard ID`/Site), 
                     family = nbinom1,
                    data=Abundance_Select_1)
anova(Select_Abund_NB_1b, Select_Abund_NB_1a)
#Round NS.

```
### 2. Apple Yield and Quality

#### 2.1 Pre- and Post-Thinning Fruit Set

```{r Apple Pre_Thinning Values Mean and se}

# Mean and SE of fruit-set pre-thinning values per plot
Apple_Values_Pre <- Apple_1 %>%
  group_by(Plot) %>%
  dplyr::summarize(Mean_Fruit_Set_Pre=mean(Fruit_Set_Pre),
                   SE_Fruit_Set_Pre=se(Fruit_Set_Pre))
Apple_Values_Pre

# Mean and SE of fruit set post-thinning per plot
Apple_Post_Values <- Apple_1 %>%
  group_by(Plot) %>%
  dplyr::summarize(Mean_Fruit_Set_Post = mean(Fruit_Set_Post), 
                   SE_Fruit_Set_Post = se(Fruit_Set_Post))
Apple_Post_Values
```

##### 2.1.1 Pre-Thinning Model

```{r Apple Yield Pre-Thinning Year 1 Model}

hist(Apple_1$Fruit_Set_Pre)
hist(Apple_1$Pre_Fruit)

# Binomial GLMER With Site as a Random Effect
Apple_Pre_Thin_GLM <- glmer(cbind(Pre_Fruit,Buds) ~ Plot + (1|`Orchard ID`/Site), 
                    family = binomial(link = "logit"),
                    data=Apple_1)
summary(Apple_Pre_Thin_GLM)

# Model diagnostics:
overdisp.glmer(Apple_Pre_Thin_GLM)
testDispersion(Apple_Pre_Thin_GLM)
simulationOutput <- simulateResiduals(fittedModel = Apple_Pre_Thin_GLM, plot = T)
#Fit isn't great.

#Ordered beta regression with Site random:
Apple_Pre_Thin_GLM <- glmmTMB(Fruit_Set_Pre ~ Plot + (1|`Orchard ID`/Site), 
                              dispformula = ~1,                  
                              family = ordbeta,
                    data=Apple_1)
summary(Apple_Pre_Thin_GLM)

# Model diagnostics:
testDispersion(Apple_Pre_Thin_GLM)
simulationOutput <- simulateResiduals(fittedModel = Apple_Pre_Thin_GLM, plot = F)
plot(simulationOutput)
#Also not great, but same conclusion.

```

##### 2.1.2 Post-Thinning Model

```{r Apple Yield Post-Thinning Model}

# With Site as a Random Effect

hist(Apple_1$Fruit_Set_Post)

# Binomial GLMER With Site as a Random Effect

Apple_Post_Thin_GLM <- glmer(Post_Fruit/Buds ~ Plot + (1|`Orchard ID`/Site), 
                    family = binomial(link = "logit"),
                    data=Apple_1)
summary(Apple_Post_Thin_GLM)

# Model diagnostics:
overdisp.glmer(Apple_Post_Thin_GLM)
testDispersion(Apple_Post_Thin_GLM)
simulationOutput <- simulateResiduals(fittedModel = Apple_Post_Thin_GLM, plot = T)
#Fit is terrible.

#Ordered beta regression with Site random:
Apple_Post_Thin_GLM <- glmmTMB(Fruit_Set_Post ~ Plot + (1|`Orchard ID`/Site), 
                              dispformula = ~1,                  
                              family = ordbeta,
                    data=Apple_1)
summary(Apple_Post_Thin_GLM)

# Model diagnostics:
testDispersion(Apple_Post_Thin_GLM)
simulationOutput <- simulateResiduals(fittedModel = Apple_Post_Thin_GLM, plot = F)
plot(simulationOutput)
#Not great, but much better.

```

#### 2.3 Apple Quality

```{r Apple Quality Mean and se}

#First summarizing by treatment within sites:
Apple_1_site_summary <- Apple_1 %>%
  group_by(Plot, Site) %>%
  dplyr::summarize(circum = mean(Mean_Circum, na.rm = T),
                   count = n())
Apple_1_site_summary

#Then summarizing by treatment across sites
Apple_1_site_summary <- Apple_1_site_summary %>%
  group_by(Plot) %>%
  dplyr::summarize(mean = mean(circum, na.rm = T),
                   SE = se(circum),
                   count = n())
Apple_1_site_summary

#How many trees with non-zero values in each treatment?
dim(subset(subset(Apple_1, `Fruits Measured`>0), Plot == "Treatment"))
dim(subset(subset(Apple_1, `Fruits Measured`>0), Plot == "Control"))

```

##### 2.3.1 Models

```{r Apple Quality Model}

#Distribution
hist(Apple_1$Mean_Circum) 

# With Site as a Random Effect
Apple_Quality_LMM <- lmer(Mean_Circum ~ Plot + (1|`Orchard ID`/Site), 
                    data=Apple_1)
summary(Apple_Quality_LMM)

# Assumptions
overdisp.glmer(Apple_Quality_LMM)
testDispersion(Apple_Quality_LMM)
simulationOutput <- simulateResiduals(fittedModel = Apple_Quality_LMM, plot = T)
#All looks good.
```

## Data Organization: Year 2

### 1. Flower-Visitor Abundance

#### 1.1 *Osmia* spp. Abundance

```{r Data Osmia Abundance Year 2, results = "hide"}

# Loading Excel file (Year 2)
Abundance_Total_Y2 <- read_excel("Total_Abundance_Year_2.xlsx")

# Merging variable names to create unique identifiers for every site
Abundance_Total_2 <- Abundance_Total_Y2 %>% unite('Site', Orchard, `Site ID`, remove = FALSE)

#Adding Year column
Abundance_Total_2 <- Abundance_Total_2 %>% add_column(Year = 2022)

# Only look at Osmia
Abundance_Osmia_2 <- Abundance_Total_2 %>% filter(Taxa == "Osmia")

# Summarizing *Osmia* counts per transect walk:
Abundance_Osmia_2 <- Abundance_Osmia_2 %>% group_by(Plot, Site, `Orchard ID`, Year) %>%
  summarise_at(vars(Count), list(Count=sum))

# Extract observations from apple bloom of Year 1 and reorganize to match Year 2:
Abundance_Osmia_Y1_AB <- Abundance_Osmia_1 %>% filter(Round == "Apple Bloom") %>% add_column(Year = 2021)
Abundance_Osmia_Y1_AB <- Abundance_Osmia_Y1_AB[,c(1,4,2,6,5)]

# Combine both years (apple bloom period only):
Abundance_Osmia_AB <- rbind(Abundance_Osmia_Y1_AB, Abundance_Osmia_2)

```

#### 1.2 All Flower Visitors 

```{r Data Pollinator Abundance Year 2, results = "hide"}

# Loading Excel file (Year 2)
Abundance_Total_2 <- read_excel("Total_Abundance_Year_2.xlsx")

# Merging variable names to create unique identifiers for every site
Abundance_Total_2 <-  Abundance_Total_2 %>% unite('Site', Orchard, `Site ID`, remove = FALSE)
  
# Adding Year column
Abundance_Total_2 <- Abundance_Total_2 %>% add_column(Year = 2022)

# Summarizing Flower visitor counts counts per transect walk
Abundance_Total_2 %>% group_by(Taxa) %>%
  summarise_at(vars(Count), list(Count=sum))

# Summarizing pollinator counts per transect walk and adding Year column:
Abundance_Total_2 <- Abundance_Total_2 %>% group_by(Plot, Site, `Orchard ID`, Year) %>%
  summarise_at(vars(Count), list(Count=sum))

# Extract observations from apple bloom of Year 1 and reorganize to match Year 2:
Abundance_Total_Y1_AB <- Abundance_Total_1 %>% filter(Round == "Apple Bloom") %>% add_column(Year = 2021)
Abundance_Total_Y1_AB <- Abundance_Total_Y1_AB[,c(1,4,2,6,5)]

# Combine both years (apple bloom period only):
Abundance_Total_AB <- rbind(Abundance_Total_Y1_AB, Abundance_Total_2)

```

#### 1.3 Small Flower Visitors

```{r Data Small Flower Abundance Year 2, results = "hide"}

# Loading Excel file (Year 2)
Abundance_Select_Y2 <- read_excel("Total_Abundance_Year_2.xlsx")

# Merging variable names to create unique identifiers for every site
Abundance_Select_2 <- Abundance_Select_Y2 %>% unite('Site', Orchard, `Site ID`, remove = FALSE)

#Adding Year column
Abundance_Select_2 <- Abundance_Select_2 %>% add_column(Year = 2022)

# Exclude Apis and Bombus
Abundance_Select_2 <- Abundance_Select_2 %>% filter(Taxa != "Apis", Taxa != "Bombus", Taxa != "Bombus Queen")

# Summarizing small insect counts per transect walk:
Abundance_Select_2 <- Abundance_Select_2 %>% group_by(Plot, Site, `Orchard ID`, Year) %>%
  summarise_at(vars(Count), list(Count=sum))

# Extract observations from apple bloom of Year 1 and reorganize to match Year 2:
Abundance_Select_Y1_AB <- Abundance_Select_1 %>% filter(Round == "Apple Bloom") %>% add_column(Year = 2021)
Abundance_Select_Y1_AB <- Abundance_Select_Y1_AB[,c(1,4,2,6,5)]

# Combine both years (apple bloom period only):
Abundance_Select_AB <- rbind(Abundance_Select_Y1_AB, Abundance_Select_2)

```

#### 1.4 All Years and Rounds (for Plotting)

##### 1.4.1 *Osmia* spp.

```{r Data Osmia Abundance All Years & Rounds, results = "hide"}

# Add columns for year and round, and reorganize for years to match:
Abundance_Osmia_1 <- Abundance_Osmia_1 %>% add_column(Year = 2021)
Abundance_Osmia_1 <- Abundance_Osmia_1 %>%
  relocate(Site, .after = Plot) 
Abundance_Osmia_2 <- Abundance_Osmia_2 %>% add_column(Round = "Apple Bloom")
Abundance_Osmia_2 <- Abundance_Osmia_2 %>%
  relocate(Round, .after = `Orchard ID`) 

# Bind Osmia Year 1 and Year 2 Data
Abundance_Osmia_All <- rbind(Abundance_Osmia_1, Abundance_Osmia_2)
Abundance_Osmia_All <- Abundance_Osmia_All %>% add_column(Taxa = "Osmia")

# Rename columns
Abundance_Osmia_All$Round <- factor(Abundance_Osmia_All$Round, levels=c("Apple Bloom", "Round 2", "Round 3"), labels=c("Apple bloom", "Pre-thinning", "Post-thinning"))
Abundance_Osmia_All$Round_Year <- paste(Abundance_Osmia_All$Round, Abundance_Osmia_All$Year, sep = " ")

```

##### 1.4.2 All Flower Visitors

```{r Data Pollinator Abundance All Years, results = "hide"}

# Add columns for year and round, and reorganize for years to match:
Abundance_Total_1 <- Abundance_Total_1 %>% add_column(Year = 2021)
Abundance_Total_1 <- Abundance_Total_1 %>%
  relocate(Site, .after = Plot) 
Abundance_Total_2 <- Abundance_Total_2 %>% add_column(Round = "Apple Bloom")
Abundance_Total_2 <- Abundance_Total_2 %>%
  relocate(Round, .after = `Orchard ID`) 

# Bind Year 1 and Year 2 Data
Abundance_Total_All <- rbind(Abundance_Total_1, Abundance_Total_2)
Abundance_Total_All <- Abundance_Total_All %>% add_column(Taxa = "All")

# Rename columns
Abundance_Total_All$Round <- factor(Abundance_Total_All$Round, levels=c("Apple Bloom", "Round 2", "Round 3"), labels=c("Apple bloom", "Pre-thinning", "Post-thinning"))
Abundance_Total_All$Round_Year <- paste(Abundance_Total_All$Round, Abundance_Total_All$Year, sep = " ")

```

### 2. Apple Yield and Quality

#### 2.1 Pre-Thinning

```{r Apple Year 2 Data, results = "hide"}

# Loading Excel File
Apple_2 <- read_excel("Apple_Year_2.xlsx", na = "N/A")

# Merging variable names (Orchard (the logistical orchard name according to orchard pairings) and Site ID) to create unique identifiers for every site in new column "Site"
Apple_2 <-  Apple_2 %>% unite('Site', `Orchard`, `Site ID`, remove = FALSE)

# Creating columns for total fruit and bud numbers per tree (N + S branches)
Apple_2$Buds <- Apple_2$`N Bud #` + Apple_2$`S Bud #`
Apple_2$Pre_Fruit <- Apple_2$`N Pre-Thinning Fruit #` + Apple_2$`S Pre-Thinning Fruit #`

# Calculating fruit set pre-thinning
Apple_2$Fruit_Set_Pre <- Apple_2$Pre_Fruit/Apple_2$Buds

```

#### 2.2 All Years

```{r Data Apple All Years, results = "hide"}

Apple_1 <- Apple_1 %>% add_column(Year = 2021)

Apple_2 <- Apple_2 %>% add_column(Year = 2022)

Apple_All <- bind_rows(Apple_1, Apple_2)
```

## Results: Year 2

### 1. Flower-Visitor Abundance

#### 1.1 *Osmia*spp. Abundance

```{r 1.1 Osmia Means and se}

# Calculating means for Osmia values for both years (apple bloom period only):
Abundance_Osmia_Summary <- Abundance_Osmia_AB %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Count), list(Mean_Count = mean, SE_Count = se))
Abundance_Osmia_Summary

# Percent Difference between treatments
Abundance_Osmia_Summary %>%
  summarise(percent_diff = (100*(Mean_Count[2]-Mean_Count[1])/mean(Mean_Count)))

# Mean and SE of Percent Difference between Site pairs, both Years:
# Note that this is unreliable because of a lot of 0 values generating NAs.
Abundance_Osmia_Percent_Diff_by_Site <- pivot_wider(Abundance_Osmia_AB[,c(-3)], id_cols=c(Site, Year), names_from = Plot, values_from = Count)
Abundance_Osmia_Percent_Diff_by_Site$Diff <- Abundance_Osmia_Percent_Diff_by_Site$Treatment - Abundance_Osmia_Percent_Diff_by_Site$Control
Abundance_Osmia_Percent_Diff_by_Site$Mean_Abd <- (Abundance_Osmia_Percent_Diff_by_Site$Treatment+Abundance_Osmia_Percent_Diff_by_Site$Control)/2 
Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff <- 100*Abundance_Osmia_Percent_Diff_by_Site$Diff/Abundance_Osmia_Percent_Diff_by_Site$Mean_Abd
Abundance_Osmia_Percent_Diff_by_Site
mean(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff, na.rm = T); sd(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff, na.rm = T)/sqrt(length(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff))

# Calculating means for Osmia values for year 2
Abundance_Osmia_Y2_Summary <- Abundance_Osmia_2 %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Count), list(Mean_Count=mean, SE_Count = se))
Abundance_Osmia_Y2_Summary

# Overall Percent Difference between treatments Year 2
Abundance_Osmia_Y2_Summary %>%
  summarise(percent_diff = (100*(Mean_Count[2]-Mean_Count[1])/mean(Mean_Count)))

# Mean and SE of Percent Difference between Site pairs, Year 2:
# Note that this is unreliable because of a lot of 0 values generating NAs.
Abundance_Osmia_Percent_Diff_by_Site <- pivot_wider(Abundance_Osmia_2[,c(-3)], names_from = Plot, values_from = Count)
Abundance_Osmia_Percent_Diff_by_Site$Diff <- Abundance_Osmia_Percent_Diff_by_Site$Treatment - Abundance_Osmia_Percent_Diff_by_Site$Control
Abundance_Osmia_Percent_Diff_by_Site$Mean_Abd <- (Abundance_Osmia_Percent_Diff_by_Site$Treatment+Abundance_Osmia_Percent_Diff_by_Site$Control)/2 
Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff <- 100*Abundance_Osmia_Percent_Diff_by_Site$Diff/Abundance_Osmia_Percent_Diff_by_Site$Mean_Abd
Abundance_Osmia_Percent_Diff_by_Site
mean(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff, na.rm = T); sd(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff, na.rm = T)/sqrt(length(Abundance_Osmia_Percent_Diff_by_Site$Percent_Diff))

```

##### 1.1.1 Models

```{r Abundance_Osmia with Year Model}

# Only Year 2

hist(Abundance_Osmia_2$Count)

Osmia_Abund_GLMM_Y2 <- glmer(Count ~ Plot + (1|`Orchard ID`/Site), 
                    family=poisson(link = "log"),
                    data=Abundance_Osmia_2)

# Model diagnostics:
overdisp.glmer(Osmia_Abund_GLMM_Y2)
testDispersion(Osmia_Abund_GLMM_Y2)
simulationOutput <- simulateResiduals(fittedModel = Osmia_Abund_GLMM_Y2, plot = F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
#No problems detected.

#Results (Year 2 only):
summary(Osmia_Abund_GLMM_Y2)

# Both Years, apple bloom period only, With Year as fixed categorical variable, including interaction:

Osmia_Abund_GLMM_Year_All <- glmer(Count ~ Plot * as.factor(Year) + (1|`Orchard ID`/Site), 
                    family=poisson(link = "log"),
                    data=Abundance_Osmia_AB, 
                    control=glmerControl(optimizer="nloptwrap"))

# Model diagnostics:
overdisp.glmer(Osmia_Abund_GLMM_Year_All)
testDispersion(Osmia_Abund_GLMM_Year_All)
simulationOutput <- simulateResiduals(fittedModel = Osmia_Abund_GLMM_Year_All, plot = F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
#No problems detected.

#Results (both years):
summary(Osmia_Abund_GLMM_Year_All)

```

#### 1.2 All Flower Visitors

```{r Summarise Pollinator Abundance Mean, SE, and % Difference During Apple Bloom, Both Years}

#Means and SEs, across the two years (48 observations):
Abundance_Total_AB %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Count), list(Mean_Count = mean, SE_Count = se))

# Percent Diff
Abundance_Total_AB_Percent_Diff <- Abundance_Total_AB %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Count), list(Mean_Count = mean, SE_Count = se))

Abundance_Total_AB_Percent_Diff %>%
  summarise(percent_diff = (100*(Mean_Count[2]-Mean_Count[1])/mean(Mean_Count)))

# Calculate mean % difference in pollinator numbers during apple bloom, both years combined.
#Orchard ID column has to be removed.
Abundance_Total_Values_Diff <- pivot_wider(Abundance_Total_AB, id_cols=c(Site, Year), names_from = Plot, values_from = Count) 
Abundance_Total_Values_Diff$Percent_Diff <- (Abundance_Total_Values_Diff$Treatment-Abundance_Total_Values_Diff$Control)*100 / ((Abundance_Total_Values_Diff$Control+Abundance_Total_Values_Diff$Treatment)/2)

#Mean percent difference in pollinator abundance:
mean(Abundance_Total_Values_Diff$Percent_Diff)
#Std. error of percent difference in bee abundance:
sd(Abundance_Total_Values_Diff$Percent_Diff, na.rm = T)/sqrt(length(Abundance_Total_Values_Diff$Percent_Diff))

```

##### 1.2.1 Models

```{r Total Pollinator Abundance with Year Model}

# With Year as Fixed categorical variable
hist(sqrt(Abundance_Total_AB$Count))

Total_Abund_GLMM_Year_All <- lmer(sqrt(Count) ~ Plot * as.factor(Year) + (1|`Orchard ID`/Site), 
                    data=Abundance_Total_AB)

# Model diagnostics:
overdisp.glmer(Total_Abund_GLMM_Year_All)
testDispersion(Total_Abund_GLMM_Year_All)
simulationOutput <- simulateResiduals(fittedModel = Total_Abund_GLMM_Year_All, plot = F)
plot(simulationOutput)

#Results (both years):
summary(Total_Abund_GLMM_Year_All)

```

#### 1.3 Small Flower Visitors

```{r Summarise Small Flower Abundance Mean, SE, and % Difference During Apple Bloom, Both Years}

#Means and SEs, across the two years (48 observations):
Abundance_Select_AB %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Count), list(Mean_Count = mean, SE_Count = se))

# Percent Diff
Abundance_Select_AB_Percent_Diff <- Abundance_Select_AB %>%
  group_by(Plot) %>%
  dplyr::summarize_at(vars(Count), list(Mean_Count = mean, SE_Count = se))

Abundance_Select_AB_Percent_Diff %>%
  summarise(percent_diff = (100*(Mean_Count[2]-Mean_Count[1])/mean(Mean_Count)))

# Calculate mean % difference in pollinator numbers during apple bloom, both years combined.
#Orchard ID column has to be removed.
Abundance_Select_Values_Diff <- pivot_wider(Abundance_Select_AB, id_cols=c(Site, Year), names_from = Plot, values_from = Count) 
Abundance_Select_Values_Diff$Percent_Diff <- (Abundance_Select_Values_Diff$Treatment-Abundance_Select_Values_Diff$Control)*100 / ((Abundance_Select_Values_Diff$Control+Abundance_Select_Values_Diff$Treatment)/2)

#Mean percent difference in pollinator abundance:
mean(Abundance_Select_Values_Diff$Percent_Diff)
#Std. error of percent difference in bee abundance:
sd(Abundance_Select_Values_Diff$Percent_Diff, na.rm = T)/sqrt(length(Abundance_Select_Values_Diff$Percent_Diff))

```

##### 1.3.1 Models

```{r Small Flower Visitor Abundance with Year Model}

# Both Years, apple bloom period only, With Year as fixed categorical variable, including interaction:

Select_Abund_GLMM_Year_All <- glmer(Count ~ Plot * as.factor(Year) + (1|`Orchard ID`/Site), 
                    family=poisson(link = "log"),
                    data=Abundance_Select_AB, 
                    control=glmerControl(optimizer="nloptwrap"))

# Model diagnostics:
overdisp.glmer(Select_Abund_GLMM_Year_All)
testDispersion(Select_Abund_GLMM_Year_All)
simulationOutput <- simulateResiduals(fittedModel = Select_Abund_GLMM_Year_All, plot = F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
#Zero-inflated but not overdispersed.

#Zero-inflated Poisson using glmmTMB:
Select_Abund_ZI_All <- glmmTMB(Count ~ Plot * as.factor(Year) + (1|`Orchard ID`/Site), 
                    ziformula = ~1 , family = "poisson",
                    data=Abundance_Select_AB)

# Model diagnostics:
testDispersion(Select_Abund_ZI_All)
simulationOutput <- simulateResiduals(fittedModel = Select_Abund_ZI_All, plot = F)
plot(simulationOutput)
testZeroInflation(simulationOutput)
#No problems.

#Selected model results:
summary(Select_Abund_ZI_All)

```

#### 1.4 All Years

##### 1.4.1 Figures

```{r All Year Figure}

pdf(file="Figure 1.pdf", width=12, height=8)

#Panel A (Osmia):
Abund_Fig_Osmia <- ggplot(data=Abundance_Osmia_All, aes(x=Plot, y=Count, colour=Plot, fill=Plot)) +
  scale_colour_brewer(palette = "Set1") +
  scale_fill_manual(values = c(adjustcolor("#E41A1C", alpha = 0.5), adjustcolor("#377EB8", alpha = 0.5)), aesthetics = "fill") +
  geom_boxplot(width = 0.5) +
  geom_jitter(height=0, width = 0.1, size = 1) +
  theme_bw() +
  labs(y = expression(paste(italic("Osmia"), " observed")), x = "") +
  scale_x_discrete(labels=c("control", "trapnests")) +
  scale_y_continuous(trans = 'sqrt') +
  facet_grid( ~ fct_relevel(Round_Year, "Apple bloom 2021", "Pre-thinning 2021", "Post-thinning 2021")) +
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text=element_text(size=14),
    axis.title=element_text(size=16),
    strip.text.x = element_text(size = 16),
    legend.position="none",
    strip.background = element_rect(color="NA", fill="NA")
    )
Abund_Fig_Osmia 

#Panel B (all pollinators):
Abund_Fig_Bees <- ggplot(data=Abundance_Total_All, aes(x=Plot, y=Count, colour=Plot, fill=Plot)) +
  scale_colour_brewer(palette = "Set1")+
  scale_fill_manual(values = c(adjustcolor("#E41A1C", alpha = 0.5), adjustcolor("#377EB8", alpha = 0.5)), aesthetics = "fill") +
  geom_boxplot(width = 0.5) +
  geom_jitter(height=0, width = 0.1, size = 1) +
  theme_bw() +
  labs(y = "pollinators observed", x = "") +
  scale_x_discrete(labels=c("control", "trapnests")) +
  scale_y_continuous(trans = 'sqrt') +
  facet_grid( ~ fct_relevel(Round_Year, "Apple bloom 2021", "Pre-thinning 2021", "Post-thinning 2021")) +
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text=element_text(size=14),
    axis.title=element_text(size=16),
    legend.position="none",
    strip.background = element_rect(color="NA",fill="NA"),
    strip.text.x = element_blank()
    )
Abund_Fig_Bees 

#Merge panels:
plot_grid(Abund_Fig_Osmia, Abund_Fig_Bees, nrow = 2, labels = c("a","b"), rel_heights = c(1.1,1))

```

### 2. Apple Yield and Quality

#### 2.1 Pre-Thinning

```{r Apple Pre_Thinning Mean and se}

# Mean and SE of fruit-set pre-thinning values per plot
Apple_Values_Pre_2 <- Apple_2 %>%
  group_by(Plot) %>%
  dplyr::summarize(Mean_Fruit_Set_Pre=mean(Fruit_Set_Pre, na.rm = T),
                   SE_Fruit_Set_Pre=se(Fruit_Set_Pre))
Apple_Values_Pre_2


```

##### 2.1.1 Model

```{r Apple Yield Pre-Thinning Model Year 2}

# Binomial GLMER With Site as a Random Effect

Apple_Pre_Thin_GLM <- glmer(Pre_Fruit/Buds ~ Plot + (1|`Orchard ID`/Site), 
                    family = binomial(link = "logit"),
                    data=Apple_2)
summary(Apple_Pre_Thin_GLM)

# Model diagnostics:
overdisp.glmer(Apple_Pre_Thin_GLM)
testDispersion(Apple_Pre_Thin_GLM)
simulationOutput <- simulateResiduals(fittedModel = Apple_Pre_Thin_GLM, plot = T)
#Fit is terrible.

#Ordered beta regression with Site random:
Apple_Pre_Thin_GLM <- glmmTMB(Fruit_Set_Pre ~ Plot + (1|`Orchard ID`/Site), 
                              dispformula = ~1,                  
                              family = ordbeta,
                    data=Apple_2)
summary(Apple_Pre_Thin_GLM)

# Model diagnostics:
testDispersion(Apple_Pre_Thin_GLM)
simulationOutput <- simulateResiduals(fittedModel = Apple_Pre_Thin_GLM, plot = T)
#Much better.

```

#### 2.2 All Years

##### 2.2.1 Figures

```{r Apple All Years Figure}

pdf(file="Figure 2.pdf", width=8, height=8)

#Panel A (fruit set pre-thinning, 2021):
Fruit_Set_Pre_Fig_2021 <- ggplot(data=subset(Apple_All, Year == 2021), aes(x=Plot, y=Fruit_Set_Pre, colour=Plot, fill=Plot)) +
  scale_colour_brewer(palette = "Set1")+
  scale_fill_manual(values = c(adjustcolor("#E41A1C", alpha = 0.5), adjustcolor("#377EB8", alpha = 0.5)), aesthetics = "fill") +
  geom_boxplot(width = 0.5) +
  geom_jitter(height=0, width = 0.1, size = 1) +
  theme_bw() +
  labs(y = "fruit set (pre-thinning)", x = "") +
  ylim(0,1) +
  scale_x_discrete(labels=c("control", "trapnests")) +
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text=element_text(size=14),
    axis.title=element_text(size=16),
    legend.position="none",
    strip.background = element_rect(color="NA",fill="NA"))
Fruit_Set_Pre_Fig_2021

# Panel B (fruit set pre-thinning, 2022):
Fruit_Set_Pre_Fig_2022 <- ggplot(data=subset(Apple_All, Year == 2022), aes(x=Plot, y=Fruit_Set_Pre, colour=Plot, fill=Plot)) +
  scale_colour_brewer(palette = "Set1")+
  scale_fill_manual(values = c(adjustcolor("#E41A1C", alpha = 0.5), adjustcolor("#377EB8", alpha = 0.5)), aesthetics = "fill") +
  geom_boxplot(width = 0.5) +
  geom_jitter(height=0, width = 0.1, size = 1) +
  theme_bw() +
  labs(y = "fruit set (pre-thinning)", x = "") +
  ylim(0,1) +
  scale_x_discrete(labels=c("control", "trapnests")) +
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text=element_text(size=14),
    axis.title=element_text(size=16),
    legend.position="none",
    strip.background = element_rect(color="NA",fill="NA"))
Fruit_Set_Pre_Fig_2022

#Panel C (fruit set post-thinning, 2021):
Fruit_Set_Post_Fig_2021 <- ggplot(data=subset(Apple_All, Year == 2021), aes(x=Plot, y=Fruit_Set_Post, colour=Plot, fill=Plot)) +
  scale_colour_brewer(palette = "Set1")+
  scale_fill_manual(values = c(adjustcolor("#E41A1C", alpha = 0.5), adjustcolor("#377EB8", alpha = 0.5)), aesthetics = "fill") +
  geom_boxplot(width = 0.5) +
  geom_jitter(height=0, width = 0.1, size = 1) +
  theme_bw() +
  labs(y = "fruit set (post-thinning)", x = "") +
  ylim(0,1) +
  scale_x_discrete(labels=c("control", "trapnests")) +
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text=element_text(size=14),
    axis.title=element_text(size=16),
    legend.position="none",
    strip.background = element_rect(color="NA",fill="NA"))
Fruit_Set_Post_Fig_2021

#Panel D (apple circumference, 2021):
Fruit_Qual_Fig_2021 <- ggplot(data=subset(Apple_All, Year == 2021), aes(x=Plot, y=Mean_Circum, colour=Plot, fill=Plot)) +
  scale_colour_brewer(palette = "Set1")+
  scale_fill_manual(values = c(adjustcolor("#E41A1C", alpha = 0.5), adjustcolor("#377EB8", alpha = 0.5)), aesthetics = "fill") +
  geom_boxplot(width = 0.5) +
  geom_jitter(height=0, width = 0.1, size = 1) +
  theme_bw() +
  labs(y = "fruit circumference (cm)", x = "") +
  scale_x_discrete(labels=c("control", "trapnests")) +
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text=element_text(size=14),
    axis.title=element_text(size=16),
    legend.position="none",
    strip.background = element_rect(color="NA",fill="NA"))
Fruit_Qual_Fig_2021

# Merge panels:
plot_grid(Fruit_Set_Pre_Fig_2021, Fruit_Set_Pre_Fig_2022, Fruit_Set_Post_Fig_2021, Fruit_Qual_Fig_2021, nrow = 2, labels = c("a","b", "c", "d"), rel_heights = c(1,1))

```

# Supporting Information

## Data Organization

### 1. Trap Nests

#### 1.1 Holes Occupied

```{r Data Trap Nests Summary}

# Load Excel File with Data from All Rounds
Trap_Nest_Total <- read_excel("Trap_Nest_Total.xlsx", na = "N/A") %>%
    unite('Site', Orchard, 'Site ID', remove = FALSE)

# Load Excel File with Verified Trap Nest Cell Counts and Numbers (Final Round Only)
Trap_Nest_Final <- read_excel("Trap_Nest_Final_Count.xlsx", na = "N/A") %>%
    unite('Site', Orchard, 'Site ID', remove = FALSE)

# Total Number of Holes in Trap Nests
Total_Holes <- Trap_Nest_Final %>% filter(Round == "Final Count") %>%
  nrow()
Total_Holes

# Total Number of Occupied Holes in Trap Nests
Total_Occupied_Holes <- Trap_Nest_Final %>% filter(Round == "Final Count") %>%
  drop_na(`# New Cells`) %>%
  nrow()
Total_Occupied_Holes

# Percent of holes occupied (overall)
Total_Occupied_Holes/Total_Holes

# Total Number of Occupied Holes in Trap Nests per Taxon per Site
Holes_Taxon_Site <- Trap_Nest_Final %>% filter(Round == "Final Count") %>%
  mutate(Taxa=replace_na(Taxa, "Unoccupied")) %>%
  group_by(Site, Taxa, .drop=FALSE) %>%
  do(data.frame(nrow=nrow(.))) %>% 
  ungroup %>%
  complete(Site, nesting(Taxa)) %>%
  arrange(Taxa) %>%
  mutate(nrow = replace_na(nrow, 0))
Holes_Taxon_Site

# Total Number of Occupied Holes in Trap Nests per Taxon 
Holes_Taxon <- Trap_Nest_Final %>% filter(Round == "Final Count") %>%
  mutate(Taxa=replace_na(Taxa, "Unoccupied")) %>%
  count("Taxa") %>%
  mutate("%" = freq/sum(freq))
Holes_Taxon

```


#### 1.2 Brood Cells Produced

Calculating trap nest offspring numbers

```{r Data Trap Nests Offspring, results = "hide"}

# Calculating total number of brood cells across all Taxa
Trap_Nest_Final %>%
  dplyr::summarize_at(vars(`# Total Cells`), list(Total_Cells_Sum=sum)) %>%
  drop_na(Total_Cells_Sum) 

# Calculating percentage of total cells produced per Taxon
Cells_Taxon_Final <- Trap_Nest_Final %>%
  group_by(Taxa) %>% 
  dplyr::summarize_at(vars(`# Total Cells`), list(Total_Cells_Sum=sum)) %>%
  mutate("%" = Total_Cells_Sum/sum(Total_Cells_Sum)) %>%
  drop_na(Total_Cells_Sum) 
Cells_Taxon_Final

# Calculating percentage of new cells produced during apple bloom, Round 2 (pre-apple thinning), and Round 3 (post-apple thinning)
Cells_Taxon_Round_Total <- Trap_Nest_Total %>%
  group_by(Taxa, Round) %>% 
  dplyr::summarize_at(vars(`# New Cells`), list(New_Cells_Sum=sum)) %>%
  filter(Round != "Round 4" & Round != "Final Count") %>%
  mutate("%" = New_Cells_Sum/sum(New_Cells_Sum)) %>%
  drop_na(New_Cells_Sum) 
Cells_Taxon_Round_Total

# Calculating percentage of new cells produced during apple bloom, Round 2 (pre-apple thinning), and Round 3 (post-apple thinning) per Site
Cells_Taxon_Round_Site_Total <- Trap_Nest_Total %>%
  mutate(Site = as.factor(Site), Taxa = as.factor(Taxa), Round = as.factor(Round)) %>%
  group_by(Site, Taxa, Round, .drop = FALSE) %>% 
  dplyr::summarize(Total_New_Cells = sum(`# New Cells`, na.rm = TRUE)) %>%
  filter(Round != "Round 4" & Round != "Final Count") %>%
  mutate("%" = Total_New_Cells/sum(Total_New_Cells)) %>%
  drop_na()
Cells_Taxon_Round_Site_Total

# Re-labeling Round 2 and Round 3
Cells_Taxon_Round_Site_Total$Round <- factor(Cells_Taxon_Round_Site_Total$Round, levels=c("Apple Bloom", "Round 2", "Round 3"), labels=c("Apple Bloom", "Pre-Apple Thinning", "Post-Apple Thinning"))

# Calculating percentage of total cells produced per Site
Cells_Site_Total <- Trap_Nest_Final %>%
  group_by(Taxa, Site) %>% 
  dplyr::summarize_at(vars(`# Total Cells`), list(Total_Cells=sum)) %>%
  mutate("%" = Total_Cells/sum(Total_Cells)) %>%
  drop_na()
Cells_Site_Total
```

#### 1.3 Edge Effects

```{r Data Edge Effects, results = "hide"}

# Loading Excel file
Edge_Effects <- read_excel("Trap_Nest_Edge_Effects.xlsx", na = "N/A")
Edge_Effects$Trap_Nest_ID <- Edge_Effects$`Trap Nest ID`

# Average number of holes per trap nest
Available_Holes_Mean <- Edge_Effects %>%
  dplyr::summarize(Mean=mean(`Trap Nest Hole #`, na.rm = T),
                   SD=sd(`Trap Nest Hole #`))
Available_Holes_Mean

# Edit database
Trap_Nest_Hole_Count <- as.data.frame(Trap_Nest_Final)
Trap_Nest_Hole_Count$`# New Cells`[is.na(Trap_Nest_Hole_Count$`# New Cells`)] <- 0
Trap_Nest_Hole_Count$`# New Cells`[(Trap_Nest_Hole_Count$`# New Cells` > 0)] <- 1
Trap_Nest_Hole_Count <- Trap_Nest_Hole_Count %>% group_by(`Trap Nest ID`, Site) %>%
  summarise_at(vars(`# New Cells`), list(Total_Holes=sum))
Trap_Nest_Hole_Count$Trap_Nest_ID <- Trap_Nest_Hole_Count$`Trap Nest ID`

Edge_Effects <- merge(Trap_Nest_Hole_Count, Edge_Effects, by="Trap_Nest_ID")
```

#### 1.4 Trap-Nest Offspring Survival

```{r Data Survival X-Ray, results = "hide"}

#Filter for Only Osmia
Osmia_Nest <- Trap_Nest_Final %>% filter(Taxa == "Osmia") %>%
  drop_na(`# New Cells`)

#How many Osmia nests per site, how many were x-rayed at each site, and how many total brood cells, surviving brood cells, and surviving female brood cells were there per x-rayed nest?
#Total brood cells calculated from full dataset includes destroyed nests (total cells recorded as 0);
#x-rayed cells excludes destroyed nests.
Osmia_Nest_Summary <- Osmia_Nest %>% group_by(Site) %>%
  dplyr::summarize(Mean_Total_Cells=mean(`# Total Cells`, na.rm = T), 
                   Sum_Total_Cells = sum(`# Total Cells`, na.rm = T),
                   N_Nests = n())

XRayed_Osmia_Nests <- filter(Osmia_Nest,`X-Ray Mount`!="NA")
XRay_Osmia_Survival <- XRayed_Osmia_Nests %>% group_by(Site) %>%
  dplyr::summarize(Total_XRayed_Cells=sum(`# Total Cells`, na.rm = T), 
                   Total_Surv_Female_Cells=sum(`# Surviving Female Osmia`, na.rm = T),
                   Total_Surv_Osmia_Cells=sum(`# Surviving Osmia`, na.rm = T),
                   N_XRayed_Nests = n())

Osmia_Nest_Survival <- join(Osmia_Nest_Summary, XRay_Osmia_Survival)
Osmia_Nest_Survival$Lost <- Osmia_Nest_Survival$N_Nests - Osmia_Nest_Survival$N_XRayed_Nests
Osmia_Nest_Survival$Cells_Per_XRayed_Nest <-
  Osmia_Nest_Survival$Total_XRayed_Cells/Osmia_Nest_Survival$N_XRayed_Nests
Osmia_Nest_Survival$Surv_Cells_Per_XRayed_Nest <-
  Osmia_Nest_Survival$Total_Surv_Osmia_Cells/Osmia_Nest_Survival$N_XRayed_Nests
Osmia_Nest_Survival$Surv_Females_Per_XRayed_Nest <-
  Osmia_Nest_Survival$Total_Surv_Female_Cells/Osmia_Nest_Survival$N_XRayed_Nests

```

#### 1.5 Pollen Usage

```{r Data Pollen, results = "hide"}

# Load Pollen Data into a dataframe
Pollen <- read_excel("Pollen.xlsx") %>%
    as.data.frame(Pollen)
```

## Results

### 2.1 Phenology of Trap Nest Offspring Production

#### 2.1.1 Figure

```{r Trap Nest Offspring Figure}

pdf(file="Figure S6.pdf", width=10, height=6)

ggplot(Cells_Taxon_Round_Site_Total, aes(x=as.factor(Round), y = as.numeric(Total_New_Cells), fill=Taxa)) +
  geom_boxplot() +
  geom_jitter(height=0, aes(y = as.numeric(Total_New_Cells)), size = 0.5) +
  scale_fill_manual(breaks = c("Chelostomoides", "Megachile", "Osmia", "Wasp"), values=c("white", "grey75", "grey45", "grey3"), labels=c("resin bees and wasps", "Megachile", "Osmia", "potter wasps")) +
  scale_color_manual(breaks = c("Chelostomoides", "Megachile", "Osmia", "Wasp"), values=c("white", "grey75", "grey45", "grey3"), labels=c("resin bees and wasps", "Megachile", "Osmia", "potter wasps")) +
  theme_bw() +
  labs(y="Number of New Cells Produced per Site", x="") +
  theme(
    axis.ticks.x=element_blank(),
    axis.text=element_text(size=14),
    axis.title=element_text(size=16, face="bold"),
    legend.text = element_text(size=14),
    legend.title = element_text(size=16))
```

### 2.2 Occupancy and Brood Cell Production by Taxon

#### 2.2.1 Holes Occupied per Taxon

##### 2.2.1.1 Figure

```{r Trap Nests per Taxon Figure}

pdf(file="Figure S4.pdf", width=8, height=6)

ggplot(Holes_Taxon_Site, aes(x=reorder(Taxa, nrow), y = nrow)) +
  geom_boxplot(fill="grey") +
  geom_jitter(height=0, size=0.25) +
  scale_x_discrete(labels = c("resin bees/wasps", "Megachile", "Osmia", "potter wasps", "unoccupied")) +
  theme_bw() +
  labs(y="Total Nest Holes per Site", x="Taxa") +
  theme(
    axis.text=element_text(size=14),
    axis.title=element_text(size=16, face="bold"))
```

#### 2.2.2 Brood Cells Produced per Taxon

##### 2.2.2.1 Figure

```{r Trap Nest Cells per Taxon Figure}

pdf(file="Figure S5.pdf", width=8, height=6)

ggplot(Cells_Site_Total, aes(x=reorder(Taxa, Total_Cells), y = Total_Cells))+
  geom_boxplot(fill="grey") +
  geom_jitter(height=0, size=0.25) +
  theme_bw() +
  labs(y="Total Cells per Site", x="Taxa") +
  theme(
    axis.text=element_text(size=14),
    axis.title=element_text(size=16, face="bold"))
```

### 2.3 Edge Effects

#### 2.3.1 Model

```{r Edge Effects on Nest Hole Occupancy (Numbers) Model }

# Re-labeling column
Edge_Type <- Edge_Effects$`Trap Nest Edge`

# Re-scaling data
Distance_Edge <- rescale(Edge_Effects$`Trap Nest Edge Distance (m)`)

# Running GLMM
Edge_Effects_Holes_GLMM <- glmer(Total_Holes ~ Edge_Type*Distance_Edge + (1|`Orchard ID`/Site), 
                                  family=poisson(link="log"), 
                                  data=Edge_Effects)
Edge_Effects_Holes_GLMM

summary(Edge_Effects_Holes_GLMM)

```

#### 2.3.2 Figures

```{r Edge Effects Figures}

pdf(file="Figure S9.pdf", width=8, height=6)

# Re-labeling column
Distance_Edge <- Edge_Effects$`Trap Nest Edge Distance (m)`

# Creating Figure
Figure_S9 <- ggplot(data = Edge_Effects, aes(x = Distance_Edge, y = Total_Holes)) +
  geom_point(aes(shape=Edge_Type, col=Edge_Type)) + 
  geom_smooth(aes(col=Edge_Type), method=lm) +
  scale_color_brewer(palette="Set1", name ="Edge Type") +
  scale_shape_manual(values = c(0,16,3), name = "") +
  theme_bw() +
  labs(y="Total Number of Holes", x="Distance to Edge (m)") +
  theme(
    axis.ticks.x=element_blank(),
    axis.text=element_text(size=13),
    axis.title=element_text(size=15, face="bold"))
Figure_S9
```

### 2.4 Trap-Nest Offspring Survival

For Year 1 of *Osmia* spp. offspring as determined by X-ray.

```{r Trap Nest Cell Survival}

#What percentage of all nests were completely destroyed?
Total_Lost <- sum(Osmia_Nest_Survival$Lost); Total_Nests <- sum(Osmia_Nest_Survival$N_Nests)
Total_Lost/Total_Nests

#What proportion of all surviving Osmia brood cells were female?
sum(Osmia_Nest_Survival$Total_Surv_Female_Cells)/sum(Osmia_Nest_Survival$Total_Surv_Osmia_Cells)

#How many total brood cells per nest on average?
mean(Osmia_Nest_Survival$Cells_Per_Nest);
sd(Osmia_Nest_Survival$Cells_Per_Nest)

#Mean and SD of number of surviving females per nest:
mean(Osmia_Nest_Survival$Surv_Females_Per_Nest); 
sd(Osmia_Nest_Survival$Surv_Females_Per_Nest)

```

#### 2.4.1 Figures

```{r X-Ray Osmia Survival Figure}

Pivot_Osmia_Nest_Survival <- pivot_longer(Osmia_Nest_Survival, cols = 10:12, names_to = "cell_type", values_to = "numbers")

pdf(file="Figure S7.pdf", width=8, height=6)

Figure_S7 <- ggplot(Pivot_Osmia_Nest_Survival, aes(x=factor(cell_type, level = c("Cells_Per_XRayed_Nest", "Surv_Cells_Per_XRayed_Nest", "Surv_Females_Per_XRayed_Nest")), y = numbers))+
  geom_boxplot(width = 0.5) +
  scale_fill_brewer(palette="YlOrBr") +
  geom_jitter(aes(cex = N_XRayed_Nests), height=0, width=0.1, pch = 21) +
  theme_bw() +
  labs(y="Mean Number of Offspring per Nest", x="") +
  scale_x_discrete(labels=c("Total", "Total Surviving", "Surviving Females")) +
  theme(
    legend.position = "None",
    axis.text=element_text(size=14),
    axis.title=element_text(size=16, face="bold")) 
Figure_S7

```

### 2.5 Pollen Usage

```{r Pollen Percentage}

# Summarizing pollen data by percentage
Pollen_Values <- Pollen %>%
  group_by(Taxa) %>% 
  dplyr::summarize_at(vars(Number_Grains), list(Number_Grains=sum)) %>%
  mutate("%" = Number_Grains/sum(Number_Grains)) 

```