##### LOCAL DIVERSITY ANALYSIS #################


spec<- read.csv("Lane et al - SpeciesData.csv")

site.round<- read.csv("Lane et al - SitebyRound.csv")


# ------------ Calculating Effective Number of Species from Probability of Interspecific Encounter

library(tidyverse) # packages for shaping data frames
library(reshape2)

spec.tally<- spec %>% # tallying number of bee species by each unique sampling round / site combination
  group_by(rd_site,bee) %>%
  tally()

ens_matrix<- acast(spec.tally, rd_site~bee, value.var="n") # creating a round/site by species matrix
ens_matrix[is.na(ens_matrix)] <- 0 # blank entries set to 0

library(devtools) # packages for installing mobr from github
library(dplyr)
install_github('MoBiodiv/mobr')
library(mobr)# package for calculating PIE and ENS

ens.1<- as.data.frame(calc_PIE(ens_matrix,ENS=TRUE)) # calculating ENS of PIE

ens.2<- tibble::rownames_to_column(ens.1, "rd_site") # turning site/round into a distinct column

colnames(ens.2)[colnames(ens.2)=="calc_PIE(ens_matrix, ENS = TRUE)"] <- "ens" # changing column name to something more tractable

ens.2$ens.r<- round(ens.2$ens,digits=0) # rounding the ens values to closest whole number

ens.3<- merge(ens.2,site.round,by="rd_site") # merging calculation with site data

ens.3$round<- as.factor(ens.3$round) # making round a factor for random effects


# Generating Model


library(lme4) # package for glmer's

model.ens<- glmer(ens.r~mat+ha+forb_rich + primary_ag_1500 + (1|round/site), data=ens.3,family="poisson") # model

summary(model.ens) # output



# Model validation

library(ade4) # package that performs Mantel tests

site.dists <- dist(cbind(ens.3$longitude, ens.3$latitude)) # Checking spatial auto-correlation
ens.dists <- dist(ens.3$ens.r)


mt<-mantel.rtest(site.dists, ens.dists, nrepet = 9999) # calculating mantel test. Will generate some warnings because it tries to compare the same site across rounds


library(DHARMa) # package for checking model fit

simulationOutput<- simulateResiduals(fittedModel = model.ens, n = 250) # simulating residuals to test model fit

testResiduals(simulationOutput) # model fit tests - model is good fit for data


# Creating prediction data frames for plotting and plots

library(ggeffects) # package for creating prediction dataframes
library(emmeans) # package for creating prediction dataframes
library(ggplot2) # for creating plots
library(gridExtra) # for packaging plots

ag.predict<- ggemmeans(model.ens,terms="primary_ag_1500[[17:87]]",full.data = FALSE) # prediction data for surrounding ag plots

colnames(ag.predict)[colnames(ag.predict)=="x"] <- "primary_ag_1500" # renaming column names
colnames(ag.predict)[colnames(ag.predict)=="predicted"] <- "ens"


rich.predict<- ggemmeans(model.ens,terms="forb_rich",full.data = FALSE) # prediction data for floral resource plots

colnames(rich.predict)[colnames(rich.predict)=="x"] <- "forb_rich" # renaming column names
colnames(rich.predict)[colnames(rich.predict)=="predicted"] <- "ens"

# Agriculture

Agrich.plot<- ggplot(ens.3,aes(x=primary_ag_1500,y=ens)) +
  theme_classic() +
  geom_point(size=3) +
  theme(axis.text.y=element_text(size=14), axis.title.y=element_text(size=18)) +
  ylab("Effective Number of Species") +
  xlab("Proportion Agriculture at 1500m") +
  theme(axis.title.x = element_text(size=18),axis.text.x=element_text(size=14)) +
  geom_line(data=ag.predict, aes(x=primary_ag_1500,y=ens),linetype="dashed",color="black",size=1) +
  labs(title = "A.") +
  theme(plot.margin=unit(c(1.5,1,1,1), "cm")) +
  theme(axis.title.y = element_text(margin=margin(0,15,0,0)), axis.title.x = element_text(margin=margin(15,0,0,0)) ) +
  theme(plot.title = element_text(hjust = -0.15, vjust=9, face="bold", size=18))

# Forb Rich

ForbRich.plot<- ggplot(ens.3,aes(x=forb_rich,y=ens)) +
  theme_classic() +
  geom_point(size=3) +
  theme(axis.text.y=element_text(size=14), axis.title.y=element_text(size=18)) +
  ylab("Effective Number of Species") +
  xlab("Floral Resource Richness") +
  theme(axis.title.x = element_text(size=18),axis.text.x=element_text(size=14)) +
  geom_line(data=rich.predict, aes(x=forb_rich,y=ens)) +
  geom_ribbon(data=rich.predict, aes(ymin=conf.low,ymax=conf.high),alpha=0.3) +
  labs(title = "B.") +
  theme(plot.margin=unit(c(1.5,1,1,1), "cm")) +
  theme(axis.title.y = element_text(margin=margin(0,15,0,0)), axis.title.x = element_text(margin=margin(15,0,0,0)) ) +
  theme(plot.title = element_text(hjust = -0.15, vjust=9, face="bold", size=18))



grid.arrange(Agrich.plot, ForbRich.plot, ncol=2) # final plots



########### BETA DIVERSITY ANALYSIS ###########

library(vegan) # for calculating dissimilarities
library(betapart) # for calculating and comparing beta diversity

site<- read.csv("Lane et al - Site.csv")

site<- site[order(site$site),] # ordering the site file alphabetically by site name for future comparison with a matrix

#------------------- setting up data

beta.tally<- spec %>% # creating a tally data frame of bee species by site
  group_by(site,bee) %>%
  tally()


beta_matrix<- acast(beta.tally, site~bee, value.var="n") #creating a matrix from tallied results

beta_matrix[is.na(beta_matrix)] <- 0 # completing matrix by inserting 0 for the NAs

site.1<- site %>% mutate(ag_group = case_when(
  primary_ag_1500 <= 33 ~ "low-ag", 
  primary_ag_1500 > 33 & primary_ag_1500 <= 66 ~ "med-ag",
  primary_ag_1500 > 66 ~ "high-ag")) # creating categories for high and low agricultural development around the restorations

site.2<- site.1 %>% mutate(forb_group = case_when(
  forb_rich <= 24  ~ "Low-rich",
  forb_rich >= 26 & forb_rich <= 29 ~ "med-rich",
  forb_rich >= 28 ~ "High-rich")) # creating categories for high and low agricultural development around the restorations

site.2$inter<- interaction(site.2$forb_group,site.2$ag_group) # creating an interaction column between proportion agriculture and floral resource richness

dist<-vegdist(beta_matrix, index.family = "bray") # creating the site pairwise dissimilarity matrix for bee communities using Bray-Curtis dissimilarity

# conducting analysis

beta.mod1 <- with(site.2, betadisper(dist, inter)) # creating beta diversity model for interaction term
anova(beta.mod1) # ANOVA of beta diversity interaction term - not significant

beta.mod2 <- with(site.2, betadisper(dist, ag_group)) # creating beta diversity model for surrounding agriculture categories term
anova(beta.mod2) # ANOVA of beta diversity amongs sites in high and low categories of surrounding agricultural cover, term - not significant

beta.mod3 <- with(site.2, betadisper(dist, forb_group)) # creating beta diversity model for floral resource richness categories term
anova(beta.mod3) # NOVA of beta diversity amongs sites in high and low categories of floral resource richness, term - not significant


# creating graphs

# generating a dataframe that ggplot2 can handle

beta.d<- as.data.frame(beta.mod2$distances)

beta.d1 <- tibble::rownames_to_column(beta.d, "site")

colnames(beta.d1)[colnames(beta.d1)=="beta.mod2$distances"] <- "distances"

beta.f<- merge(site.2,beta.d1,by="site")

beta.f %>% 
  group_by(ag_group) %>%
  summarize(distance=mean(distances)) # mean of the two categories

beta.f$ag_group <- factor(beta.f$ag_group,
                          levels = c('low-ag','med-ag','high-ag'),ordered = TRUE)

b1<- ggplot(beta.f, aes(x=ag_group,y=distances)) + # ag group plot
  theme_classic() +
  geom_boxplot(size=1.5) +
  theme(axis.text.y=element_text(size=14), axis.title.y=element_text(size=18)) +
  ylab("Bray-Curtis Dissimilarity") +
  xlab("Surrounding Agricultural Cover") +
  theme(axis.title.x = element_text(size=18),axis.text.x=element_text(size=16)) +
  theme(plot.margin=unit(c(1.5,1.1,1.1,1.1), "cm")) +
  theme(axis.title.y = element_text(margin=margin(0,15,0,0)), axis.title.x = element_text(margin=margin(15,0,0,0)) ) +
  theme(legend.position = "none") +
  labs(title = "A.") +
  theme(plot.title = element_text(hjust = -0.15, vjust=9, face="bold", size=18)) +
  scale_x_discrete(labels = c('Low','Medium','High'))

# generating a dataframe that ggplot2 can handle

beta.d.f<- as.data.frame(beta.mod3$distances)

beta.d2 <- tibble::rownames_to_column(beta.d.f, "site")

colnames(beta.d2)[colnames(beta.d2)=="beta.mod3$distances"] <- "distances"

beta.f2<- merge(site.2,beta.d2,by="site")

beta.f2 %>% 
  group_by(forb_group) %>%
  summarize(distance=mean(distances)) # mean of the two categories

beta.f2$forb_group <- factor(beta.f2$forb_group,
                             levels = c('Low-rich','med-rich','High-rich'),ordered = TRUE)



b2<-ggplot(beta.f2, aes(x=forb_group,y=distances)) + # floral resource richness plot
  theme_classic() +
  geom_boxplot(size=1.5) +
  theme(axis.text.y=element_text(size=14), axis.title.y=element_text(size=18)) +
  ylab("Bray-Curtis Dissimilarity") +
  xlab("Floral Resource Richness") +
  theme(axis.title.x = element_text(size=18),axis.text.x=element_text(size=16)) +
  theme(plot.margin=unit(c(1.5,1.1,1.1,1.1), "cm")) +
  theme(axis.title.y = element_text(margin=margin(0,15,0,0)), axis.title.x = element_text(margin=margin(15,0,0,0)) ) +
  theme(legend.position = "none")+
  labs(title = "B.") +
  theme(plot.title = element_text(hjust = -0.15, vjust=9, face="bold", size=18)) +
  scale_x_discrete(labels = c('Low','Medium','High'))


grid.arrange(b1, b2, ncol=2) # final figure