# This script conducts the analyses and creates the plots presented in Ritchie et al 2020 Restoration Ecology.
# meta data are available from the Digital Repository for the University of Minnesota (DRUM)
# in seed_land
# each row represents a fruit (fruit ID is unique identifier)
# ID identifies each unique plant
# prop.c = proportion surrounding agriculture
# 


library(geosphere)
library(lme4)
library(nlme)
library(glmmTMB)
library(MASS)
library(car)
library(tidyverse)
library(sjPlot)
library(ggeffects)
library(sjlabelled)
library(snakecase)
library(DHARMa)

#upload seed_land dataframe; seed counts for each fruit collected from each plant; a row represents a fruit.
seed_land<-read.csv("seed_land.csv")



##########################
# General data summaries

#make site data frame to look at distances of sites

sites<-seed_land%>%select(site, prop.c,longitude,latitude,ha)%>%group_by(site,prop.c,longitude,latitude,ha)%>%summarize()
###mean, max, min distances between sites using geosphere
#distm(c(site$latitude,sites$longitude),fun=distHaversine)

site_coord<-sites%>%select(longitude,latitude)

site_dist<-distm(site_coord,fun=distHaversine)
###
##data presented in ms under study design section
mean(site_dist)
max(site_dist)
min(site_dist)
max(sites$prop.c)
min(sites$prop.c)
mean(sites$prop.c)
sd(sites$prop.c)
max(sites$ha)
min(sites$ha)
mean(sites$ha)
#correlation test for size and surrounding ag
cor.test(sites$ha,sites$prop.c)


############################
#General summary info  for seed collection

# # of distinct plants we collected seed from
n_distinct(seed_land$ID)


# seed and fruit data summaryS


# Table of total fruit and seeds produced per plant 
# (same framework as plt_nr table in models code)
tot_seed_tbl<-seed_land%>%
  group_by(trmnt,ID)%>%
  summarise(sum_seeds=sum(total_seeds), n_fruit=n(), seeds_fruit=sum_seeds/n_fruit)%>%
  group_by(trmnt)%>%
  summarise(n_plants=n(),total_fruit=sum(n_fruit),sd_fruit=sd(n_fruit),mean_seeds=mean(seeds_fruit),
            mean_fruit=mean(n_fruit), sum_seeds=sum(sum_seeds),sd_mean_seeds=sd(seeds_fruit),
            avg_sum_seeds_per_plant=mean(sum_seeds), med_sum_seeds_per_plant=median(sum_seeds))
# mean_seeds, n_plants, total_fruit, and mean_fruit 
# are reported in first paragraph of results

## 95CIs for mean # fruit per treatment; presented in results section of MS
### These are reported  in first paragraph of results
# sp 
u_sp <- 5.17
sd_s <- 3.67
n_s<- 62
error_s <- qnorm(0.975)*sd_s/sqrt(n_s)
left <- u_sp-error_s
right <- u_sp+error_s
left
error_s
right
#OP
u_op <- 4.25
sd <- 2.81
n <- 51
error <- qnorm(0.975)*sd/sqrt(n)
left <- u_op-error
right <- u_op+error
left
error
right

################################################################ 
#95Confidence interval for mean seeds/fruit 


u_sp <- 8.68
sd_s <- 2.17
n_s<- 62
error_s <- qnorm(0.975)*sd_s/sqrt(n_s)
left <- u_sp-error_s
right <- u_sp+error_s
left
error_s
right
#OP
u_op <- 8.78
sd <- 2.65
n <- 51
error <- qnorm(0.975)*sd/sqrt(n)
left <- u_op-error
right <- u_op+error
left
error
right 

############
###Seed Set Models
plt_seed<- seed_land %>%mutate( trmnt = fct_recode(trmnt, 
                          "Supplemental Pollination" = "hp", "Open Pollination"="op"))%>%
  group_by(site,prop.c,trmnt,ID,round)%>%#calculate fruit removed
  mutate(frt_removed_rd=sum(frt_removed)/n(), trt_flw_rd=sum(flowers)/n())%>%
  ungroup()%>%
  group_by(site,prop.c,trmnt,ID)%>%
  summarise(frt=n(),trt_flw=sum(trt_flw_rd),tot_frt_removed=sum(frt_removed_rd),
            seeds=sum(total_seeds),
            mean_seed=mean(total_seeds),
            seeds_frt=seeds/frt,mean_sf=mean(seeds/frt))



#Post-hoc t.test for difference in mean untreated frt remvoal between treatment

t.test(tot_frt_removed~trmnt,data=plt_seed)
t.test(trt_flw~trmnt,data=plt_seed)
#

############
#   models #
###########

####
###Models 1&2:
#mean seed per fruit response

#model 1:
#using glmmTMB
seed_mean_mod_TMB<-glmmTMB(mean_seed~trmnt*prop.c +(1|site), data=plt_seed) 
summary(seed_mean_mod_TMB)#reported in table 1

#using lme4
seed_mean_mod<-lmer(mean_seed~trmnt*prop.c +(1|site), data=plt_seed)

summary(seed_mean_mod)


###model 2:
#no interaction term
seed_mean_mod2_TMB<-glmmTMB(mean_seed~trmnt+prop.c +(1|site), data=plt_seed) 
summary(seed_mean_mod2_TMB)#reported in table 2

seed_mean_mod2<-lmer(mean_seed~trmnt+prop.c +(1|site), data=plt_seed) 
summary(seed_mean_mod2)



### Confidence intervals


seed_mod_CI_TMB_wald<- confint(seed_mean_mod_TMB)
seed_mod_CI_TMB_wald # used in table 1

seed_mod_CI2<-confint(seed_mean_mod2_TMB)
seed_mod_CI2# used in table 2


# Model diagnostics using DHARMa

# Diagnostiscts on model 1: 
sim_mod1<-simulateResiduals(fittedModel = seed_mean_mod_TMB, n = 250)# warning message because glmmTMB was used
plot(sim_mod1) # no strange patterns in predicted vs residuals
testDispersion(sim_mod1)
testUniformity(sim_mod1)


#compare variation among trmnts and sites
#control for random effect
new_plt<-plt_seed
# set ID and site to NA in new DF
new_plt$ID=NA
new_plt$site=NA

#calculate new predicted values that average over the REs  
pred = predict(seed_mean_mod_TMB, newdata = new_plt)
#
plotResiduals(pred, sim_mod1$scaledResiduals) #doesn't differ from direct approach

#plot against other predictors
par(mfrow=c(1,2))
plotResiduals(plt_seed$prop.c, sim_mod1$scaledResiduals)
plotResiduals(plt_seed$trmnt, sim_mod1$scaledResiduals)
par(mfrow=c(1,1))
## looks good!

# Diagnostiscts on model 2: 
sim_mod2<-simulateResiduals(fittedModel = seed_mean_mod2_TMB, n = 250)# warning message because glmmTMB was used
plot(sim_mod2) # no strange patterns in predicted vs residuals
testDispersion(sim_mod2)
testUniformity(sim_mod2)


#compare variation among trmnts and sites
#control for re
new_plt<-plt_seed
# set ID and site to NA in new DF
new_plt$ID=NA
new_plt$site=NA

#calculate new predicted values that average over the REs  
pred = predict(seed_mean_mod2_TMB, newdata = new_plt)
#
plotResiduals(pred, sim_mod2$scaledResiduals) #doesn't differ from direct approach
#plot against other predictors
par(mfrow=c(1,2))
plotResiduals(plt_seed$prop.c, sim_mod2$scaledResiduals)
plotResiduals(plt_seed$trmnt, sim_mod2$scaledResiduals)
par(mfrow=c(1,1))
## looks good!


### 
#Supplemental models

#model s1a
#Frt matured as a response
frt_mod<-glmmTMB(frt~trmnt*prop.c+(1|site), data=plt_seed,family=nbinom2)
#output used in supplemental table 1a
summary(frt_mod)

frt_mod_CI<- confint(frt_mod)
frt_mod_CI

#model s1b
#no interaction frt model
frt_mod2<-glmmTMB(frt~trmnt+prop.c+(1|site), data=plt_seed,family=nbinom2)
#output used in supplemental table 1b
summary(frt_mod2)
#frt unimpacted by prop.c or trmtn

frt_mod2_CI<- confint(frt_mod2)
frt_mod2_CI

# Diagnostiscts on model s1a: 
sim_frt<-simulateResiduals(fittedModel = frt_mod, n = 250)# warning message because glmmTMB was used
plot(sim_frt) # no strange patterns in predicted vs residuals
testDispersion(sim_frt)
testUniformity(sim_frt)


#compare variation among trmnts and sites
#control for re
new_plt<-plt_seed
# set ID and site to NA in new DF
new_plt$ID=NA
new_plt$site=NA

#calculate new predicted values that average over the REs  
pred = predict(frt_mod, newdata = new_plt)
#
plotResiduals(pred, sim_frt$scaledResiduals) #doesn't differ from direct approach
#plot against other predictors
par(mfrow=c(1,2))
plotResiduals(plt_seed$prop.c, sim_frt$scaledResiduals)
plotResiduals(plt_seed$trmnt, sim_frt$scaledResiduals)
par(mfrow=c(1,1))
## frt mod looks okay! Modify the above code for the second model for diagnostics


##models s2a&b

#Seed mod, no offset
seed_mod<-glmmTMB(seeds~trmnt*prop.c+(1|site), data=plt_seed,family=nbinom2)
summary(seed_mod)#data for table s2a


seed_mod_CI<- confint(seed_mod)
seed_mod_CI

seed_mod2<-glmmTMB(seeds~trmnt+prop.c+(1|site), data=plt_seed,family=nbinom2)
summary(seed_mod2)#data for table s1a


seed_mod2_CI<- confint(seed_mod2)
seed_mod2_CI

####diagnostics on models s2a&b
sim_seed2<-simulateResiduals(fittedModel = seed_mod2, n = 250)# warning message because glmmTMB was used
plot(sim_seed2) # no strange patterns in predicted vs residuals
testDispersion(sim_seed2)
testUniformity(sim_seed2)

#now for no interaction model
sim_seed<-simulateResiduals(fittedModel = seed_mod, n = 250)
plot(sim_seed)
testDispersion(sim_seed)#both this and interaction model have similar results, 
#Few other tests: uniformity 
testUniformity(sim_seed)


#################
###Bee abundance
bee_abun<-read.csv("bee_abun.csv")

#model 5- bee abundance

abun_mod<-lm(bee_abun~prop.c,data=bee_abun)
summary(abun_mod)
confint(abun_mod)
qqPlot(residuals(abun_mod))
plot(abun_mod)

# total bees
sum(bee_abun$bee_abun)
################################################################ 
#95Confidence interval for mean
mean(bee_abun$bee_abun)
sd(bee_abun$bee_abun)
u_b<- 458.14
sd_b <- 130.32
n_b<- 7
error_b<- qnorm(0.975)*sd_b/sqrt(n_b)
left <- u_b-error_b
right <- u_b+error_b
left
error_b
right




############################
#Fig 2: LMM predicted response vs raw data

##########################

####note!!!: model object is generated in model script! Run it first!

#########################
#Plotting predicted effect and actual data
################

theme_set(theme_bw())
theme_set(theme_classic())
#run get_model_data to extract ggplot usable output
s<-get_model_data(seed_mean_mod_TMB,ci.lvl= .95, type="pred",terms=c("prop.c","trmnt"), 
                  pred.type="re") 
s<-s%>%mutate(group = fct_recode(group, 
                                 "Supplemental Pollination" = "SP","Open Pollination"="OP"))
### make separate dataframes for CI for plot
shp<-s%>%filter(group=="Supplemental Pollination")
sop<-s%>%filter(group=="Open Pollination")

###make plot for seed model
ggplot(data=s, aes(x, predicted),linetype=group)+
  geom_line(aes(linetype=group),size=2)+
  geom_point(data=plt_nr,aes(prop.c, (seeds/frt),shape=trmnt),position="jitter", inherit.aes = FALSE,size=5)+
  geom_ribbon(data=shp,aes(x=x,ymin=conf.low, ymax=conf.high),alpha = 0.3,inherit.aes = FALSE)+
  geom_ribbon(data=sop,aes(x=x,ymin=conf.low, ymax=conf.high),alpha = 0.3,inherit.aes = FALSE)+
  labs(shape="Treatment",linetype="Predicted Response")+
  xlab("% Agriculture Surrounding Restoration")+ylab("Mean Seeds per Treated Fruit")+
  ggtitle("Fig. 1: Effect of Surrounding Agriculture on Seed Set",
          subtitle=" LMM Predicted Response vs. Data")+
  theme(axis.text=element_text(size=20,color="black"),
        axis.title=element_text(size=24),
        title=element_text(size=26),
        legend.text=element_text(size=24),
        axis.line = element_line(colour = "black",size=1),panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),legend.position = "none")


###
###Fig 2: 
#mean seeds produced per fruits collected for a plant 
seed_land%>% 
  mutate(trmnt = fct_recode(trmnt, "Supplemental Pollination" = "hp", "Open Pollination"="op"))%>%
  group_by(trmnt, ID)%>%
  summarize(mean_seed=mean(total_seeds))%>%
  ggplot(aes(trmnt ,mean_seed))+geom_boxplot(lwd=1.5,outlier.size=3,color="black")+ theme_bw()+
  labs(x="Treatment" ,y="Mean Seeds per Treated Fruit")+
  ggtitle("Fig. 2: Mean Seed Set per Treatment")+
  theme(axis.text=element_text(size=24,color="black"),
        axis.title=element_text(size=26),title=element_text(size=26),
        axis.line = element_line(colour = "black", size=2),panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),legend.position = "none")

#plot 3: bee abundance
theme_set(theme_bw())
theme_set(theme_classic())
#run get_model_data to extract ggplot usable output
b<-get_model_data(abun_mod,ci.lvl= .95, type="pred",terms=c("prop.c")) 
confint(abun_mod)

###make plot for bee model
ggplot(data=b, aes(x, predicted),linetype=group)+
  geom_line(aes(linetype=group),size=2)+
  geom_point(data=abun_sites,aes(prop.c, sum_abun), inherit.aes = FALSE,size=6)+
  geom_ribbon(data=b,aes(x=x,ymin=conf.low, ymax=conf.high),alpha = 0.3,inherit.aes = FALSE)+
  labs(linetype="Predicted Response")+
  xlab("% Agriculture Surrounding Restoration")+ylab("Bee Abundance")+
  ggtitle("Fig. 3: Effect of Surrounding Agriculture on Bee Abundance",
          subtitle=" LM Predicted Response vs. Data")+
  theme(axis.text=element_text(size=20,color="black"),
        axis.title=element_text(size=24),
        title=element_text(size=26),
        legend.text=element_text(size=24),
        axis.line = element_line(colour = "black",size=1),panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),legend.position = "none")