## Code associated with analyses of species-specific responses to noise and light and traits in Senzaki et al. (2020) Sensory pollutants alter bird phenology and fitness across a continent. Nature. ## See original paper, plus the ReadMe file for more information ## necessary packages install.packages(c("nlme","ape","geiger", "MuMIn")) library(nlme) # for gls function library(ape) library(geiger) library(MuMIn) # for AICc values dat <- read.csv("Senzaki_traits&responses.csv", header=T) row.names(dat) <- dat$Genus_species dat$cavity <- ifelse(dat$cavity == 1, "Yes", "No") dat$cavity <- as.factor(dat$cavity) tree_out<- read.tree("Jetz_ConsensusPhy.tre") #### link to tree et <- treedata(tree_out, dat, sort=T) ###################################################################### ###################################################################### ### Changes in clutch initiation date (FLD) in response to light ###################################################################### ###################################################################### # FLD and eye FLD_phy_err <- gls(FLD_Light_Est ~ corneal.transverse, data =dat, correlation = corPagel(0.5, phy=et$phy, fixed=F), weights = varFixed(~1/sqrt(FLD_Light_SE)), method = "ML") summary(FLD_phy_err) confint(FLD_phy_err) AICc(FLD_phy_err) ### ## obtain residuals following protocol here: http://nunn.rc.fas.harvard.edu/groups/pica/wiki/d8009/751_Running_PGLS_in_R_using_caper.html ### create studentized residuals and remove those that are greater than absolute value of 3 ### this code can be altered and applied to all subsequent models ind_res <- residuals(FLD_phy_err, phylo = TRUE) ind_res <- ind_res/sqrt(var(ind_res))[1] ind_res ## max(abs(ind_res)) # to see if any are over absolute value of three # FLD and diet ## lambda estimated to be negative, fixed at zero FLD_phy_errD <- gls(FLD_Light_Est ~ Diet , data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(FLD_Light_SE)), method = "ML") summary(FLD_phy_errD) confint(FLD_phy_errD, level =0.95) AICc(FLD_phy_errD) ### cavity and FLD FLD_phy_errC <- gls(FLD_Light_Est ~ cavity, data =dat, correlation = corPagel(0.1, phy=et$phy, fixed=F), weights = varFixed(~1/sqrt(FLD_Light_SE)), method = "ML") summary(FLD_phy_errC) confint(FLD_phy_errC, level=0.95) AICc(FLD_phy_errC) ###################################################################### ###################################################################### ## FLD and noise ###################################################################### ###################################################################### ## estimated lambda slightly higher than 1 so fixed to 1 FLD_phy_N_F <- gls(FLD_L50_Est ~ Freq, data =dat, correlation = corPagel(1, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(FLD_L50_SE)), method = "ML") summary(FLD_phy_N_F) confint(FLD_phy_N_F, level=0.95) AICc(FLD_phy_N_F) # diet ## estimated lambda slightly higher than 1 so fixed to 1 FLD_phy_N_D <- gls(FLD_L50_Est ~ Diet, data =dat, correlation = corPagel(1, phy=et$phy, fixed =T), weights = varFixed(~1/sqrt(FLD_L50_SE)), method = "ML") summary(FLD_phy_N_D) confint(FLD_phy_N_D, level=0.95) AICc(FLD_phy_N_D) # cavity ## estimated lambda slightly higher than 1 so fixed to 1 FLD_phy_N_C <- gls(FLD_L50_Est ~ cavity, data =dat, correlation = corPagel(1, phy=et$phy, fixed =1), weights = varFixed(~1/sqrt(FLD_L50_SE)), method = "ML") summary(FLD_phy_N_C) confint(FLD_phy_N_C, level=0.95) AICc(FLD_phy_N_C) ###################################################### ###################################################### #### clutch size and light ###################################################### ###################################################### # lambda estimated lower that zero so fixed at zero, started at 0.9 CS_phy_L_E <- gls(CS_Light_Est ~ corneal.transverse, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(CS_Light_SE)), method = "ML") summary(CS_phy_L_E) confint(CS_phy_L_E, level=0.95) AICc(CS_phy_L_E) # Diet # lambda estimated lower that zero so fixed at zero CS_phy_L_D <- gls(CS_Light_Est ~ Diet, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(CS_Light_SE)), method = "ML") summary(CS_phy_L_D) confint(CS_phy_L_D, level=0.95) AICc(CS_phy_L_D) #cavity # lambda estimated lower that zero so fixed at zero CS_phy_L_C <- gls(CS_Light_Est ~ cavity, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(CS_Light_SE)), method = "ML") summary(CS_phy_L_C) confint(CS_phy_L_C, level=0.95) AICc(CS_phy_L_C) ###################################################################### ###################################################################### ## Clutch size and noise ###################################################### ###################################################### ### fixed at 0 due to lowest AIC in range # One outlier detected (Baeolophus_inornatus). Removal did not change inference. CS_phy_N_F <- gls(CS_L50_Est ~ Freq, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(CS_L50_SE)), method = "ML") summary(CS_phy_N_F) confint(CS_phy_N_F, level=0.95) AICc(CS_phy_N_F) # diet ### lambda fixed at 0 due to estimate lower than zero # One outlier detected (Baeolophus_inornatus). Removal did not change inference. CS_phy_N_D <- gls(CS_L50_Est ~ Diet, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(CS_L50_SE)), method = "ML") summary(CS_phy_N_D) confint(CS_phy_N_D, level=0.95) AICc(CS_phy_N_D) # cavity ### lambda fixed at 0 due to estimate lower than zero # One outlier detected (Baeolophus_inornatus). Removal did not change inference. CS_phy_N_C <- gls(CS_L50_Est ~ cavity, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(CS_L50_SE)), method = "ML") summary(CS_phy_N_C) confint(CS_phy_N_C, level=0.95) AICc(CS_phy_N_C) ###################################################################### ###################################################################### #### partial hatch and light ###################################################### ###################################################### ### light gathering ## lambda < 0, fixed =0, start value = 0.7 PH_phy_L <- gls(PH_Light_Est ~ corneal.transverse, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(PH_Light_SE)), method = "ML") summary(PH_phy_L) confint(PH_phy_L, level=0.95) AICc(PH_phy_L) # Diet ## lambda < 0, fixed =0 PH_phy_L_D <- gls(PH_Light_Est ~ Diet, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(PH_Light_SE)), method = "ML") summary(PH_phy_L_D) confint(PH_phy_L_D, level=0.95) AICc(PH_phy_L_D) # cavity ## lambda < 0, fixed =0 PH_phy_L_C <- gls(PH_Light_Est ~ cavity, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), weights = varFixed(~1/sqrt(PH_Light_SE)), method = "ML") summary(PH_phy_L_C) confint(PH_phy_L_C, level=0.95) AICc(PH_phy_L_C) ###################################################################### ###################################################################### ## PH and noise ###################################################### ###################################################### ### lambda < 0, fixed = 0 PH_phy_N_F <- gls(PH_L50_Est ~ Freq, data =dat, correlation = corPagel(0, phy=et$phy,fixed =T), weights = varFixed(~1/sqrt(PH_L50_SE)), method = "ML") summary(PH_phy_N_F) confint(PH_phy_N_F, level=0.95) AICc(PH_phy_N_F) ## Diet # lambda < 0, fixed at 0 PH_phy_N_D <- gls(PH_L50_Est ~ Diet, data =dat, correlation = corPagel(0, phy=et$phy, fixed =T), weights = varFixed(~1/sqrt(PH_L50_SE)), method = "ML") summary(PH_phy_N_D) confint(PH_phy_N_D, level = 0.95) AICc(PH_phy_N_D) ## Cavity # lambda < 0, fixed at 0 PH_phy_N_C <- gls(PH_L50_Est ~ cavity, data =dat, correlation = corPagel(0, phy=et$phy, fixed =T), weights = varFixed(~1/sqrt(PH_L50_SE)), method = "ML") summary(PH_phy_N_C) confint(PH_phy_N_C, level = 0.95) AICc(PH_phy_N_C) ###################################################################### ###################################################################### ## light and clutch failure ###################################################### ###################################################### # light gathering CF_phy_L_E <- gls(CF_Light_Est ~ corneal.transverse, data =dat, correlation = corPagel(.5, phy=et$phy,fixed =F), weights = varFixed(~1/sqrt(CF_Light_SE)), method = "ML") summary(CF_phy_L_E) confint(CF_phy_L_E, level = 0.95) AICc(CF_phy_L_E) ## Diet CF_phy_L_D <- gls(CF_Light_Est ~ Diet, data =dat, correlation = corPagel(0.1, phy=et$phy,fixed =F), weights = varFixed(~1/sqrt(CF_Light_SE)), method = "ML") summary(CF_phy_L_D) confint(CF_phy_L_D, level = 0.95) AICc(CF_phy_L_D) ## Cavity CF_phy_L_C <- gls(CF_Light_Est ~ cavity, data =dat, correlation = corPagel(0.2, phy=et$phy,fixed =F), weights = varFixed(~1/sqrt(CF_Light_SE)), method = "ML") summary(CF_phy_L_C) confint(CF_phy_L_C, level = 0.95) AICc(CF_phy_L_C) ###################################################################### ###################################################################### ## CF and noise ###################################################### ###################################################### # Freq # lambda < 0, fixed at zero CF_phy_N_F <- gls(CF_L50_Est ~ Freq, data =dat, correlation = corPagel(0, phy=et$phy, fixed =T), weights = varFixed(~1/sqrt(CF_L50_SE)), method = "ML") summary(CF_phy_N_F) confint(CF_phy_N_F, level = 0.95) AICc(CF_phy_N_F) # Diet # lambda < 0, fixed at zero CF_phy_N_D <- gls(CF_L50_Est ~ Diet, data =dat, correlation = corPagel(0, phy=et$phy, fixed =T), weights = varFixed(~1/sqrt(CF_L50_SE)), method = "ML") summary(CF_phy_N_D) confint(CF_phy_N_D, level = 0.95) AICc(CF_phy_N_D) # Cavity # lambda < 0, fixed at zero CF_phy_N_C <- gls(CF_L50_Est ~ cavity, data =dat, correlation = corPagel(0, phy=et$phy, fixed =T), weights = varFixed(~1/sqrt(CF_L50_SE)), method = "ML") summary(CF_phy_N_C) confint(CF_phy_N_C, level = 0.95) AICc(CF_phy_N_C) ###################################################################### ###################################################################### ## light and success ###################################################### ###################################################### ## remove house sparrow d2 <- subset(dat, dat$species != "houspa") etS <- treedata(tree_out, d2, sort=T) # light gathering ability of eye # lambda < 0, fixed at 0, start value = 0.25 S_phy_L_E <- gls(Succ_Light_Est ~ corneal.transverse, data =d2, correlation = corPagel(0, phy=etS$phy,fixed=T), weights = varFixed(~1/sqrt(Succ_Light_SE)), method = "ML") summary(S_phy_L_E) confint(S_phy_L_E, level = 0.95) AICc(S_phy_L_E) ## diet # lambda < 0, fixed at 0 S_phy_L_D <- gls(Succ_Light_Est ~ Diet, data =d2, correlation = corPagel(0, phy=etS$phy,fixed=T), weights = varFixed(~1/sqrt(Succ_Light_SE)), method = "ML") summary(S_phy_L_D) confint(S_phy_L_D, level = 0.95) AICc(S_phy_L_D) ## Cavity # lambda < 0, fixed at 0 S_phy_L_C <- gls(Succ_Light_Est ~ cavity, data =d2, correlation = corPagel(0, phy=etS$phy,fixed=T), weights = varFixed(~1/sqrt(Succ_Light_SE)), method = "ML") summary(S_phy_L_C) confint(S_phy_L_C, level = 0.95) AICc(S_phy_L_C) ###################################################################### ###################################################################### ## Succ and noise ###################################################### ###################################################### ### Freq # lambda < 0, fixed at zero # One outlier detected (Baeolophus_inornatus). Removal did not change inference. S_phy_N <- gls(Succ_L50_Est ~ Freq, data =d2, correlation = corPagel(0, phy=etS$phy, fixed =T), weights = varFixed(~1/sqrt(Succ_L50_SE)), method = "ML") summary(S_phy_N) confint(S_phy_N, level = 0.95) AICc(S_phy_N) # Diet # lambda < 0, fixed at zero # One outlier detected (Baeolophus_inornatus). Removal did not change inference. S_phy_N_D <- gls(Succ_L50_Est ~ Diet, data =d2, correlation = corPagel(0, phy=etS$phy, fixed =T), weights = varFixed(~1/sqrt(Succ_L50_SE)), method = "ML") summary(S_phy_N_D) confint(S_phy_N_D, level = 0.95) AICc(S_phy_N_D) # cavity # lambda < 0, fixed at zero # One outlier detected (Baeolophus_inornatus). Removal did not change inference. S_phy_N_C <- gls(Succ_L50_Est ~ cavity, data =d2, correlation = corPagel(0, phy=etS$phy, fixed = T), weights = varFixed(~1/sqrt(Succ_L50_SE)), method = "ML") summary(S_phy_N_C) confint(S_phy_N_C, level = 0.95) AICc(S_phy_N_C) ####################################### ####################################### ## sensitivity analysis restricted to only species without species where corneal.transverse was imputed ####################################### ####################################### # first, success without house sparrow d2NoImp <- subset(d2, d2$Impute != "Yes") etS2 <- treedata(tree_out, d2NoImp, sort=T) ## success S_phy_L_E <- gls(Succ_Light_Est ~ corneal.transverse, data =d2NoImp, correlation = corPagel(0, phy=etS2$phy,fixed=T), weights = varFixed(~1/sqrt(Succ_Light_SE)), method = "ML") summary(S_phy_L_E) confint(S_phy_L_E, level = 0.95) AICc(S_phy_L_E) ## from whole dataset dat2 <- subset(dat, dat$Impute != "Yes") et2 <- treedata(tree_out, dat2, sort=T) CF_phy_L_E <- gls(CF_Light_Est ~ corneal.transverse, data =dat2, correlation = corPagel(0, phy=et2$phy,fixed =T), weights = varFixed(~1/sqrt(CF_Light_SE)), method = "ML") summary(CF_phy_L_E) confint(CF_phy_L_E, level = 0.95) AICc(CF_phy_L_E) ## lambda > 1, fixed =1 PH_phy_L <- gls(PH_Light_Est ~ corneal.transverse, data =dat2, correlation = corPagel(1, phy=et2$phy, fixed=T), weights = varFixed(~1/sqrt(PH_Light_SE)), method = "ML") summary(PH_phy_L) confint(PH_phy_L, level=0.95) AICc(PH_phy_L) # lambda estimated lower that zero so fixed at zero CS_phy_L_E <- gls(CS_Light_Est ~ corneal.transverse, data =dat2, correlation = corPagel(1, phy=et2$phy, fixed=F), weights = varFixed(~1/sqrt(CS_Light_SE)), method = "ML") summary(CS_phy_L_E) confint(CS_phy_L_E, level=0.95) # FLD and eye FLD_phy_err <- gls(FLD_Light_Est ~ corneal.transverse, data =dat2, correlation = corPagel(1, phy=et2$phy, fixed=T), weights = varFixed(~1/sqrt(FLD_Light_SE)), method = "ML") summary(FLD_phy_err) confint(FLD_phy_err) AICc(FLD_phy_err) ##### corneal transverse model # converged < 0 so fixed at zero dat$H2 <- ifelse(dat$Habitat == "o", "aO", as.character(dat$Habitat)) CTmod <- gls(corneal.transverse ~ H2, data =dat, correlation = corPagel(0, phy=et$phy, fixed=T), method = "ML") summary(CTmod) AICc(CTmod) confint(CTmod)