##This is code for the analyses in the manuscript "Prescribed fire increases the number of ground-nesting bee nests in tallgrass prairie remnants"

#Loading packages used----
pacman::p_load(lme4, DHARMa, lmerTest, tidyverse, 
               vegan, ggplot2, GGally)

#reading in data frame of master file ----
master = read.csv("masterfile.csv")
#changing sample round from integer to factor
master$round<-as.factor(master$round)

#reading in data frames for PERMANOVA analyses of community similarity ----
#Bee community
#master file of sites and treatments
bn2nmdsm = read.csv("bn2_nmds_m.csv")
#bee community matrix
beematrix1<-read.csv("beematrix1.csv", row.names=1)
#floral community
flrmatrix1<-read.csv("flowermatrix1.csv", row.names=1)

#reading in dataframe for analysis of slope----
slope = read.csv("slopetable.csv", header=T)

#Correlation analysis ----
#plot of correlation
master %>% ggpairs(columns = c("mean_bare", "mean_veg", "mean_thatch", "mean_litter", "mean_rocks", "mean_thatchdepth"),
                             upper = list(continuous = wrap('cor', size=5)))

#Analysis of sum number of active nests----
MOD_negbin<-glmer.nb(sum.nests ~ treatment + (1|site), data=master)
simulationOutput <- simulateResiduals(fittedModel = MOD_negbin, plot = T)
summary(MOD_negbin) 

#checking for treatment effects with the model to see if the bare ground found in the burned areas was driving the result
resid = resid(MOD_negbin)
plot(master$mean_bare, resid) #appears evenly dispersed, no correlation between increasing mean bare ground and model residuals


#Analysis of the ENS of bees ----

m1_bens = lmer(bee.ens~treatment + (1|site), data=master)
simulationOutput <- simulateResiduals(fittedModel = m1_bens, plot = T) 
summary(m1_bens)

#Analysis of the average number of flowers ----

m2_nflrs = lmer(log(avg.no.flrs+1)~treatment +(1|site), data=master)
simulationOutputm1 <- simulateResiduals(fittedModel = m2_nflrs, plot = T)
summary(m2_nflrs) 

#Analysis of the number of flower species ----

m3_nflrsp = glmer(sum.flr.sp~treatment +(1|site), data=master, family = poisson(link = "log")) 
simulationOutputm1 <- simulateResiduals(fittedModel = m3_nflrsp, plot = T)
summary(m3_nflrsp)

#Analysis of the ENS of of flower species ----

m4_flrens=lmer(flr.ens~treatment+(1|site), data=master) 
simulationOutput <- simulateResiduals(fittedModel = m4_flrens, plot = T)
summary(m4_flrens)

#Analysis of the average percent bare ground ----

m5_bg = lmer(asin(sqrt(mean_bare/100)) ~ treatment + (1|site), data=master)
simulationOutputm3 <- simulateResiduals(fittedModel = m5_bg, plot = T)
summary(m5_bg)

#Analysis of the average percent vegetative cover ----

m6_vc = lmer((mean_veg/100)~treatment + (1|site), data=master)
simulationOutputm1 <- simulateResiduals(fittedModel = m6_vc, plot = T)
summary(m6_vc)

#Analysis of average thatch depth ----

m7_td = lmer(log(mean_thatchdepth+1)~treatment + (1|site), data=master)
simulationOutputm2 <- simulateResiduals(fittedModel = m7_td, plot = T)
summary(m7_td)

#Analysis of slope ----

#Just east, southeast, southwest, and south facing slopes

chisq.test(rbind(c(52, 90), c(86,103)))

### Community similarity analysis of bees ----
#permanova analysis
bn2.dist <-vegdist(beematrix1, method="bray")
perm <- how(nperm = 1000)
setBlocks(perm) <- with(bn2nmdsm, pair)
adonis2(bn2.dist~trt, data=bn2nmdsm, permutations = perm)
#note that this gives an error because there aren't enough pairs/sites to do as many permutations as the default
#error just means the analysis did as many permutations as was possible (15), likely is OK for this analysis bc so few sites

### Community similarity analysis of flowers ----
flrdist = vegdist(flrmatrix1, method="bray")
adonis2(flrdist~trt, data=bn2nmdsm, permutations = perm)