#R Script for RWBB Clutch Size statistical analyses

#Robinson SK, McChesney HM [in prep] Nesting success of Red-winged Blackbirds (Agelaius phoeniceus) in marshes in an anthropogenic landscape.
#Please refer to the manuscript for data collection methods and statistical analyses. For questions or to notify the authors if any errors are identified in the data, please contact Holly McChesney (holly.mcchesney@arcadis.com) or Scott Robinson (srobinson@flmnh.ufl.edu).

###########################################
#data input
###########################################
library("lme4")
library("MuMIn")
library("jtools")
library("ggplot2")
library("tidyverse")
library("Hmisc")
library("boot")
library(effects)
library(sjPlot)
library(sjmisc)
library(ggplot2)
library(plotly)
library(scales)

#Set to appropriate directory
setwd("~/")

# Analyses below look at clutch mass in 3 complementary datasets (A=all data, B=predation excluded, C=predation and flooding excluded)


clutch.A <- na.omit(read.csv("1clutch.data.csv")) 

#Set distances less than 1 to 1
clutch.A$dist2edge[clutch.A$dist2edge < 1] <- 1
clutch.A$dist2water[clutch.A$dist2water < 1] <- 1


###########################################
#Transformations and centering
###########################################
#treating year as categorical, for random effect (creating an object named year):
year.a<-as.factor(clutch.A$Year)

#Transforming and/or scaling all of the numeric variables:
DoYr_ctr.a<-scale(clutch.A$IncDoYr)
PCB_ctr.a<-scale(clutch.A$PCB_class)
distWtld_ctr.a<-scale(log10(clutch.A$dist2wtlnd))
density_ctr.a<-scale(sqrt(clutch.A$density))
distEdge_ctr.a<-scale(log10(clutch.A$dist2edge))
distWater_ctr.a<-scale(log10(clutch.A$dist2water))
NestID<-clutch.A$NestID
nestHt.a<-clutch.A$nestHt
Year<-clutch.A$Year
loc<-clutch.A$loc

###########################################
#Examined pairwise correlation of fixed effects.
###########################################
fixedeffects.matrix.a<-cbind(DoYr_ctr.a,PCB_ctr.a,distWater_ctr.a,density_ctr.a, distEdge_ctr.a)
pairwise1.a <- rcorr(fixedeffects.matrix.a)
pairwise1.a

###########################################
#Clutch models with no Random Effects
###########################################

#Apriori models:
null.A<-glm(clutch.A$clutch.category ~ 1, na.action="na.fail", family="binomial")

global.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+distWater_ctr.a*(nestHt.a+PCB_ctr.a)+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

DOY.A<-glm(clutch.A$clutch.category~DoYr_ctr.a, na.action="na.fail", family="binomial")

PCB.A<-glm(clutch.A$clutch.category~PCB_ctr.a, na.action="na.fail", family="binomial")

Diet.A<-glm(clutch.A$clutch.category~(distWater_ctr.a), na.action="na.fail", family="binomial")

Flood.A<-glm(clutch.A$clutch.category~nestHt.a*distWater_ctr.a, na.action="na.fail", family="binomial")

Predation.A<-glm(clutch.A$clutch.category~nestHt.a+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

DOY_PCB.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+PCB_ctr.a, na.action="na.fail", family="binomial")

DOY_Diet.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+distWater_ctr.a, na.action="na.fail", family="binomial")

DOY_Flood.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+nestHt.a*distWater_ctr.a, na.action="na.fail", family="binomial")

DOY_Predation.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

DOY_Flood_Predation.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

DOY_Diet_Predation.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+distWater_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

PCB_Diet.A<-glm(clutch.A$clutch.category~PCB_ctr.a*distWater_ctr.a, na.action="na.fail", family="binomial")

DOY_PCB_Diet.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+PCB_ctr.a*distWater_ctr.a, na.action="na.fail", family="binomial")

PCB_Flood.A<-glm(clutch.A$clutch.category~PCB_ctr.a+nestHt.a*distWater_ctr.a, na.action="na.fail", family="binomial")

DOY_PCB_Flood.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+PCB_ctr.a+nestHt.a*distWater_ctr.a, na.action="na.fail", family="binomial")

PCB_Predation.A<-glm(clutch.A$clutch.category~PCB_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

DOY_PCB_Predation.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+PCB_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

PCB_Diet_Predation.A<-glm(clutch.A$clutch.category~PCB_ctr.a*distWater_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

DOY_PCB_Diet_Predation.A<-glm(clutch.A$clutch.category~DoYr_ctr.a+PCB_ctr.a*distWater_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

Diet_Predation.A<-glm(clutch.A$clutch.category~nestHt.a+(distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

Flood_Predation.A<-glm(clutch.A$clutch.category~(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")

PCB_Flood_Predation.A<-glm(clutch.A$clutch.category~PCB_ctr.a+(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a), na.action="na.fail", family="binomial")


#Model selection:  
  
models.A<-model.sel(null.A, global.A, DOY.A, PCB.A, Diet.A, Flood.A, Predation.A, DOY_PCB.A, DOY_Diet.A, DOY_Flood.A, DOY_Predation.A, DOY_Flood_Predation.A, DOY_Diet_Predation.A, PCB_Diet.A, PCB_Flood.A, PCB_Predation.A, DOY_PCB_Diet.A, DOY_PCB_Flood.A, PCB_Diet_Predation.A, DOY_PCB_Predation.A, DOY_PCB_Diet_Predation.A, Diet_Predation.A, Flood_Predation.A, PCB_Flood_Predation.A)

View(models.A)

print(subset(models.A, delta<2), abbrev = FALSE) 
  
#Top-ranked model, used in plots and for comparison with model averaged results:
clutchtopmod.a<-get.models(models.A,1)

summ(clutchtopmod.a[[1]]) 


###########################################
#Plot DOY Figure (Figure 3a)
###########################################


#Gets Transformed predictors and CIs
dat<-get_model_data(clutchtopmod.a[[1]], type="eff", show.values = TRUE)
x<-dat$DoYr_ctr.a$x
a<-attributes(DoYr_ctr.a)
x.orig = x*a$`scaled:scale` + a$'scaled:center'
DOY_plot.data<-as.data.frame(cbind(x.orig,(dat$DoYr_ctr.a$predicted), (dat$DoYr_ctr.a$conf.low), (dat$DoYr_ctr.a$conf.high)))
colnames(DOY_plot.data)<-c("x", "y", "lower", "upper")

labs<-c("≤3", "≥4")

p <- ggplot(DOY_plot.data) + geom_line(aes(y=y, x=x))+
    geom_ribbon(aes(ymin=lower, ymax=upper, x=x), alpha = 0.3)+
    scale_colour_manual("",values="black")+
    scale_fill_manual("",values="grey12")+
    labs(x = "Day of Year", y = "Probability that Cluch Size \nis at least 4 Eggs")+
    theme_classic()+
    theme(axis.text=element_text(size=16), axis.title=element_text(size=16))


fig <- ggplotly(p)
fig


###########################################
#Plot Distance to Water Figure (Figure3b)
###########################################

x3<-dat$distWater_ctr.a$x
a3<-attributes(distWater_ctr.a)
x3.orig = x3*a3$`scaled:scale` + a3$'scaled:center'
x3.orig_bt<-10^x3.orig
dWater_plot.data<-as.data.frame(cbind(x3.orig_bt,dat$distWater_ctr.a$predicted, dat$distWater_ctr.a$conf.low, dat$distWater_ctr.a$conf.high))
colnames(dWater_plot.data)<-c("x", "y", "lower", "upper")

p2 <- ggplot(dWater_plot.data) + geom_line(aes(y=y, x=x))+
    #scale_x_log10(breaks = c(0.01,0.1,1,10,100),labels = comma) +
    geom_ribbon(aes(ymin=lower, ymax=upper, x=x), alpha = 0.3)+
    scale_colour_manual("",values="black")+
    scale_fill_manual("",values="grey12")+
    labs(x = "Distance to Water (m)", y = "Probability that Cluch Size \nis at least 4 Eggs")+
    theme_classic()+
    theme(axis.text=element_text(size=16), axis.title=element_text(size=16))
    

fig2 <- ggplotly(p2)

fig2