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


# Script supporting analyses of parental behavior in the context of continuous noise 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(DHARMa)

library(ggeffects)
library(ggplot2)
library(MuMIn)

## set option for model ranking via AICc
options(na.action = "na.fail")

# load provisioning data
wdat <- read.csv("Ozkan_ContinuousNoiseData.csv", header = T)


#######################
# variable descriptions
#######################
# year = study year in which nest was monitored
# nestbox = unique identifying ID for nest box
# site = ID for study site location
# noise = time-averaged sound level at the nest box in unweighted decibels
# feedrate = rate of nestling provisioning events per hour
# enterrate = rate of adult nest box entrances per hour including successful provisioning and events when adults exited nest box without feeding nestlings
# missedfeed = rate of failed provisioning attempts per hour 
# latency = latency for parents to resume provisioning after nest disturbance
# chicks = brood size of the nest at time of experiment
# laydate = ordinal date when the first egg of the clutch was laid
# time = time of day of provisioning monitoring, expressed in decimal hours


######################################################################
######################################################################
## provisioning rate model
######################################################################
######################################################################

## first consider model with site as a random intercept
mrate <- lmer(feedrate ~ noise + chicks + laydate + time + (1|site), data = wdat, REML =F)
summary(mrate)


### singular, removed random effect because no variance explained by site as a random intercept
mrate <- lm(feedrate ~ noise + chicks + laydate + time, data = wdat)


# model selection
dr_rate <- dredge(mrate)
dr_rate

# model with noise is single competitive model
mrateTop <- lm(feedrate ~ noise, data = wdat)

# model diagnostic plots
plot(mrateTop)

summary(mrateTop)
confint(mrateTop)
AICc(mrateTop)

### box TP05B5 had high leverage, test again without 
rate_df <- subset(wdat, wdat$nestbox != "TP05B5")

mrate2 <- lm(feedrate ~ noise, data = rate_df)

summary(mrate2)  # removal of high leverage point does not change interpretation
plot(mrate2)


##### plotting using full model
df_m2 <- ggpredict(mrateTop, terms = c("noise"))

PratePlot <- plot(df_m2, add.data =T, colors = "social") + xlab("Amplitude (dB)") + 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 = "A")+ theme(plot.subtitle = element_text(face="bold", size = 14))+
  coord_cartesian(ylim = c(0,19))
P2 <- PratePlot + geom_point(data = wdat, aes(x = noise, y =feedrate), color = "tomato3")

P2





######################################################################
######################################################################
## latency to resume provisioning
######################################################################
######################################################################

## first consider model with site as a random intercept
mlag <- lmer(latency ~ noise + chicks + laydate + time + (1|site), data = wdat, REML =F)
summary(mlag)


### singular, removed random effect because no variance explained by site as a random intercept

mLag <- lm(latency ~ noise + chicks + laydate + time, data = wdat)

drLag <- dredge(mLag)
drLag


mLagTop <- lm(latency ~ noise, data = wdat)
plot(mLagTop)
summary(mLagTop)
confint(mLagTop)
AICc(mLagTop)


### box CP03B10 identified as having high leverage, removed and run model again
Lag_df <- subset(wdat, wdat$nestbox != "CP03B10")

mLag2 <- lm(latency ~ noise, data = Lag_df)
summary(mLag2) ### qualitative interpretation of results is unchanged



#### plot original
df_mLag <- ggpredict(mLagTop, terms = c("noise"))

LagPlot <- plot(df_mLag, add.data =T, colors = "social") + xlab("Amplitude (dB)") + ylab("Provisioning latency (s)") + 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 = "B")+ theme(plot.subtitle = element_text(face="bold", size = 14))
LagPlot

P3 <- LagPlot + geom_point(data = wdat, aes(x = noise, y =latency), color = "tomato3")

P3

######################################################################
######################################################################
## failed provisioning rate
######################################################################
######################################################################

m4 <- lmer(missedfeed ~ noise + chicks + laydate + time + (1|site), data = wdat, REML=F)
summary(m4)


### singular, removed random effect because no variance explained by site as a random intercept

mFrate <- lm(missedfeed ~ noise + chicks + laydate + time, data = wdat)



dr_Frate <- dredge(mFrate)
dr_Frate

## highest supported model was the intercept only model


library("gridExtra")
grid.arrange(P2, P3,  
             ncol = 2, nrow = 1)
g <- arrangeGrob(P2, P3,  
                 ncol = 2, nrow = 1)

ggsave("Figure3.pdf",g,width = 20, height = 9, units = "cm")