#R Script for RWBB Nestling 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)
library(scales)

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

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


nestling.A <- na.omit(read.csv("5nestling.data.csv")) 


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



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

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


#Transforming and/or scaling all of the numeric variables:
DoYr_ctr.a<-scale(nestling.A$IncDofYr)
age_ctr.a<-scale(nestling.A$EstNstlAge)
PCB_ctr.a<-scale(nestling.A$PCBclass)
density_ctr.a<-scale(sqrt(nestling.A$density))
distEdge_ctr.a<-scale(log10(nestling.A$dist2edge))
distWater_ctr.a<-scale(log10(nestling.A$dist2water))
NestID<-nestling.A$NestID
nestHt.a<-nestling.A$nestHt
clutch.a<-nestling.A$CLUTCH.SIZE
loc<-nestling.A$loc


###########################################
#Examined pairwise correlation of fixed effects.
###########################################

fixedeffects.matrix.a<-cbind(age_ctr.a, DoYr_ctr.a,PCB_ctr.a,distWater_ctr.a,density_ctr.a, distEdge_ctr.a)
pairwise1.a <- rcorr(fixedeffects.matrix.a)
pairwise1.a
```




###########################################
#Nestling Mass models with nest-ID as random effect
###########################################

#use nest ID as random effect

#Apriori models:
null.A<-lmer(nestling.A$nstl.mass.g~(1|NestID) , na.action="na.fail", REML=F)

global.A<-lmer(nestling.A$nstl.mass.g~age_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", REML=F)

DOY.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

PCB.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+PCB_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

Diet.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+(distWater_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

Flood.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+nestHt.a*distWater_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

DOY_PCB.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+PCB_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

DOY_Diet.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+distWater_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

DOY_Flood.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+nestHt.a*distWater_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

DOY_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

DOY_Flood_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

DOY_Diet_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+distWater_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

PCB_Diet.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+PCB_ctr.a*distWater_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

DOY_PCB_Diet.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+PCB_ctr.a*distWater_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

PCB_Flood.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+PCB_ctr.a+nestHt.a*distWater_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

DOY_PCB_Flood.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+PCB_ctr.a+nestHt.a*distWater_ctr.a+(1|NestID) , na.action="na.fail", REML=F)

PCB_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+PCB_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

DOY_PCB_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+DoYr_ctr.a+PCB_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

PCB_Diet_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+PCB_ctr.a*distWater_ctr.a+nestHt.a+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

DOY_PCB_Diet_Predation.A<-lmer(nestling.A$nstl.mass.g~age_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", REML=F)

Diet_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+nestHt.a+(distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

Flood_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)

PCB_Flood_Predation.A<-lmer(nestling.A$nstl.mass.g~age_ctr.a+PCB_ctr.a+(nestHt.a*distWater_ctr.a)+(density_ctr.a*distEdge_ctr.a)+(1|NestID) , na.action="na.fail", REML=F)


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

summ(nestlingtopmod.a[[1]]) 

modAveTop.a<-model.avg(models.A, delta<2)
summary(modAveTop.a)

#age, PCB, density, and edges significant


###########################################
#Plot PCB Figure (Figure 7b)
###########################################
#Gets Transformed predictors and CIs
dat<-get_model_data(nestlingtopmod.a[[1]], type="eff", group.terms= "age_ctr.a", show.values = TRUE)
x<-dat$PCB_ctr.a$x
a<-attributes(PCB_ctr.a)
x.orig = x*a$`scaled:scale` + a$'scaled:center'
PCB_plot.data<-as.data.frame(cbind(x.orig,(dat$PCB_ctr.a$predicted), (dat$PCB_ctr.a$conf.low), (dat$PCB_ctr.a$conf.high)))
colnames(PCB_plot.data)<-c("x", "y", "lower", "upper")

p <- ggplot(PCB_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 = "PCB Category", y = "Nestling Mass (g)")+
    theme_classic()+
    theme(axis.text=element_text(size=16), axis.title=element_text(size=16), plot.margin = margin(10,10,20,10))+
    ylim(29,32)

fig1 <- ggplotly(p)

fig1

###########################################
#Plot Distance to Edge Figure (Figure 7c)
###########################################
x2<-dat$distEdge_ctr.a$x
a2<-attributes(distEdge_ctr.a)
x2.orig = x2*a2$`scaled:scale` + a2$'scaled:center'
x2.orig_bt<-10^x2.orig
#x.orig<-10^x.orig_log
dEdge_plot.data<-as.data.frame(cbind(x2.orig_bt,dat$distEdge_ctr.a$predicted, dat$distEdge_ctr.a$conf.low, dat$distEdge_ctr.a$conf.high))
colnames(dEdge_plot.data)<-c("x", "y", "lower", "upper")

p2 <- ggplot(dEdge_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 Edge (m)", y = "Nestling Mass (g)")+
    theme_classic()+
    theme(axis.text=element_text(size=16), axis.title=element_text(size=16), plot.margin = margin(10,10,20,10))
    

fig2 <- ggplotly(p2)

fig2

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

p3 <- ggplot(age_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 = "Nestling Age", y = "Nestling Mass (g)")+
    theme_classic()+
    theme(axis.text=element_text(size=16), axis.title=element_text(size=16), plot.margin = margin(10,10,20,10))
    

fig3 <- ggplotly(p3)

fig3

###########################################
#Plot Density Figure (Figure 7d)
###########################################
#Gets Transformed predictors and CIs
#dat<-get_model_data(nestlingtopmod.a[[1]], type="eff", terms = c("density_ctr.a", "age_ctr.a [-1.57,1.90]"))
x4<-dat$density_ctr.a$x
a4<-attributes(density_ctr.a)
x4.orig = (x4*a4$'scaled:scale' + a4$'scaled:center')^2
density_plot.data<-as.data.frame(cbind(x4.orig,(dat$density_ctr.a$predicted), (dat$density_ctr.a$conf.low), (dat$density_ctr.a$conf.high)))
colnames(density_plot.data)<-c("x", "y", "lower", "upper")

p4 <- ggplot(density_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 = "Nest Density (number per hectare)", y = "Nestling Mass (g)")+
    theme_classic()+
    theme(axis.text=element_text(size=16), axis.title=element_text(size=16), plot.margin = margin(10,10,20,10))
    

fig4 <- ggplotly(p4)

fig4