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


# Script supporting main analyses of short-term noise playback experiment within “Divergent effects of short-term and continuous anthropogenic noise exposure on Western Bluebird parental care behavior”

# Contact Authors:

#Name: Kerstin Ozkan
#Institution: Department of Biological Sciences, California Polytechnic State University
#Address: 1 Grand Ave, San Luis Obispo, California, USA 93407
#Email: kerstin.ozkan@gmail.com

#Name: Clinton D. Francis
#Institution: Department of Biological Sciences, California Polytechnic State University
#Address: 1 Grand Ave, San Luis Obispo, California, USA 93407
#Email: cdfranci@calpoly.edu


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

## load necessary packages
library(lme4)
library(lmerTest)
library(ggplot2)
library(lubridate)
library(MuMIn)
library(dplyr)
options(na.action = "na.fail")

library(blmeco) ## model diagnostics
library(performance) ## model diagnostics
library(DHARMa) ## model diagnostics
library(ggeffects)

####################################################################################################
####################################################################################################
# set working directory

## load data

dat <- read.csv("Ozkan_ShortNoiseData.csv")

################################################
# variable descriptions
################################################
# nestid = nest box number
# parent = identifies record to activity from the male or female parent
# trial = indicates whether the record was the first or second trial
# treatment = indicates whether the trial was a noise playback or silent control
# visits = rate of nestling provisioning events per hour
# date = date the trial was conducted
# number_chicks = brood size of the nest at time of experiment
# track = number reflects 1 of six traffic noise playback tracks that were randomly selected for use at a nest
# tw10m = time (in seconds) parent spent perched within 10m of the nest box
# failed_attempt = rate of failed provisioning attempts per hour





dat$trial<- as.factor(as.character(dat$trial)) #changing integer to categorical for trial

## add ordinal date

as.Date(dat$date) #, format = "%Y-%m-%d")

dat$ordinal<- yday(dat$date)

### how are visits distributed?

hist(dat$visits) # right skew

# set optimizer for glmer
control <- glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5))  

################################################################
################################################################
# provisioning rate models
################################################################
################################################################


m1 <- glmer(visits ~ treatment*trial + number_chicks + scale(ordinal) + scale(tw10m) + parent + (1|nestid), data = dat, family = poisson, control = control) 

## model diagnostics
## check dispersion, just for poisson
dispersion_glmer(m1) #value should be between 0.75 and 1.4 are ok. looks good at 1.094

### check residuals with qq plots and residual vs. predicted
simulationOutput <- simulateResiduals(fittedModel = m1, plot = F)

plot(simulationOutput) 

# model selection to determine subset of variables that best explain the data
dr <- dredge(m1)

# see top five models
head(dr) 

# top model

mTtop1 <- glmer(visits ~  treatment*trial + parent + scale(tw10m) + (1|nestid), data = dat, family = poisson, control = control)

## check dispersion. This is good (1.09)
dispersion_glmer(mTtop1)


### check residuals
simulationOutput1 <- simulateResiduals(fittedModel = mTtop1, plot = F)

plot(simulationOutput1) ## this look good. 

summary(mTtop1)

confint(mTtop1)

################################################################
################################################################
# time spent with 10m 
################################################################
################################################################
m2 <- lmer(log(tw10m+1) ~ treatment*trial + number_chicks + scale(ordinal)  + parent + (1|nestid), data = dat, REML = F) 


### check residuals with qq plots and residual vs. predicted
simulationOutput2 <- simulateResiduals(fittedModel = m2, plot = F)

plot(simulationOutput2) 
check_model(m2)

# model selection to determine subset of variables that best explain the data
dr2 <- dredge(m2)

head(dr2,10) 

### second model
m2v2 <- lmer(log(tw10m+1) ~ treatment*trial  + (1|nestid), data = dat, REML = F) 
summary(m2v2)
confint(m2v2)

dfTime <- ggpredict(m2v2, terms= c("trial", "treatment"))
plot(dfTime)

################################################################
################################################################
# failed provisioning attempts 
################################################################
################################################################


mFA <- glmer.nb(failed_attempt ~ treatment*trial + number_chicks + scale(ordinal) + scale(tw10m) + parent + (1|nestid), data = dat, control = control) #

## model diagnostics
# now look at qqplots. We want to see that the points generally fall on the line
### check residuals
simulationOutputFA <- simulateResiduals(fittedModel = mFA, plot = F)

plot(simulationOutputFA) ## these look ok

check_zeroinflation(mFA) ### this model seems ok


# model selection to determine subset of variables that best explain the data
drFA <- dredge(mFA)

head(drFA) 


candSet <- model.avg(get.models(drFA,subset = delta <=2))
summary(candSet)

mFADF <- ggpredict(mFA, terms=c("parent","treatment","trial"))
plot(mFADF)


### model selection

drFA <- dredge(mFA)

head(drFA,35)


## top model

mFAtop <- glmer.nb(failed_attempt ~ treatment + parent + number_chicks+ (1|nestid), data = dat, control = control)
summary(mFAtop)

simulationOutputFAtop <- simulateResiduals(fittedModel = mFAtop, plot = F)

plot(simulationOutputFAtop) ## these look ok

check_zeroinflation(mFAtop) ### this model seems ok

summary(mFAtop) ##
confint(mFAtop)


mFADFtop <- ggpredict(mFAtop, terms=c("treatment", "parent"))
plot(mFADFtop, add.data =T)


# null
mFAnull <- glmer.nb(failed_attempt ~ 1 + (1|nestid), data = dat, control = control)
AICc(mFAnull)


#### sensitivity analysis removing female from SR065
dat$nestidParent <- paste(dat$nestid,dat$parent)
filterdat <- dat %>% filter(nestidParent != "SR065 female")

# Top model from above with female removed
control_bobyqa <- glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))

mFAsubTop <- glmer.nb(failed_attempt ~ treatment + number_chicks +  parent + (1|nestid), data = filterdat, control = control_bobyqa)

### model diagnostics
simulationOutputFAsubTop <- simulateResiduals(fittedModel = mFAsubTop, plot = F)

plot(simulationOutputFAsubTop) ## these look ok

check_zeroinflation(mFAsubTop) ### this model seems ok


summary(mFAsubTop)

# confint function failing to converge. For this model, 95% CI were estimated using the estimate +/- (SE*1.96)
confint(mFAsubTop)

dfFHsub <- ggpredict(mFAsubTop, "treatment")

#######################
#######################
# Figure 2 plot
#######################
#######################

dat$treattrial <- paste(dat$treatment,dat$trial)
dat$nestidPar <- paste(dat$nestid,dat$parent)
dat$treatOrder <- ifelse(dat$treattrial == "Control 1" | dat$treattrial == "Noise 2", "Control First", "Treatment First")
dat$treatOrder <- as.factor(dat$treatOrder)
#######################


### first panel control first and noise second
### reformat data
datC1 <- dat %>%
  filter(treatOrder == "Control First")


b <- ggplot(datC1, aes(x = trial, y = visits, color =parent, linetype = parent, shape = parent)) + 
  geom_point(position = position_dodge(width =0.025), size=1.5) +
  geom_line(aes(group=nestidPar), lwd=0.5, alpha=0.25) +
  stat_summary(aes(group = parent), fun.y = "mean", geom = "line", size = 2) +
  stat_summary(aes(group = parent), fun.data = "mean_se", 
               geom = "pointrange", size = 1, position = position_dodge(.1)) + 
  scale_color_manual(values = c("tomato3","steelblue4")) +
  scale_linetype_manual(values = c("solid","dotted")) +
  scale_shape_manual(values = c(16, 15)) +
  theme_classic()  +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1),                        legend.position = "none", axis.line = element_line(colour = "black"),
        axis.text.x = element_text(color = "black", size = 12),axis.text.y =                         element_text(color = "black", size = 12),axis.title.x = element_text                         (color = "black", size = 12),
        axis.title.y = element_text(color = "black", size = 12))+
  labs(subtitle = "A")+ theme(plot.subtitle = element_text(face="bold",                        size = 14)) + xlab("")+
  ylab(expression("Visits" ~h^-1))+scale_x_discrete(labels=c("Ambient", "Noise"))+
  coord_cartesian(ylim = c(0,13))+ theme(plot.margin = unit(c(0.1, 0.2, 0, 0.1), "cm"))
b


### Now noise first and control second
datT1 <- dat %>%
  filter(treatOrder == "Treatment First")


c <- ggplot(datT1, aes(x = trial, y = visits, color =parent, linetype = parent, shape =parent)) + 
  geom_point(position = position_dodge(width =0.025), size=1.5) +
  geom_line(aes(group=nestidPar), lwd=0.5, alpha=0.25) +
  stat_summary(aes(group = parent), fun.y = "mean", geom = "line", size = 2) +
  stat_summary(aes(group = parent), fun.data = "mean_se", 
               geom = "pointrange", size = 1, position = position_dodge(.1)) + 
  scale_color_manual(values = c("tomato3","steelblue4")) +
  scale_linetype_manual(values = c("solid","dotted")) +
  scale_shape_manual(values = c(16, 15)) +
  theme_classic()  +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1),                        legend.position = "none", axis.line = element_line(colour = "black"),
        axis.text.x = element_text(color = "black", size = 12),axis.text.y =                         element_text(color = "black", size = 12),axis.title.x = element_text                         (color = "black", size = 12),
        axis.title.y = element_text(color = "black", size = 12))+
  labs(subtitle = "B")+ theme(plot.subtitle = element_text(face="bold",                        size = 14)) + xlab("")+
  ylab(expression("Visits" ~h^-1))+scale_x_discrete(labels=c("Noise", "Ambient")) +
  coord_cartesian(ylim = c(0,13))+ theme(plot.margin = unit(c(0, 0.2, 0, 0.1), "cm"))
c



#### plot time spent within 10 marginal effect on provisioning
dftime <- ggpredict(mTtop1, terms = "tw10m")

Pt10_1 <- plot(dftime, show_data=T, colors = "simply")

Pt10_2 <- Pt10_1 + geom_point(data = dat, aes(x = tw10m, y =visits), color = "tomato3") + xlab("Time spent within 10m (s)") + ylab(expression("Visits" ~h^-1)) + theme_classic()  +theme(panel.border = element_rect(colour = "black", fill=NA, size=1), panel.grid.major = element_blank(),  panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),  axis.text.x = element_text(color = "black", size = 12),axis.text.y = element_text(color = "black", size = 12),axis.title.x = element_text(color = "black", size = 12),axis.title.y = element_text(color = "black", size = 12)) + labs(title ="") +
  labs(subtitle = "C")+ theme(plot.subtitle = element_text(face="bold", size = 14))+ theme(plot.margin = unit(c(0, 0.2, 0.2, 0.1), "cm"))

Pt10_2




library("gridExtra")
grid.arrange(b, c, Pt10_2, 
             ncol = 1, nrow = 3, top = NULL, bottom = NULL)
g <- arrangeGrob(b, c,  Pt10_2,
                 ncol = 1, nrow = 3, top = NULL, bottom = NULL)

ggsave("Figure2.pdf",g,width = 10, height = 30, units = "cm")