## Set working directory first

# Clear plots
if(!is.null(dev.list())) dev.off()
# Clear console
cat("\014") 
# Clean workspace
rm(list=ls())
# Set working directory
setwd("C:/Users/")


#load packages
library(ggplot2)
library(dplyr)
library(glmmTMB)
library(pscl)
library(DHARMa)
library(emmeans)
library(bbmle)
library(ggeffects)
library(lmerTest)
library(ggpubr)
library(AICcmodavg)
library(see)
library(performance)
library(data.table)
library(lazyWeave)
library(tidyr)

################################################################################
################################################################################
##################                                         #####################  
##################            total bat eaten              #####################
##################                                         #####################
################################################################################
################################################################################

## attach data 
total_orchard <- read.csv("frugivore.foraging.updated.csv",na.strings=c("", "NA"),header=T,sep=",") # read text file
bat_backyard <- read.csv("backyard.foraging.pattern.csv",na.strings=c("", "NA"),header=T,sep=",") 
#total_orchard <- total_orchard[, -c(8)] #remove 'X' column
total_orchard <- total_orchard[, -c(2)] #remove 'date' column
#total_orchard <- total_orchard[-c(959),] #remove duplicate row (two bird sums for that tree/cycle)
attach(total_orchard)

#explore
total_orchard %>% group_by(site, damage.type) %>% summarise(total = sum(total.damage.type))
#only bats and birds are relevant, rest barely ate anything
#make separate datasets for bats and birds as birds were not recorded in backyards
#for bat dataset, combine with orchards. Have to remove and adjust a few columns
bat_orchard <- total_orchard[total_orchard$damage.type %in% c("Bat"),]
bat_orchard <- bat_orchard[bat_orchard$time %in% c("18:00 - 00:00", "00:00 - 06:00"),]
bat_orchard$total.bat <- bat_orchard$total.damage.type 
bat_orchard <- bat_orchard[, -c(5:6)] #remove 'damage.type' and 'total.damage.type'
#adjust bat_backyard and combine with orchards
bat_backyard <- bat_backyard[, -c(3)] #remove timeslot
bat_orchard <- rbind(bat_orchard, bat_backyard) #combine datasets

#add observation-level random effect (ID)
bat_orchard$ID <- c(1:682)

attach(bat_orchard)

#we are dealing with a count variable, check if distribution is poisson shaped 
hist(total.bat) #poisson shaped

#see how bat damage varies among weeks and time slots between sites
ggplot(bat_orchard, aes(factor(cycle), total.bat, fill = factor(time))) +
  geom_boxplot(outlier.colour="black", outlier.shape=16, outlier.size=1.5, notch=FALSE) +
  labs(fill = "time", y = "total eaten (n)", x = "cycle") +
  facet_wrap( ~ site)
#seems like total bat eaten follows a similar trend over the weeks in every timeslot,
#adding interaction to model may not be necessary

#variation in response among trees for each site 
ggplot(bat_orchard, aes(factor(tree), total.bat, fill = factor(site))) +
  geom_boxplot(outlier.colour="black", outlier.shape=16, outlier.size=1.5, notch=FALSE) +
  labs(fill = "site", y = "total eaten (n)", x = "tree ID") 
#based on variation in mean total bat eaten among trees, it seems appropriate to add tree
#as random effect

#test whether random effect actually improves the regression model we want to fit
#fit models with REML to compare random effect structures or model with and without
#we are testing on the boundary, p-values are a bit unreliable.
#advice is to divide p value by 2 when comparing model with and without random effect ()
#model without random effect
GLMMbat1 <- glmmTMB(total.bat ~ cycle*site + site*time, data = bat_orchard, family = poisson, REML = TRUE)
#now model with random intercept
GLMMbat2 <- glmmTMB(total.bat ~ cycle*site + time*site + (1|tree), data = bat_orchard, family = poisson, REML = TRUE)
anova(GLMMbat1,GLMMbat2, test = "Chisq") #random effect improves model
summary(GLMMbat2)
#model also manages to estimate variances among all levels, perfect!

#fit global model without REML and check if not overparameterized
GLMMbat2 <- glmmTMB(total.bat ~ cycle*site + time*site + (1|tree), data = bat_orchard, family = poisson)
AICcmodavg::AICc(GLMMbat2, return.K = TRUE) # k = 10 on n = 682, so n/k = ~68 = perfect!

#now lets validate the global model
simulationOutput <- simulateResiduals(fittedModel = GLMMbat2, n = 1000) 
testDispersion(simulationOutput) #looks good
plotSimulatedResiduals(simulationOutput = simulationOutput) 
testZeroInflation(simulationOutput) #indicates the observed data has more zeros than expected

#slightly zero-inflated (maybe caused by overdispersion), try nbinom2
GLMMbat3 <- glmmTMB(total.bat ~ cycle*site + time*site + (1|tree), data = bat_orchard, family = nbinom2)
AICctab(GLMMbat3, GLMMbat2, weights=TRUE,delta=TRUE) #big improvement
simulationOutput <- simulateResiduals(fittedModel = GLMMbat3, n = 1000) 
plotSimulatedResiduals(simulationOutput = simulationOutput) #still some deviations
testZeroInflation(simulationOutput) #no longer zero-inflated, but still some dispersion issues

#try fitting zero-inflation parameter
GLMMbat4 <- glmmTMB(total.bat ~ cycle*site + time*site + (1|tree), ziformula = ~1, data = bat_orchard, family = nbinom2)
GLMMbat5 <- glmmTMB(total.bat ~ cycle*site + time*site + (1|tree), ziformula = ~cycle*time, data = bat_orchard, family = nbinom2)
AICctab(GLMMbat4, GLMMbat3, GLMMbat5, weights=TRUE,delta=TRUE) 
#zero-inflation parameter depending on fixed effects improves model greatly
#validate again
simulationOutput <- simulateResiduals(fittedModel = GLMMbat5, n = 1000) 
testDispersion(simulationOutput) #good
plotSimulatedResiduals(simulationOutput = simulationOutput) #perfect
testZeroInflation(simulationOutput) #perfect

if (require("see")) {
  check_predictions(GLMMbat5)  #looks great
}

#check if global model is better than null model (variation in data not explained by just random variation)
GLMMbatnull <- glmmTMB(total.bat ~ 1 + (1|tree), data = bat_orchard, family = nbinom2)
anova(GLMMbatnull, GLMMbat5, test = "Chisq")

#now inference with GLMMbat5
summary(GLMMbat5)
confint(emmeans(GLMMbat5, pairwise ~ time|site))
  
#put test() in front of emtrends to get p value for trends

#visualize trend in number of bat eaten fruits over the cycles
pred1 <- ggpredict(GLMMbat5, c("cycle", "site"), type = 're')
pred1$site <- pred1$group #this line is necessary to get a nice legend

#test plot
ggplot() + geom_smooth(aes(x, predicted, ymin = conf.low, ymax = conf.high, color =site, shape = site), data = pred1, se = FALSE) +
  geom_point(aes(cycle, total.bat, group = site, color = site, shape = site), data = bat_orchard, size = 2, inherit.aes = FALSE) +
  facet_wrap(~site)

#now nice plot, also relevel site to improve it for ggarrange later
bat_orchard$site <- factor(bat_orchard$site, levels=c('Beaux Songes', 'Calebasses', 'Backyard'))
p1 <- ggplot() + geom_point(aes(cycle, total.bat, group = site, color = site), data = bat_orchard, size = 1.5, inherit.aes = FALSE, position = position_jitterdodge(dodge.width = -0.7, jitter.width = 0.1)) +
  geom_smooth(aes(x, predicted, ymin = conf.low, ymax = conf.high, fill = site, shape = site), data = pred1, color = "black") +
  facet_wrap(. ~ site)+
  scale_color_manual(values = c("grey55", "grey55", "grey55")) +
  xlab("cycle") + ylab("fruits eaten by bats (n)") + 
  theme_linedraw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=11),
        legend.text=element_text(size= 11),
        legend.title = element_blank())+
  guides(color = guide_legend(override.aes = list(fill = NA)),
         linetype = guide_legend(override.aes = list(fill = NA))) +
  theme(panel.background = element_rect(fill = "white"),
        strip.background = element_rect(colour="white",fill="white"),
        strip.text = element_text(colour = 'black', size = 11),
        legend.position = 'none') +
  theme(panel.spacing = unit(1, "lines")) #for better spacing in between facets
  
  
#visualize difference between time slots in every site
pred2 <- ggpredict(GLMMbat5, c("site", "time"), type = 're')
pred2$time <- pred2$group #this line is necessary to get a nice legend

p2 <- ggplot() + geom_crossbar(aes(x, predicted, ymin = conf.low, ymax = conf.high, color = group, shape = group), data = pred2, width = 0.55, position = position_dodge(width = 0.7)) +
  geom_point(aes(site, total.bat, group = time, color = time, shape = time), data = bat_orchard, size = 2, inherit.aes = FALSE, position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.1)) +
  scale_color_manual(values = c("black", "grey55")) +
  xlab("") + ylab("total bat eaten (n)") + 
  theme_linedraw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"), 
        axis.text=element_text(size=11),
        axis.title=element_text(size=12),
        legend.text=element_text(size= 11),
        legend.title = element_blank()) +
  theme(panel.background = element_rect(fill = "white"), 
        strip.background = element_rect(colour="white",fill="white"),
        strip.text = element_text(colour = 'black', size = 11))

#alternative plot for ggarrange later  
bat_orchard$time <- factor(bat_orchard$time, levels=c('18:00 - 00:00', '00:00 - 06:00'))
pred22 <- ggpredict(GLMMbat5, c("time", "site"), type = 're')
pred22$site <- pred22$group #this line is necessary to get a nice legend
  
p22 <- ggplot() + geom_point(aes(time, total.bat, group = site, color = site), data = bat_orchard, size = 2, inherit.aes = FALSE, position = position_jitterdodge(dodge.width = -0.7, jitter.width = 0.1)) +
  geom_crossbar(aes(x, predicted, ymin = conf.low, ymax = conf.high), color = 'black', data = pred22, width = 0.6, position = position_dodge(width = 0.7)) +
  scale_color_manual(values = c("grey55", "grey55", "grey55")) +
  facet_wrap(~ site) +
  xlab("") + ylab("fruit eaten by bats (n)")+ 
  theme_linedraw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=11),
        legend.text=element_text(size= 11),
        legend.title = element_blank()) +
 theme(panel.background = element_rect(fill = "white"),
        legend.position = "none", 
        strip.background = element_rect(colour="white",fill="white"),
        strip.text = element_text(colour = 'black', size = 11))
 
############################################
########                            ########
########       now only birds       ########
########                            ########
############################################
bird_orchard <- total_orchard[total_orchard$damage.type %in% c("Bird"),]
attach(bird_orchard)
  
GLMMbird1 <- glmmTMB(total.damage.type ~ cycle*site + time*site, data = bird_orchard, family = poisson, REML = TRUE)
#now model with random intercept
GLMMbird2 <- glmmTMB(total.damage.type ~ cycle*site + site*time + (1|tree), data = bird_orchard, family = poisson, REML = TRUE)
GLMMbird3 <- glmmTMB(total.damage.type ~ cycle*site + time*site + (1 + cycle|tree), data = bird_orchard, family = poisson, REML = TRUE)
anova(GLMMbird1,GLMMbird3, test = "Chisq") #random effect improves model
anova(GLMMbird3,GLMMbird2, test = "Chisq") #random slopes model is best
summary(GLMMbird2)#mtrouble estimating variances, go with simpler model!
#fit global model without REML and check if not overparameterized
GLMMbird2 <- glmmTMB(total.damage.type ~ cycle*site + time*site + (1|tree), data = bird_orchard, family = poisson)
AICcmodavg::AICc(GLMMbird2, return.K = TRUE) # k = 11 on n = 1093, more than enough

#fit negative binomial straight away and compare
GLMMbird3 <- glmmTMB(total.damage.type ~ cycle*site + time*site + (1|tree), data = bird_orchard, family = nbinom2)
AICctab(GLMMbird3, GLMMbird2, weights=TRUE,delta=TRUE) #big improvement
simulationOutput <- simulateResiduals(fittedModel = GLMMbird3, n = 1000) 
plotSimulatedResiduals(simulationOutput = simulationOutput) #perfect
testZeroInflation(simulationOutput) #perfect

#check if global model is better than null model (variation in data not explained by just random variation)
GLMMbirdnull <- glmmTMB(total.damage.type ~ 1 + (1|tree), data = bird_orchard, family = nbinom2)
anova(GLMMbirdnull, GLMMbird3, test = "Chisq") #better than null model
  
#now inference
confint(emmeans(GLMMbird3, pairwise ~ time|site))
emtrends(GLMMbird3, pairwise ~ site, var = "cycle") 

#visualize trend in number of bird eaten fruits over the cycles
predbird <- ggpredict(GLMMbird3, c("cycle", "site"), type = 're')
predbird$site <- predbird$group #this line is necessary to get a nice legend

#test plot
ggplot() + geom_smooth(aes(x, predicted, ymin = conf.low, ymax = conf.high, color =site, shape = site), data = predbird, se = FALSE) +
  geom_point(aes(cycle, total.damage.type, group = site, color = site, shape = site), data = bird_orchard, size = 2, inherit.aes = FALSE) +
  facet_wrap(~site)

pbird <- ggplot() + geom_point(aes(cycle, total.damage.type, group = site, color = site), data = bird_orchard, size = 1.5, inherit.aes = FALSE, position = position_jitterdodge(dodge.width = -0.7, jitter.width = 0.1)) +
  geom_smooth(aes(x, predicted, ymin = conf.low, ymax = conf.high, fill = site, shape = site), data = predbird, color = "black") +
  facet_wrap(. ~ site)+
  scale_color_manual(values = c("grey55", "grey55", "grey55")) +
  xlab("cycle") + ylab("fruits eaten by birds (n)") + 
  theme_linedraw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=11),
        legend.text=element_text(size= 11),
        legend.title = element_blank())+
  guides(color = guide_legend(override.aes = list(fill = NA)),
         linetype = guide_legend(override.aes = list(fill = NA))) +
  theme(panel.background = element_rect(fill = "white"),
        strip.background = element_rect(colour="white",fill="white"),
        strip.text = element_text(colour = 'black', size = 11),
        legend.position = 'none') +
  theme(panel.spacing = unit(1, "lines")) #for better spacing in between facets

#visualize difference between time slots in every site
predbird2 <- ggpredict(GLMMbird3, c("time", "site"), type = 're')
predbird2$site <- predbird2$group #this line is necessary to get a nice legend

#also relevel timeslots
bird_orchard$time <- factor(bird_orchard$time, levels=c('06:00 - 12:00', '12:00 - 18:00', '18:00 - 00:00', '00:00 - 06:00'))
#now you have to run the model GLMMbird31 + predbird2 again to make it work
pbird2 <- ggplot() + geom_point(aes(time, total.damage.type, group = site, color = site), data = bird_orchard, size = 2, inherit.aes = FALSE, position = position_jitterdodge(dodge.width = -0.7, jitter.width = 0.1)) +
  geom_crossbar(aes(x, predicted, ymin = conf.low, ymax = conf.high), color = 'black', data = predbird2, width = 0.6, position = position_dodge(width = 0.7)) +
  scale_color_manual(values = c("grey55", "grey55", "grey55", "grey55")) +
  facet_wrap(~ site) +
  xlab("") + ylab("fruit eaten by birds (n)") + scale_y_continuous(limits = c(0, 40)) +
  theme_linedraw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=11),
        legend.text=element_text(size= 11),
        legend.title = element_blank()) +
theme(panel.background = element_rect(fill = "white"),
      legend.position = "none", 
      strip.background = element_rect(colour="white",fill="white"),
      strip.text = element_text(colour = 'black', size = 11)) 

#alternative plot difference between time slots in every site
#first relevel factor
bird_orchard$time <- factor(bird_orchard$time, levels=c('06:00 - 12:00', '12:00 - 18:00', '18:00 - 00:00', '00:00 - 06:00'))
#NOW RERUN GLMMbird3 BEFORE CONTINUING!
predbird22 <- ggpredict(GLMMbird3, c("site", "time"), type = 're')
predbird22$time <- predbird22$group #this line is necessary to get a nice legend



pbird22 <- ggplot() + geom_crossbar(aes(x, predicted, ymin = conf.low, ymax = conf.high, color = group, shape = group), data = predbird22, width = 0.75, position = position_dodge(width = 0.9)) +
  geom_point(aes(site, total.damage.type, group = time, color = time, shape = time), data = bird_orchard, size = 2, inherit.aes = FALSE, position = position_jitterdodge(dodge.width = 0.9, jitter.width = 0.1)) +
  scale_color_manual(values = c("black", "grey30", "grey60", "grey80")) +
  xlab("") + ylab("total bird eaten (n)") + scale_y_continuous(limits = c(0, 40)) +
  theme_linedraw()+ 
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"), 
        axis.text=element_text(size=11),
        axis.title=element_text(size=12),
        legend.text=element_text(size= 11),
        legend.title = element_blank()) +
  theme(panel.background = element_rect(fill = "white"), 
        strip.background = element_rect(colour="white",fill="white"),
        strip.text = element_text(colour = 'black', size = 11))

#arrange the plots
ggarrange(p1, ggarrange(pbird, ncol = 2, widths = c(1.1, 0.5)), nrow = 2)

ggarrange(p2, pbird22, nrow = 2)

ggarrange(ggarrange(p22, ncol = 2, widths = c(2, 0.5)), pbird2, nrow = 2)




################################################################################
################################################################################
##################                                         #####################  
##################    fruit ripeness proportion eaten      #####################
##################                                         #####################
################################################################################
################################################################################
  
## attach data 
bat_orchard_prop <- read.csv("amount.eaten.bat.bird.csv",na.strings=c("", "NA"),header=T,sep=",") # read text file
#add observation-level random effect (ID)
bat_orchard_prop$ID <- c(1:1563)
#add fruit ripeness as factor (we need this for the plots later)
bat_orchard_prop$factorripe <- as.factor(bat_orchard_prop$fruit.ripeness)
#it only has one observation and fucks with analysis
attach(bat_orchard_prop)
#our response variable is a continuous proportion (not arising from counts)
#usually you would model this with a beta error distribution
#check distribution of response variable to see if appropriate
hist(amount.eaten) #clumping of values at lower and less at upper end of the scale
#looks like beta regression is appropriate (Kubinec 2022)

#check whether including both cycle and fruit ripeness in model would lead to
#to problems with autocorrelation i.o.w. VIF values > 5
Z <- cbind(bat_orchard_prop$cycle, bat_orchard_prop$fruit.ripeness)
colnames(Z) <- c("cycle", "fruit.ripeness")
corvif(Z) #no problems with VIF values 

#see how proportion eaten varies among ripe/unripe, weeks and time slots between sites
ggplot(bat_orchard_prop, aes(factor(cycle), amount.eaten, fill = factor(time))) +
  geom_boxplot(outlier.colour="black", outlier.shape=16, outlier.size=1.5, notch=FALSE) +
  labs(fill = "time", y = "proportion flesh eaten", x = "cycle") +
  facet_wrap( ~ orchard + damage.type)
#seems like proportion flesh eaten by bats follows a similar trend over the weeks in every timeslot,
#adding interaction to model may not be necessary
#however, for bird and parakeet the trend seems different between timeslots
d1<-bat_orchard_prop %>% group_by(damage.type, orchard, time) %>% summarise(mean=mean(amount.eaten))
d1 #also means in flesh eaten is the same before and after midnight for bats
#for birds there are some differences between the 06:00 - 12:00 and 12:00 - 18:00
#maybe include timeslot in model?
bat_orchard_prop$count <- rep(1, 1563)
bat_orchard_prop %>% group_by(damage.type) %>% summarise(sum=sum(count))

#lets see how response varies among trees for each site 
ggplot(bat_orchard_prop, aes(factor(tree), amount.eaten, fill = factor(orchard))) +
  geom_boxplot(outlier.colour="black", outlier.shape=16, outlier.size=1.5, notch=FALSE) +
  labs(fill = "orchard", y = "proportion flesh eaten", x = "tree ID") +
  facet_wrap( ~ damage.type)
#based on variation in mean proportion flesh eaten among trees, it there is minimal
#variation in mean proportion of flesh eaten among trees. Test with REML to see whether
#random effect is appropriate

#first choose an error distribution. As we are dealing with a non-count proportion 
#(fraction of a fruit that was eaten) we need to employ a beta error distribution
#the problem with the beta-distribution is that it does not allow zeros and ones,
#therefore you would usually run a zero-one augmented beta GLMM, but this is tedious
#nowadays, a nice alternative exists -> the ordered beta distribution
#this provides better results than zero-one augmented models (Kubinec 2022) and 
#is also easy to fit with glmmTMB
#now test whether random effect is appropriate
GLMMprop1 <- glmmTMB(amount.eaten ~ orchard*fruit.ripeness*damage.type + cycle*orchard*damage.type, data = bat_orchard_prop, family = ordbeta, REML = TRUE)
#now model with random intetercept
GLMMprop2 <- glmmTMB(amount.eaten ~ orchard*fruit.ripeness*damage.type + cycle*orchard*damage.type + (1|tree), data = bat_orchard_prop, family = ordbeta, REML = TRUE)
GLMMprop3 <- glmmTMB(amount.eaten ~ orchard*fruit.ripeness*damage.type + cycle*orchard*damage.type + (1 + fruit.ripeness|tree), data = bat_orchard_prop, family = ordbeta, REML = TRUE)
GLMMprop4 <- glmmTMB(amount.eaten ~ orchard*fruit.ripeness*damage.type + cycle*orchard*damage.type + (1 + cycle|tree), data = bat_orchard_prop, family = ordbeta, REML = TRUE)
anova(GLMMprop1,GLMMprop2, test = "Chisq") #random effect improves model
anova(GLMMprop2,GLMMprop3, test = "Chisq") #random slope model of fruit ripeness improves model
anova(GLMMprop2,GLMMprop4, test = "Chisq") #random slope model of cycle doesnt
#cannot leave out random effect! random intercept + slope model is best
#however, model with random slopes cannot compute variances: a lot are very close to 0
#use random intercept model instead
summary(GLMMprop2) #no difficulties estimating the variances with random intercept model

#fit global model without REML and check if not overparameterized
GLMMprop2 <- glmmTMB(amount.eaten ~ orchard*fruit.ripeness*damage.type + cycle*orchard*damage.type + (1|tree), data = bat_orchard_prop, family = ordbeta)
AICcmodavg::AICc(GLMMprop2, return.K = TRUE) # k = 22 on n = 1564, so n/k = 71 = perfect

#now lets validate the global model
simulationOutput <- simulateResiduals(fittedModel = GLMMprop2, n = 1000) 
testDispersion(simulationOutput) #overdispersed
plotSimulatedResiduals(simulationOutput = simulationOutput) 
testZeroInflation(simulationOutput)
#model slightly overdispersed, might help to specify zero-inflation parameter

#specify zero-inflation parameter
GLMMprop3 <- glmmTMB(amount.eaten ~ orchard*fruit.ripeness*damage.type + cycle*damage.type*orchard + (1|tree), data = bat_orchard_prop, ziformula = ~orchard + damage.type + (1|tree), family = ordbeta)
simulationOutput <- simulateResiduals(fittedModel = GLMMprop3, n = 1000) 
testDispersion(simulationOutput) #good
plotSimulatedResiduals(simulationOutput = simulationOutput) #has improved a little bit

if (require("see")) {
  check_predictions(GLMMprop3)  #model underestimates 0.25 and 0.75 values, while overestimating 0.5 values
}

#check if global model is better than null model (variation in data not explained by just random variation)
GLMMpropnull <- glmmTMB(amount.eaten ~ 1 + (1|tree), ziformula = ~1, dispformula = ~1, data = bat_orchard_prop, family = ordbeta)
anova(GLMMpropnull, GLMMprop3, test = "Chisq") #better than null model

#now inference
summary(GLMMprop3)
emtrends(GLMMprop3, pairwise ~ orchard*damage.type, var = "fruit.ripeness")
#here it shows that generally higher fruit ripeness leads to higher proportion eaten by bats,
#this effect is highly significant within all orchards
emtrends(GLMMprop3, pairwise ~ orchard*damage.type, var = "cycle")
#here it shows that bats eat larger proportions in Beaux Songes as cycle increases
mylist <- list(cycle=c((mean(cycle)-sd(cycle)),mean(cycle),(mean(cycle)+sd(cycle))), orchard=c("Beaux Songes", "Calebasses"), damage.type=c("Bat","Parakeet","Other birds"), fruit.ripeness = c(0, 0.25, 0.5, 0.75, 1))
confint(emmeans(GLMMprop3, pairwise ~ damage.type|orchard|cycle, at = mylist))
#here it shows that bats eat more than parakeets in one site, and parakeets less
#than other birds in all sites 
confint(emmeans(GLMMprop3, pairwise ~ damage.type|orchard|fruit.ripeness, at = mylist))

#visualize trend in proportion of fruit eaten for the ripeness gradient
pred4 <- ggpredict(GLMMprop3, c("fruit.ripeness", "damage.type", "orchard"), type = 're')
pred4$damage.type <- pred4$group #this line is necessary to get a nice legend
pred4$orchard <- pred4$facet #this line is necessary to get a nice legend

#test plot
ggplot() + geom_boxplot(aes(fruit.ripeness, amount.eaten, group = fruit.ripeness), data = bat_orchard_prop, size = 0.5, inherit.aes = FALSE) +
  geom_smooth(aes(x, predicted, ymin = conf.low, ymax = conf.high, color =damage.type, shape = damage.type), data = pred4, se = FALSE) +
  facet_wrap(~orchard+damage.type)

#now nice plot
bat_orchard_prop$damage.type <- factor(bat_orchard_prop$damage.type, levels=c('Bat', 'Parakeet', 'Other birds'))
bat_orchard_prop$yloc = max(bat_orchard_prop$amount.eaten) + .1
bat_orchard_prop <- bat_orchard_prop %>% 
  mutate(obslabel = case_when(orchard == 'Beaux Songes' & damage.type == 'Bat' & fruit.ripeness == 0 ~ '171',
                              orchard == 'Beaux Songes' & damage.type == 'Bat' & fruit.ripeness == 0.25 ~ '417',
                              orchard == 'Beaux Songes' & damage.type == 'Bat' & fruit.ripeness == 0.5 ~ '99',
                              orchard == 'Beaux Songes' & damage.type == 'Bat' & fruit.ripeness == 0.75 ~ '50',
                              orchard == 'Beaux Songes' & damage.type == 'Bat' & fruit.ripeness == 1 ~ '18',
                              orchard == 'Beaux Songes' & damage.type == 'Parakeet' & fruit.ripeness == 0 ~ '5',
                              orchard == 'Beaux Songes' & damage.type == 'Parakeet' & fruit.ripeness == 0.25 ~ '17', 
                              orchard == 'Beaux Songes' & damage.type == 'Parakeet' & fruit.ripeness == 0.5 ~ '8', 
                              orchard == 'Beaux Songes' & damage.type == 'Parakeet' & fruit.ripeness == 0.75 ~ '5',
                              orchard == 'Beaux Songes' & damage.type == 'Parakeet' & fruit.ripeness == 1 ~ '1',
                              orchard == 'Beaux Songes' & damage.type == 'Other birds' & fruit.ripeness == 0 ~ '1',
                              orchard == 'Beaux Songes' & damage.type == 'Other birds' & fruit.ripeness == 0.25 ~ '19',
                              orchard == 'Beaux Songes' & damage.type == 'Other birds' & fruit.ripeness == 0.5 ~ '17',
                              orchard == 'Beaux Songes' & damage.type == 'Other birds' & fruit.ripeness == 0.75 ~ '11',
                              orchard == 'Beaux Songes' & damage.type == 'Other birds' & fruit.ripeness == 1 ~ '5',
                              orchard == 'Calebasses' & damage.type == 'Bat' & fruit.ripeness == 0 ~ '160',
                              orchard == 'Calebasses' & damage.type == 'Bat' & fruit.ripeness == 0.25 ~ '282',
                              orchard == 'Calebasses' & damage.type == 'Bat' & fruit.ripeness == 0.5 ~ '109',
                              orchard == 'Calebasses' & damage.type == 'Bat' & fruit.ripeness == 0.75 ~ '18',
                              orchard == 'Calebasses' & damage.type == 'Bat' & fruit.ripeness == 1 ~ '1',
                              orchard == 'Calebasses' & damage.type == 'Parakeet' & fruit.ripeness == 0 ~ '27',
                              orchard == 'Calebasses' & damage.type == 'Parakeet' & fruit.ripeness == 0.25 ~ '56', 
                              orchard == 'Calebasses' & damage.type == 'Parakeet' & fruit.ripeness == 0.5 ~ '12', 
                              orchard == 'Calebasses' & damage.type == 'Parakeet' & fruit.ripeness == 0.75 ~ '1',
                              orchard == 'Calebasses' & damage.type == 'Parakeet' & fruit.ripeness == 1 ~ ' ',
                              orchard == 'Calebasses' & damage.type == 'Other birds' & fruit.ripeness == 0 ~ '2',
                              orchard == 'Calebasses' & damage.type == 'Other birds' & fruit.ripeness == 0.25 ~ '34',
                              orchard == 'Calebasses' & damage.type == 'Other birds' & fruit.ripeness == 0.5 ~ '13',
                              orchard == 'Calebasses' & damage.type == 'Other birds' & fruit.ripeness == 0.75 ~ '2',
                              orchard == 'Calebasses' & damage.type == 'Other birds' & fruit.ripeness == 1 ~ ' '))
ggplot() + geom_boxplot(aes(fruit.ripeness, amount.eaten, group = fruit.ripeness, color = damage.type), data = bat_orchard_prop, size = 0.5, inherit.aes = FALSE) +
  geom_smooth(aes(x, predicted, ymin = conf.low, ymax = conf.high, fill =damage.type, shape = damage.type), data = pred4, se = FALSE, color = "black") +
  facet_wrap(~orchard+damage.type) +
  scale_color_manual(values = c("grey55", "grey55", "grey55")) + scale_y_continuous(labels = c(0, 25, 50, 75, 100)) +
  scale_x_continuous(breaks=seq(0, 1, 0.25), labels = c(0, 25, 50, 75, 100)) +
  xlab("fruit ripeness (%)") + ylab("proportion of flesh eaten (%)") + 
  theme_linedraw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=11),
        legend.text=element_text(size= 11),
        legend.title = element_blank()) +
theme(panel.background = element_rect(fill = "white"),
      legend.position = "none", 
      strip.background = element_rect(colour="white",fill="white"),
      strip.text = element_text(colour = 'black', size = 11, vjust = 2),
      panel.spacing = unit(1, "lines"),
      strip.text.x = element_text(margin = margin(2, 0, 2, 0))) +
  geom_text(data = bat_orchard_prop, aes(y = yloc, x = fruit.ripeness, label = obslabel), size = 3, check_overlap = TRUE)


#visualize trend in proportion of fruit eaten over the cycles
pred5 <- ggpredict(GLMMprop3, c("cycle", "damage.type", "orchard"), type = 're')
pred5$damage.type <- pred5$group #this line is necessary to get a nice legend
pred5$orchard <- pred5$facet #this line is necessary to get a nice legend

#test plot
ggplot() + geom_boxplot(aes(cycle, amount.eaten, group = cycle), data = bat_orchard_prop, size = 0.5, inherit.aes = FALSE) +
  geom_smooth(aes(x, predicted, ymin = conf.low, ymax = conf.high, color =damage.type, shape = damage.type), data = pred5, se = FALSE) +
  facet_wrap(~orchard+damage.type)

#visualize fruit ripeness over the cycles 
bat_orchard_prop <- bat_orchard_prop %>% 
  mutate(obslabel2 = case_when(orchard == 'Beaux Songes' & cycle == 1 ~ '77',
                               orchard == 'Beaux Songes' & cycle == 2 ~ '30',
                               orchard == 'Beaux Songes' & cycle == 3 ~ '61',
                               orchard == 'Beaux Songes' & cycle == 4 ~ '107',
                               orchard == 'Beaux Songes' & cycle == 5 ~ '87',
                               orchard == 'Beaux Songes' & cycle == 6 ~ '95',
                               orchard == 'Beaux Songes' & cycle == 7 ~ '103',
                               orchard == 'Beaux Songes' & cycle == 8 ~ '68',
                               orchard == 'Beaux Songes' & cycle == 9 ~ '75',
                               orchard == 'Beaux Songes' & cycle == 10 ~ '36',
                               orchard == 'Beaux Songes' & cycle == 11 ~ '33',
                               orchard == 'Beaux Songes' & cycle == 12 ~ '36',
                               orchard == 'Beaux Songes' & cycle == 13 ~ '36',
                               orchard == 'Calebasses' & cycle == 1 ~ '11',
                               orchard == 'Calebasses' & cycle == 2 ~ '22',
                               orchard == 'Calebasses' & cycle == 3 ~ '71',
                               orchard == 'Calebasses' & cycle == 4 ~ '87',
                               orchard == 'Calebasses' & cycle == 5 ~ '93',
                               orchard == 'Calebasses' & cycle == 6 ~ '90',
                               orchard == 'Calebasses' & cycle == 7 ~ '83',
                               orchard == 'Calebasses' & cycle == 8 ~ '93',
                               orchard == 'Calebasses' & cycle == 9 ~ '59',
                               orchard == 'Calebasses' & cycle == 10 ~ '52',
                               orchard == 'Calebasses' & cycle == 11 ~ '58'))
                              
ggplot() + geom_boxplot(aes(cycle, fruit.ripeness, group = cycle), data = bat_orchard_prop, inherit.aes = FALSE) +
  facet_wrap(~orchard) + geom_text(data = bat_orchard_prop, aes(y = yloc, x = cycle, label = obslabel2), size = 3, check_overlap = TRUE)

##################################################################
##################################################################
### now also plot ripe and unripe over the course of the study ###
##################################################################
##################################################################
adjust2 <- bat_orchard_prop %>% group_by(orchard, cycle, fruit.ripeness, factorripe) %>% summarise(sum=sum(count))
adjust2$fruit.ripeness <- adjust2$fruit.ripeness*100
adjust2$factorripe <- as.factor(adjust2$fruit.ripeness)
#first lines
ggplot(adjust2, aes(x=cycle, y=sum, group=factorripe, color=factorripe)) +
  geom_line() + facet_wrap(~orchard) #doesnt look great
#now stacked bar chart with percentages
#relevel factorripe for nicer plot
adjust2 <- adjust2 %>% mutate(factorripe = factor(factorripe, levels = c("100","75","50","25","0"))) %>% arrange(factorripe)
ggplot(adjust2, aes(fill=factorripe, y=sum, x=cycle)) + 
  geom_bar(position="fill", stat="identity") + facet_wrap(~orchard) +
  labs(fill = "fruit ripeness (%)", y = "Percentage of eaten fruits", x = "cycle") +
  scale_fill_manual(values = viridis(5, direction = -1)) +
  scale_y_continuous(labels = c(0, 25, 50, 75, 100)) +
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=11),
        legend.text=element_text(size= 10)) +
  theme(panel.background = element_rect(fill = "white"), 
        strip.background = element_rect(colour="white",fill="white"),
        strip.text = element_text(colour = 'black', size = 11))