#R Script for RWBB Egg Mass 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)

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

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

egg.A <- na.omit(read.csv("4egg.data.csv")) 


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


###########################################
#Transformations and centering
###########################################

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


#Transforming and/or scaling all of the numeric variables:
DoYr_ctr.a<-scale(egg.A$DayofYr)
PCB_ctr.a<-scale(egg.A$PCBclass)
incDay_ctr.a<-scale(egg.A$incDay)
density_ctr.a<-scale(sqrt(egg.A$density))
distEdge_ctr.a<-scale(log10(egg.A$dist2edge))
distWater_ctr.a<-scale(log10(egg.A$dist2water))
NestID<-egg.A$NestID
nestHt.a<-egg.A$nestHt
loc<-egg.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
```


###########################################
#Egg Mass models with Nest ID as Random Variable
###########################################


#Apriori models:
null.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+(1|NestID), na.action="na.fail")

global.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+distWater_ctr.a*(nestHt.a+PCB_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

DOY.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+(1|NestID), na.action="na.fail")

PCB.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+PCB_ctr.a+(1|NestID), na.action="na.fail")

Diet.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+(distWater_ctr.a)+(1|NestID), na.action="na.fail")

Flood.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+nestHt.a*distWater_ctr.a+(1|NestID), na.action="na.fail")

Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

DOY_PCB.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+PCB_ctr.a+(1|NestID), na.action="na.fail")

DOY_Diet.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+distWater_ctr.a+(1|NestID), na.action="na.fail")

DOY_Flood.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+nestHt.a*distWater_ctr.a+(1|NestID), na.action="na.fail")

DOY_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

DOY_Flood_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

DOY_Diet_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+distWater_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

PCB_Diet.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+PCB_ctr.a*distWater_ctr.a+(1|NestID), na.action="na.fail")

DOY_PCB_Diet.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+PCB_ctr.a*distWater_ctr.a+(1|NestID), na.action="na.fail")

PCB_Flood.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+PCB_ctr.a+nestHt.a*distWater_ctr.a+(1|NestID), na.action="na.fail")

DOY_PCB_Flood.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+PCB_ctr.a+nestHt.a*distWater_ctr.a+(1|NestID), na.action="na.fail")

PCB_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+PCB_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

DOY_PCB_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+PCB_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

PCB_Diet_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+PCB_ctr.a*distWater_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

DOY_PCB_Diet_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+DoYr_ctr.a+PCB_ctr.a*distWater_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

Diet_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+nestHt.a+(distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

Flood_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")

PCB_Flood_Predation.A<-lmer(egg.A$MASS_g^2~incDay_ctr.a+PCB_ctr.a+(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID), na.action="na.fail")


#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:
eggtopmod.a<-get.models(models.A,1)

summ(eggtopmod.a[[1]]) 

modAveTop.a<-model.avg(DOY.A, DOY.A)
summary(modAveTop.a)

###########################################
#Plot DOY Figure (Figure 4a)
###########################################

#Gets Transformed predictors and CIs
dat<-get_model_data(eggtopmod.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,sqrt(dat$DoYr_ctr.a$predicted), sqrt(dat$DoYr_ctr.a$conf.low), sqrt(dat$DoYr_ctr.a$conf.high)))
colnames(DOY_plot.data)<-c("x", "y", "lower", "upper")

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 = "Egg Mass (g)")+
    theme_classic()+
    theme(axis.text=element_text(size=16), axis.title=element_text(size=16))

fig <- ggplotly(p)

fig

###########################################
#Plot Incubation Dayr Figure (Figure4b)
###########################################


dat<-get_model_data(eggtopmod.a[[1]], type="eff", show.values = TRUE)
x2<-dat$incDay_ctr.a$x
a2<-attributes(incDay_ctr.a)
x2.orig = x2*a2$`scaled:scale` + a2$'scaled:center'
incDay_plot.data<-as.data.frame(cbind(x2.orig,sqrt(dat$incDay_ctr.a$predicted), sqrt(dat$incDay_ctr.a$conf.low), sqrt(dat$incDay_ctr.a$conf.high)))
colnames(incDay_plot.data)<-c("x", "y", "lower", "upper")

p2 <- ggplot(incDay_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 = "Incubation Day", y = "Egg Mass (g)")+
    theme_classic()+
    theme(axis.text=element_text(size=16), axis.title=element_text(size=16))

fig2 <- ggplotly(p2)

fig2