############################################################################ ########### ###### ########### Code for modelling MOCH fecundity from BC Nest Web Project ###### ########### ###### ############################################################################# ## Nest data from Andrea Norris & Kathy Martin ## Weather data from William's Lake Weather Station ## Bark Beetle data from station veg (Kathy Martin & Andrea Norris) ## Lepidopteran data from integration of station veg with BC Aerial Insect Surveys ## Code last updated 13 September 2022 by Kristina Cockle rm(list = ls()) dev.off() ####################################################### ########### Modelling approach ################# ####################################################### ### Predictor variables # Year is a random effect # Some predictor variables are measured at the nest*year level (DFE, clutch size, etc.) # Some predictor variables are measured at site*year level (site-level beetle, site-level leps, MOCH, TAHU) # Some predictor variables are measured at year level across all sites (regional beetle and weather) ### Response variables # DFE (laying date) is a continuous response variable and uses a glmm (gaussian) # Clutch size (CS) is a count variable but it is underdispersed (no zeros by default), so uses COM-Poisson GLMM # Number fledged (Brood Size) is also a count variable and also underdispersed (no zeros by default), so uses COM-Poisson GLMM # Daily Nest Survival is modelled using logistic exposure to account for exposure days ########################################## ######## Packages ####################### ######################################### library(dplyr) library(nlme) #for glmm predicting laying date library(MuMIn) #for model selection library(glmmTMB) #for COM-Poisson glmm predicting clutch size and n-fledged (underdispersed count data) library(viridis) #colour palette library(scales) library(ggplot2) #graphs library(gridExtra) #arrange ggplot figs ############################################################################ ################## Data ################################################### ############################################################################ TCL <- read.csv("allmoch13.csv", header = TRUE, fileEncoding = 'UTF-8-BOM') MOCH <- data.frame ( Year = TCL$Year, Nest.ID = TCL$NestID, #unique ID for each nesting attempt Site = TCL$Site, # To calculate incubation & nestling periods dt.inc.init = TCL$Dt_inc_init, #incubation initiation date date.hatch = TCL$Dt_hatch, #hatch date dt.hatch.cert = TCL$Dt_hatch_cert, #hatch date certainty date.fledge = TCL$Dt_flg, #fledging date dt.fledge.cert = TCL$Dt_flg_cert, #fledging date certainty # Fecundity variables DFE = TCL$DFE, #laying date = date of first egg PSC = TCL$Prob2ndClutch, # probable 2nd clutch based on DFE, hatchdate, datefound, stagefound - see below CS = TCL$Cl_size, #clutch size Cert = TCL$Cl_cert, #clutch size certainty Flgcert = TCL$Fate_cert, #fate certainty FledgeQ = TCL$Flg, #fledged? (Yes/No) exp_days = TCL$exp_days, # exposure days, for Daily Nest Survival model (nest-level) NoFledge = TCL$N_flg, #Brood Size Rsn_fail = TCL$Rsn_fail, #why did the nest fail # Nest-level extrinsic predictors CavSize = TCL$CavSizeEdworthy, # Cavity size (S, M, L) by definition in Table 1 of Edworthy et al 2019 Ht = TCL$Cav_ht_abv_grnd, #cavity height above ground # Site-level extrinsic predictors TAHU = TCL$TAHU, #Squirrels/ha from point counts (predators) MOCH = TCL$MOCH, #Mountain Chickadees/ha from point counts (conspecifics) BeetlePHA = TCL$beetle.ha, # site*yr-level beetle (beetle-infected pines/ha; spring food) LepTrees = TCL$LepTreesPerHa, # estimated lepidopteran-infected trees/ha (includes Aspen defoliators + Western Spruce Budworm; spring food) # Study area-level extrinsic predictors BeetleSA = TCL$BeetleSA, # Beetle-infected pines/ha across whole study area (winter food) SpringTemp = TCL$Temp, # weather; study area*Yr level StormsJanApr = TCL$StormsJanApr, StormsMayJune=TCL$StormsMayJune, WinterSnow = TCL$WinterSnow, # mean snow depth, 1 Jan to 15 March Templag2 = TCL$Templag2, Rainlag2=TCL$Rainlag2, SpringSnow = TCL$SnowDepthJanApr, #mean snow depth, 1 Jan to 30 Apr # Nests per ha FNestsPha = TCL$FrstNestPha, AllnestsPha=TCL$AllnestsPha) #Templag2 = Mean temperature from April1 through Aug30 2 years prior #Rainlag2 = Sum annual rainfall Sept1 in year x-3 through Aug30yearx-2, where x is MOCH$Year MOCH$Site <- as.factor(MOCH$Site) MOCH$Year <- as.factor(MOCH$Year) MOCH$Cert<-as.factor(MOCH$Cert) MOCH$CavSize<-as.factor(MOCH$CavSize) MOCH$PSC<-as.factor(MOCH$PSC) #reduce Cavity size to 2 categories (small, medium/large) MOCH$CavSize2<-ifelse(MOCH$CavSize=="M", "ML", ifelse(MOCH$CavSize=="L", "ML", ifelse(MOCH$CavSize=="S", "S", NA))) MOCH$CavSize2<-as.factor(MOCH$CavSize2) #make sure the variables below were brought in numeric #summary(MOCH$DFE) #summary(MOCH$CS) #summary(MOCH$NoFledge) #summary(MOCH$MOCH) #summary(MOCH$TAHU) #summary(MOCH$BeetleSA) #summary(MOCH$SpringTemp) #summary(MOCH$StormsJanApr) #summary(MOCH$StormsMayJune) #summary(MOCH$WinterSnow) #summary(MOCH$Ht) ######################################################################################## ############### Results Part I. Patterns of variation in fecundity ################# ################### (Descriptive) ############################## ######################################################################################## MOCH4<-MOCH MOCH4<-MOCH4[!is.na(MOCH4$TAHU),] MOCH4<-MOCH4[!is.na(MOCH4$MOCH),] MOCH4<-MOCH4[!is.na(MOCH4$BeetleSA),] MOCH4<-MOCH4[!is.na(MOCH4$BeetlePHA),] MOCH4<-MOCH4[!is.na(MOCH4$LepTrees),] #Incubation Period MOCHinc<-MOCH MOCHinc<-subset(MOCHinc, MOCHinc$dt.hatch.cert=="Y") MOCHinc<-MOCHinc[!is.na(MOCHinc$dt.inc.init),] MOCHinc<-subset(MOCHinc, MOCHinc$dt.inc.init!="U") MOCHinc<-droplevels(MOCHinc) MOCHinc$dt.inc.init<-as.character(MOCHinc$dt.inc.init) MOCHinc$dt.inc.init<-as.numeric(MOCHinc$dt.inc.init) summary(MOCHinc$dt.inc.init) summary(MOCHinc$date.hatch) MOCHinc$incperiod<-MOCHinc$date.hatch-MOCHinc$dt.inc.init summary(MOCHinc$incperiod) sd(MOCHinc$incperiod) #Nestling Period MOCHnp<-MOCH MOCHnp<-subset(MOCHnp, MOCHnp$dt.hatch.cert=="Y") MOCHnp<-MOCHnp[!is.na(MOCHnp$date.fledge),] MOCHnp<-subset(MOCHnp, MOCHnp$dt.fledge.cert == "1") MOCHnp<-MOCHnp[!is.na(MOCHnp$date.hatch),] MOCHnp<-droplevels(MOCHnp) MOCHnp$date.fledge<-as.character(MOCHnp$date.fledge) MOCHnp$date.fledge<-as.numeric(MOCHnp$date.fledge) summary(MOCHnp$date.fledge) summary(MOCHnp$date.hatch) MOCHnp$nestlingperiod<-MOCHnp$date.fledge-MOCHnp$date.hatch summary(MOCHnp$nestlingperiod) sd(MOCHnp$nestlingperiod) # Figure 2. Annual variation in Fecundity #Make a 4 panel plot of fecundity variables vs Year, with raw data (points), annual means (circles), and mean of annual means (horizontal line) MOCH.plot<-MOCH MOCH.plot<-MOCH.plot[!is.na(MOCH.plot$PSC),] summary(MOCH.plot$PSC) MOCH.plot$Year.n<-as.character(MOCH.plot$Year) MOCH.plot$Year.n<-as.numeric(MOCH.plot$Year.n) MOCH.plot$Year.n<-ifelse(MOCH.plot$PSC=="1", MOCH.plot$Year.n+0.1, MOCH.plot$Year.n) summary(MOCH.plot$Year.n) table(MOCH.plot$Year.n) dfe.plot<-MOCH.plot[!is.na(MOCH.plot$DFE),] dfe.plot$DFEann<-ifelse(dfe.plot$Year.n=="2000", mean(dfe.plot$DFE[dfe.plot$Year.n=="2000"]), ifelse(dfe.plot$Year.n=="2001", mean(dfe.plot$DFE[dfe.plot$Year.n=="2001"]), ifelse(dfe.plot$Year.n=="2002", mean(dfe.plot$DFE[dfe.plot$Year.n=="2002"]), ifelse(dfe.plot$Year.n=="2003", mean(dfe.plot$DFE[dfe.plot$Year.n=="2003"]), ifelse(dfe.plot$Year.n=="2004", mean(dfe.plot$DFE[dfe.plot$Year.n=="2004"]), ifelse(dfe.plot$Year.n=="2005", mean(dfe.plot$DFE[dfe.plot$Year.n=="2005"]), ifelse(dfe.plot$Year.n=="2005.1", mean(dfe.plot$DFE[dfe.plot$Year.n=="2005.1"]), ifelse(dfe.plot$Year.n=="2006", mean(dfe.plot$DFE[dfe.plot$Year.n=="2006"]), ifelse(dfe.plot$Year.n=="2007", mean(dfe.plot$DFE[dfe.plot$Year.n=="2007"]), ifelse(dfe.plot$Year.n=="2007.1", mean(dfe.plot$DFE[dfe.plot$Year.n=="2007.1"]), ifelse(dfe.plot$Year.n=="2008", mean(dfe.plot$DFE[dfe.plot$Year.n=="2008"]), ifelse(dfe.plot$Year.n=="2009", mean(dfe.plot$DFE[dfe.plot$Year.n=="2009"]), ifelse(dfe.plot$Year.n=="2010", mean(dfe.plot$DFE[dfe.plot$Year.n=="2010"]), ifelse(dfe.plot$Year.n=="2011", mean(dfe.plot$DFE[dfe.plot$Year.n=="2011"]), NA)))))))))))))) MeanDFE.df<-data.frame(Year.n = c(2000:2011)) #overall mean DFE for early nests MeanDFE.df$meanDFE<-(c(mean(dfe.plot$DFE[dfe.plot$Year.n==2000]), mean(dfe.plot$DFE[dfe.plot$Year.n==2001]), mean(dfe.plot$DFE[dfe.plot$Year.n==2002]), mean(dfe.plot$DFE[dfe.plot$Year.n==2003]), mean(dfe.plot$DFE[dfe.plot$Year.n==2004]), mean(dfe.plot$DFE[dfe.plot$Year.n==2005]), mean(dfe.plot$DFE[dfe.plot$Year.n==2006]), mean(dfe.plot$DFE[dfe.plot$Year.n==2007]), mean(dfe.plot$DFE[dfe.plot$Year.n==2008]), mean(dfe.plot$DFE[dfe.plot$Year.n==2009]), mean(dfe.plot$DFE[dfe.plot$Year.n==2010]), mean(dfe.plot$DFE[dfe.plot$Year.n==2011]))) MeanOverallDFE<-mean(MeanDFE.df$meanDFE) MeanOverallDFE CS.plot<-MOCH.plot CS.plot<-CS.plot[!is.na(CS.plot$CS),] summary(CS.plot$CS) CS.plot$CSann<-ifelse(CS.plot$Year.n=="2000", mean(CS.plot$CS[CS.plot$Year.n=="2000"]), ifelse(CS.plot$Year.n=="2001", mean(CS.plot$CS[CS.plot$Year.n=="2001"]), ifelse(CS.plot$Year.n=="2002", mean(CS.plot$CS[CS.plot$Year.n=="2002"]), ifelse(CS.plot$Year.n=="2003", mean(CS.plot$CS[CS.plot$Year.n=="2003"]), ifelse(CS.plot$Year.n=="2004", mean(CS.plot$CS[CS.plot$Year.n=="2004"]), ifelse(CS.plot$Year.n=="2005", mean(CS.plot$CS[CS.plot$Year.n=="2005"]), ifelse(CS.plot$Year.n=="2005.1", mean(CS.plot$CS[CS.plot$Year.n=="2005.1"]), ifelse(CS.plot$Year.n=="2006", mean(CS.plot$CS[CS.plot$Year.n=="2006"]), ifelse(CS.plot$Year.n=="2007", mean(CS.plot$CS[CS.plot$Year.n=="2007"]), ifelse(CS.plot$Year.n=="2007.1", mean(CS.plot$CS[CS.plot$Year.n=="2007.1"]), ifelse(CS.plot$Year.n=="2008", mean(CS.plot$CS[CS.plot$Year.n=="2008"]), ifelse(CS.plot$Year.n=="2009", mean(CS.plot$CS[CS.plot$Year.n=="2009"]), ifelse(CS.plot$Year.n=="2010", mean(CS.plot$CS[CS.plot$Year.n=="2010"]), ifelse(CS.plot$Year.n=="2011", mean(CS.plot$CS[CS.plot$Year.n=="2011"]), NA)))))))))))))) MeanCS.df<-data.frame(Year.n = c(2000:2011)) #overall mean CS for early nests MeanCS.df$meanCS<-(c(mean(CS.plot$CS[CS.plot$Year.n==2000]), mean(CS.plot$CS[CS.plot$Year.n==2001]), mean(CS.plot$CS[CS.plot$Year.n==2002]), mean(CS.plot$CS[CS.plot$Year.n==2003]), mean(CS.plot$CS[CS.plot$Year.n==2004]), mean(CS.plot$CS[CS.plot$Year.n==2005]), mean(CS.plot$CS[CS.plot$Year.n==2006]), mean(CS.plot$CS[CS.plot$Year.n==2007]), mean(CS.plot$CS[CS.plot$Year.n==2008]), mean(CS.plot$CS[CS.plot$Year.n==2009]), mean(CS.plot$CS[CS.plot$Year.n==2010]), mean(CS.plot$CS[CS.plot$Year.n==2011]))) MeanOverallCS<-mean(MeanCS.df$meanCS) MeanOverallCS NF.plot<-MOCH.plot NF.plot<-NF.plot[!is.na(NF.plot$NoFledge),] NF.plot<-subset(NF.plot, NF.plot$NoFledge>0) #restrict to nests that were successful summary(NF.plot$NoFledge) NF.plot$NFann<-ifelse(NF.plot$Year.n=="2000", mean(NF.plot$NoFledge[NF.plot$Year.n=="2000"]), ifelse(NF.plot$Year.n=="2001", mean(NF.plot$NoFledge[NF.plot$Year.n=="2001"]), ifelse(NF.plot$Year.n=="2002", mean(NF.plot$NoFledge[NF.plot$Year.n=="2002"]), ifelse(NF.plot$Year.n=="2003", mean(NF.plot$NoFledge[NF.plot$Year.n=="2003"]), ifelse(NF.plot$Year.n=="2004", mean(NF.plot$NoFledge[NF.plot$Year.n=="2004"]), ifelse(NF.plot$Year.n=="2005", mean(NF.plot$NoFledge[NF.plot$Year.n=="2005"]), ifelse(NF.plot$Year.n=="2005.1", mean(NF.plot$NoFledge[NF.plot$Year.n=="2005.1"]), ifelse(NF.plot$Year.n=="2006", mean(NF.plot$NoFledge[NF.plot$Year.n=="2006"]), ifelse(NF.plot$Year.n=="2007", mean(NF.plot$NoFledge[NF.plot$Year.n=="2007"]), ifelse(NF.plot$Year.n=="2007.1", mean(NF.plot$NoFledge[NF.plot$Year.n=="2007.1"]), ifelse(NF.plot$Year.n=="2008", mean(NF.plot$NoFledge[NF.plot$Year.n=="2008"]), ifelse(NF.plot$Year.n=="2009", mean(NF.plot$NoFledge[NF.plot$Year.n=="2009"]), ifelse(NF.plot$Year.n=="2010", mean(NF.plot$NoFledge[NF.plot$Year.n=="2010"]), ifelse(NF.plot$Year.n=="2011", mean(NF.plot$NoFledge[NF.plot$Year.n=="2011"]), NA)))))))))))))) MeanNF.df<-data.frame(Year.n = c(2000:2011)) #overall mean NF (brood size) for early nests MeanNF.df$meanNF<-(c(mean(NF.plot$NF[NF.plot$Year.n==2000]), mean(NF.plot$NF[NF.plot$Year.n==2001]), mean(NF.plot$NF[NF.plot$Year.n==2002]), mean(NF.plot$NF[NF.plot$Year.n==2003]), mean(NF.plot$NF[NF.plot$Year.n==2004]), mean(NF.plot$NF[NF.plot$Year.n==2005]), mean(NF.plot$NF[NF.plot$Year.n==2006]), mean(NF.plot$NF[NF.plot$Year.n==2007]), mean(NF.plot$NF[NF.plot$Year.n==2008]), mean(NF.plot$NF[NF.plot$Year.n==2009]), mean(NF.plot$NF[NF.plot$Year.n==2010]), mean(NF.plot$NF[NF.plot$Year.n==2011]))) MeanOverallNF<-mean(MeanNF.df$meanNF) MeanOverallNF #nests survived vs fledged NS.plot<-MOCH.plot NS.plot$Flgcert<-as.character(NS.plot$Flgcert) NS.plot$FledgeQ<-as.character(NS.plot$FledgeQ) NS.plot$Flgcert<-ifelse(NS.plot$FledgeQ=="N", "know", NS.plot$Flgcert) NS.plot$FledgeQ<-ifelse(NS.plot$Flgcert=="POSS", "U", NS.plot$FledgeQ) NS.plot <- subset(NS.plot, NS.plot$FledgeQ != "U") NS.plot$FledgeQQ <- ifelse(NS.plot$FledgeQ=="Y", 1, 0) NS.1<-data.frame(table(NS.plot$Year.n, NS.plot$FledgeQQ)) colnames(NS.1) <- c("Year.f", "Fledged", "Count") NS.1$Year.n<-as.character (NS.1$Year.f) NS.1$Year.n<-as.numeric(NS.1$Year.n) NS.1$Year.n2<-ifelse(NS.1$Fledged == 0, NS.1$Year.n+.05, NS.1$Year.n) NS.1$PSC<-c(0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0) # Plot the raw data and means in base R dev.off() par(mfrow=c(4,1), oma=c(1.5,2,3,1), cex.lab=1.2, mar=c(4,3,0.1,0.5)) #A: Laying date plot(dfe.plot$Year.n, dfe.plot$DFE, col=ifelse(dfe.plot$PSC==0,"black","blue"), cex = 0.3, las =1, xlab="", ylab="") abline(h=MeanOverallDFE, col='grey') points(dfe.plot$Year.n, dfe.plot$DFEann, col=ifelse(dfe.plot$PSC==0,"black","blue"), cex=2.2, pch=21) mtext("Laying date", side =2, line = 3, cex = 1) mtext("A", side =3, line=0, cex=0.8, adj=0.01) #B: Clutch size plot(CS.plot$Year.n, jitter(CS.plot$CS, factor=0.7), col=ifelse(CS.plot$PSC==0,"black","blue"), cex = 0.3, las =1, xlab="", ylab="", ylim=c(2,10)) abline(h=MeanOverallCS, col='grey') points(CS.plot$Year.n, CS.plot$CSann, col=ifelse(CS.plot$PSC==0,"black","blue"), cex=2.2, pch=21) mtext("Clutch size", side =2, line = 3, cex = 1) mtext("B", side =3, line=0, cex=0.8, adj=0.01) #C: Survived vs failed plot(NS.1$Year.n2, NS.1$Count, col = ifelse(NS.1$PSC==0,"black","blue"), las=1, xlab = "", ylab= "", bg=ifelse(NS.1$Fledged==0,ifelse(NS.1$PSC==0, "black", "blue"), NA), pch=ifelse(NS.1$Fledged==1, 21, 25)) mtext("Number of nests", side =2, line = 3, cex = 1) mtext("C", side =3, line=0, cex=0.8, adj=0.01) #D: Brood size plot(NF.plot$Year.n, jitter(NF.plot$NoFledge, factor=0.7), col=ifelse(NF.plot$PSC==0,"black","blue"), cex = 0.3, las =1, xlab="", ylab="", xlim=c(2000, 2011), ylim=c(0,10)) abline(h=MeanOverallNF, col='grey') points(NF.plot$Year.n, NF.plot$NFann, col=ifelse(NF.plot$PSC==0,"black","blue"), cex=2.2, pch=21) mtext("Brood size", side =2, line = 3, cex = 1) mtext("Year", side =1, line =3, cex =1) mtext("D", side =3, line=0, cex=0.8, adj=.01) # Some summary numbers for first section of Results #number of clutches laid after 3 June (day 154) in 2002 Day154<-subset(MOCH4, MOCH4$Year=="2002") Day154<-subset(Day154, Day154$DFE>154) length(Day154$DFE) summary(Day154$DFE) #laying and hatch dates mean(dfe.plot$DFE[dfe.plot$PSC==0]) #mean laying date of first nests length(dfe.plot$DFE[dfe.plot$PSC==0]) #n mean(dfe.plot$DFE[dfe.plot$PSC==1]) #mean laying date of second nests length(dfe.plot$DFE[dfe.plot$PSC==1]) #n hatch<-subset(MOCH4, !is.na(MOCH4)) hatch<-subset(hatch, hatch$dt.hatch.cert=="Y") hatchfirst<-subset(hatch, hatch$PSC==0) mean(hatchfirst$date.hatch) # mean hatch day of first clutches sd(hatchfirst$date.hatch) length(hatchfirst$date.hatch) #n hatch2nd<-subset(hatch, hatch$PSC==1) mean(hatch2nd$date.hatch) # mean hatch day of second clutches sd(hatch2nd$date.hatch) length(hatch2nd$date.hatch) #n #Clutch size mean(CS.plot$CS[CS.plot$PSC==0]) #mean clutch size of first nests length(CS.plot$CS[CS.plot$PSC==0]) #n sd(CS.plot$CS[CS.plot$PSC==0]) min(CS.plot$CS[CS.plot$PSC==0]) max(CS.plot$CS[CS.plot$PSC==0]) mean(CS.plot$CS[CS.plot$PSC==1]) #mean clutch size of late nests length(CS.plot$CS[CS.plot$PSC==1]) #n sd(CS.plot$CS[CS.plot$PSC==1]) min(CS.plot$CS[CS.plot$PSC==1]) max(CS.plot$CS[CS.plot$PSC==1]) #overall nest success sum(NS.1$Count[NS.1$Fledged==1]) # total number of successful nests sum(NS.1$Count) # total number of nests with known fate sum(NS.1$Count[NS.1$Fledged==1])/sum(NS.1$Count)*100 # % nest success length(NS.plot$FledgeQQ) length(NS.plot$FledgeQQ[NS.plot$FledgeQQ==1]) NS.first<-subset(NS.plot, NS.plot$PSC==0) length(NS.first$FledgeQQ) length(NS.first$FledgeQQ[NS.first$FledgeQQ==1]) NS.2nd<-subset(NS.plot, NS.plot$PSC==1) length(NS.2nd$FledgeQQ) length(NS.2nd$FledgeQQ[NS.2nd$FledgeQQ==1]) #percent success by yr NS.yr<-data.frame(Year=NS.1$Year.n[1:14]) NS.yr$Failed<-NS.1$Count[1:14] NS.yr$Fledged<-NS.1$Count[15:28] NS.yr$PctSuccess<-NS.yr$Fledged/(NS.yr$Fledged+NS.yr$Failed) # nest predation NS.plot$Rsn_fail<-as.factor(NS.plot$Rsn_fail) #Reason for nest failure - P = Predated, D= cavity Destroyed (eg. tree fell over, not if nest ripped out), C = Competitor, A = Abandoned, U = if reason for failure not known. NA if nest successfully fledged ("Flg"=Y) or fate is unknown. More than one code can be used, enter in aphabetical order. NS.fail<-subset(NS.plot, NS.plot$FledgeQQ==0) summary(NS.fail$Rsn_fail) length(NS.fail$FledgeQQ) #Brood size mean(NF.plot$NoFledge[NF.plot$PSC==0]) #mean brood size of first nests sd(NF.plot$NoFledge[NF.plot$PSC==0]) #sd brood size first nests min(NF.plot$NoFledge[NF.plot$PSC==0]) max(NF.plot$NoFledge[NF.plot$PSC==0]) mean(NF.plot$NoFledge[NF.plot$PSC==1]) #mean brood size of 2nd nests sd(NF.plot$NoFledge[NF.plot$PSC==1]) #sd brood size 2nd nests min(NF.plot$NoFledge[NF.plot$PSC==1]) max(NF.plot$NoFledge[NF.plot$PSC==1]) #Yearly mean MOCH point counts vs N nests dens<-data.frame(Year=c(2000:2011)) dens$MOCHpc<-c(mean(MOCH4$MOCH[MOCH4$Year==2000]),mean(MOCH4$MOCH[MOCH4$Year==2001]), mean(MOCH4$MOCH[MOCH4$Year==2002]),mean(MOCH4$MOCH[MOCH4$Year==2003]), mean(MOCH4$MOCH[MOCH4$Year==2004]),mean(MOCH4$MOCH[MOCH4$Year==2005]), mean(MOCH4$MOCH[MOCH4$Year==2006]),mean(MOCH4$MOCH[MOCH4$Year==2007]), mean(MOCH4$MOCH[MOCH4$Year==2008]),mean(MOCH4$MOCH[MOCH4$Year==2009]), mean(MOCH4$MOCH[MOCH4$Year==2010]),mean(MOCH4$MOCH[MOCH4$Year==2011])) dens$MOCHnests<-c(mean(MOCH4$FNestsPha[MOCH4$Year==2000]), mean(MOCH4$FNestsPha[MOCH4$Year==2001]), mean(MOCH4$FNestsPha[MOCH4$Year==2002]), mean(MOCH4$FNestsPha[MOCH4$Year==2003]), mean(MOCH4$FNestsPha[MOCH4$Year==2004]), mean(MOCH4$FNestsPha[MOCH4$Year==2005]), mean(MOCH4$FNestsPha[MOCH4$Year==2006]), mean(MOCH4$FNestsPha[MOCH4$Year==2007]), mean(MOCH4$FNestsPha[MOCH4$Year==2008]), mean(MOCH4$FNestsPha[MOCH4$Year==2009]), mean(MOCH4$FNestsPha[MOCH4$Year==2010]), mean(MOCH4$FNestsPha[MOCH4$Year==2011])) cor.test(dens$MOCHnests, dens$MOCHpc) ######################################################################################## ############### Results Part II. Drivers of variation in fecundity ################# ################### (Model selection) ############################## ######################################################################################## ############################################################# # LAYING DATE ~ ############################################################# #Laying date dataset MOCHdfe<-MOCH4 MOCHdfe<-MOCHdfe[!is.na(MOCHdfe$DFE),] MOCHdfe<-subset(MOCHdfe, MOCHdfe$PSC=="0") #restrict to probable first nests length(MOCHdfe$DFE) #n for the DFE model = 115 nests #scale Laying Date predictor variables to use MOCHdfe$SpringTemp.s<-scale(MOCHdfe$SpringTemp, center=TRUE, scale = TRUE) MOCHdfe$WinterSnow.s<-scale(MOCHdfe$WinterSnow, center=TRUE, scale=TRUE) MOCHdfe$SpringSnow.s<-scale(MOCHdfe$SpringSnow, center=TRUE, scale=TRUE) MOCHdfe$StormsJanApr.s<-scale(MOCHdfe$StormsJanApr, center=TRUE, scale = TRUE) #winter storms because laying decision comes before May-June storms MOCHdfe$Templag2.s<-scale(MOCHdfe$Templag2, center=TRUE, scale = TRUE) MOCHdfe$Rainlag2.s<-scale(MOCHdfe$Rainlag2, center=TRUE, scale = TRUE) MOCHdfe$MOCH.s<-scale(MOCHdfe$MOCH, center=TRUE, scale = TRUE) MOCHdfe$TAHU.s<-scale(MOCHdfe$TAHU, center=TRUE, scale = TRUE) MOCHdfe$LepTrees.s<-scale(MOCHdfe$LepTrees, center=TRUE, scale = TRUE) MOCHdfe$BeetleSA.s<-scale(MOCHdfe$BeetleSA, center=TRUE, scale = TRUE) MOCHdfe$BeetlePHA.s<-scale(MOCHdfe$BeetlePHA, center = TRUE, scale = TRUE) MOCHdfe$FNestsPha.s<-scale(MOCHdfe$FNestsPha, center = TRUE, scale = TRUE) MOCHdfe$AllnestsPha.s<-scale(MOCHdfe$AllnestsPha, center = TRUE, scale = TRUE) #Laying Date models dfe.reg1<-lme(DFE~SpringTemp.s*StormsJanApr.s + SpringSnow.s, random=~1|Year, data = MOCHdfe) #regional level, weather summary(dfe.reg1) intervals(dfe.reg1, which = "fixed") #95% confidence intervals ranef(dfe.reg1) plot(dfe.reg1) #looks fine dfe.reg2<-lme(DFE~ Templag2.s + Rainlag2.s + BeetleSA.s, random=~1|Year, data = MOCHdfe) #regional level, winter food summary(dfe.reg2) dfe.N<- lme(DFE~ 1 , random=~1|Year, data = MOCHdfe) #Null summary(dfe.N) plot(dfe.N) #looks fine model.sel(dfe.reg1, dfe.reg2, dfe.N) dfe.P<-lme(DFE~ TAHU.s + MOCH.s, random=~1|Year, data = MOCHdfe) #site level, predators and competitors summary(dfe.P) plot(dfe.P) #looks fine model.sel(dfe.P, dfe.N) ###################################################### #CLUTCH SIZE ~ ###################################################### MOCHcs<-MOCH4 MOCHcs<-subset(MOCHcs, MOCHcs$PSC=="0") #restrict to probable first nests MOCHcs<-MOCHcs[!is.na(MOCHcs$Cert), ] #clutchsize MOCHcs<-subset(MOCHcs, MOCHcs$Cert!="3") #clutchsize MOCHcs<-MOCHcs[!is.na(MOCHcs$DFE),] MOCHcs<-MOCHcs[!is.na(MOCHcs$Ht),] #exclude obs if no ht. MOCHcs<-MOCHcs[!is.na(MOCHcs$CavSize2),] length(MOCHcs$CS) # n = 78 observations #Scale Clutch Size Predictors: MOCHcs$DFE.s<-scale(MOCHcs$DFE, center=TRUE, scale = TRUE) MOCHcs$Ht.s<-scale(MOCHcs$Ht, center=TRUE, scale = TRUE) MOCHcs$SpringSnow.s<-scale(MOCHcs$SpringSnow, scale=TRUE, center=TRUE) MOCHcs$SpringTemp.s<-scale(MOCHcs$SpringTemp, center=TRUE, scale = TRUE) MOCHcs$StormsJanApr.s<-scale(MOCHcs$StormsJanApr, center =TRUE, scale=TRUE) #Winter Storms, as above MOCHcs$Templag2.s<-scale(MOCHcs$Templag2, center=TRUE, scale = TRUE) MOCHcs$Rainlag2.s<-scale(MOCHcs$Rainlag2, center=TRUE, scale = TRUE) MOCHcs$LepTrees.s<-scale(MOCHcs$LepTrees, center=TRUE, scale = TRUE) MOCHcs$BeetleSA.s<-scale(MOCHcs$BeetleSA, center=TRUE, scale = TRUE) MOCHcs$BeetlePHA.s<-scale(MOCHcs$BeetlePHA, center = TRUE, scale = TRUE) MOCHcs$MOCH.s<-scale(MOCHcs$MOCH, center=TRUE, scale = TRUE) MOCHcs$TAHU.s<-scale(MOCHcs$TAHU, center=TRUE, scale = TRUE) # Clutch Size models cs.reg1<-glmmTMB(CS ~ SpringTemp.s*StormsJanApr.s + + SpringSnow.s + (1| Year), family = compois, data=MOCHcs) summary(cs.reg1) cov2cor(vcov(cs.reg1)$cond) #check correlations #regional-level effects (weather only) cs.reg2<-glmmTMB(CS ~ Templag2.s + Rainlag2.s + BeetleSA.s + (1| Year), family = compois, data=MOCHcs) summary(cs.reg2) cov2cor(vcov(cs.reg2)$cond)#check correlations #regional-level effects (winter food with food-related weather) cs.N<- glmmTMB(CS ~ 1 + (1| Year), family = compois, data=MOCHcs) summary(cs.N) model.sel(cs.reg1, cs.reg2, cs.N) #neither regional model is better dAICc<2 but only SpringTemp*Storms is sig #clutch size increases with spring temp when low storms at regional level cs.FP<-glmmTMB(CS ~ LepTrees.s + BeetlePHA.s + MOCH.s + TAHU.s + (1| Year), family = compois, data=MOCHcs) summary(cs.FP) cov2cor(vcov(cs.FP)$cond)#check correlations model.sel(cs.FP, cs.N) cs.I<-glmmTMB(CS ~ DFE.s + Ht.s + CavSize2 + (1| Year), family = compois, data=MOCHcs) summary(cs.I) cov2cor(vcov(cs.I)$cond)#check correlations confint(cs.I) #95% confidence intervals to report in parameter table (paper) model.sel(cs.I, cs.N) ################################################### # Daily Nest Survival ~ ################################################### #Nest survival dataset MOCHns<-MOCH4 MOCHns$Flgcert<-as.character(MOCHns$Flgcert) MOCHns$FledgeQ<-as.character(MOCHns$FledgeQ) MOCHns$Flgcert<-ifelse(MOCHns$FledgeQ=="N", "know", MOCHns$Flgcert) MOCHns$FledgeQ<-ifelse(MOCHns$Flgcert=="POSS", "U", MOCHns$FledgeQ) MOCHns <- subset(MOCHns, MOCHns$FledgeQ != "U") MOCHns$FledgeQQ <- ifelse(MOCHns$FledgeQ=="Y", 1, 0) MOCHns <- MOCHns[!is.na(MOCHns$exp_days),] MOCHns <- MOCHns[!is.na(MOCHns$Ht),] MOCHns$trials<-1 MOCHns<-subset(MOCHns, MOCHns$exp_days!=0) MOCHns<-MOCHns[!is.na(MOCHns$CS),] MOCHns<-subset(MOCHns, MOCHns$PSC=="0") #Early nests only MOCHns<-MOCHns[!is.na(MOCHns$CavSize2),] MOCHns<-droplevels(MOCHns) length(MOCHns$CS) #logistic exposure link logexp <- function(exposure = 1) { linkfun <- function(mu) qlogis(mu^(1/exposure)) ## FIXME: is there some trick we can play here to allow ## evaluation in the context of the 'data' argument? linkinv <- function(eta) plogis(eta)^exposure mu.eta <- function(eta) exposure * plogis(eta)^(exposure-1) * .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats") valideta <- function(eta) TRUE link <- paste("logexp(", deparse(substitute(exposure)), ")", sep="") structure(list(linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta, valideta = valideta, name = link), class = "link-glm") } #Scale Nest Survival Variables MOCHns$CS.s<-scale(MOCHns$CS, scale=TRUE, center=TRUE) MOCHns$dfe.s<-scale(MOCHns$DFE, scale=TRUE, center=TRUE) MOCHns$Ht.s<-scale(MOCHns$Ht, scale=TRUE, center=TRUE) MOCHns$LepTrees.s<-scale(MOCHns$LepTrees, scale=TRUE, center=TRUE) MOCHns$BeetleSA.s<-scale(MOCHns$BeetleSA, scale=TRUE, center=TRUE) MOCHns$TAHU.s<-scale(MOCHns$TAHU, scale=TRUE, center=TRUE) MOCHns$MOCH.s<-scale(MOCHns$MOCH, scale = TRUE, center=TRUE) MOCHns$StormsMayJune.s<-scale(MOCHns$StormsMayJune, scale = TRUE, center=TRUE) MOCHns$SpringSnow.s<-scale(MOCHns$SpringSnow, scale=TRUE,center=TRUE) MOCHns$SpringTemp.s<-scale(MOCHns$SpringTemp, center=TRUE, scale = TRUE) MOCHns$Templag2.s<-scale(MOCHns$Templag2, center=TRUE, scale = TRUE) MOCHns$Rainlag2.s<-scale(MOCHns$Rainlag2, center=TRUE, scale = TRUE) MOCHns$BeetlePHA.s<-scale(MOCHns$BeetlePHA, center = TRUE, scale = TRUE) # Nest Survival Models #glm ns.N <- glm(FledgeQQ ~ 1, family= binomial(link = logexp(MOCHns$exp_days)), data=MOCHns) summary(ns.N) ns.I <- glm(FledgeQQ ~ CS.s + Ht.s + relevel(CavSize2, ref="S"), family= binomial(link = logexp(MOCHns$exp_days)), data=MOCHns) summary(ns.I) confint(ns.I, method = "Wald") model.sel(ns.I, ns.N) ns.reg1 <- glm(FledgeQQ ~ SpringTemp.s + SpringSnow.s*StormsMayJune.s, family= binomial(link = logexp(MOCHns$exp_days)), data=MOCHns) summary(ns.reg1) ns.reg2 <- glm(FledgeQQ ~ Templag2.s + Rainlag2.s + BeetleSA.s, family= binomial(link = logexp(MOCHns$exp_days)), data=MOCHns) summary(ns.reg2) model.sel(ns.reg1, ns.reg2, ns.N) # ns.FP <- glm(FledgeQQ ~ TAHU.s + LepTrees.s + BeetlePHA.s + MOCH.s, family= binomial(link = logexp(MOCHns$exp_days)), data=MOCHns) summary(ns.FP) model.sel(ns.N, ns.FP) ############################################################################### ### Brood Size ~ (for successful, probable 1st nests only) #################### ############################################################################### #DSR was for fail/fledge; this is for Brood Size in the nests that were successful (most nests) #Brood Size data set MOCHnf<-MOCH4 MOCHnf<-subset(MOCHnf, MOCHnf$NoFledge>0) #restrict to nests that were successful MOCHnf<-subset(MOCHnf, MOCHnf$PSC=="0") # Restrict to probable first clutches MOCHnf<-MOCHnf[!is.na(MOCHnf$NoFledge),] MOCHnf<-MOCHnf[!is.na(MOCHnf$Cert), ] #clutchsize MOCHnf<-subset(MOCHnf, MOCHnf$Cert!="3") #clutchsize MOCHnf<-MOCHnf[!is.na(MOCHnf$CavSize2), ] MOCHnf<-MOCHnf[!is.na(MOCHnf$Ht),] length(MOCHnf$NoFledge) # n #scale Brood Size variables to use MOCHnf$CS.s<-scale(MOCHnf$CS, center=TRUE, scale = TRUE) MOCHnf$Ht.s<-scale(MOCHnf$Ht, center = TRUE, scale = TRUE) MOCHnf$SpringTemp.s<-scale(MOCHnf$SpringTemp, center=TRUE, scale = TRUE) MOCHnf$SpringSnow.s<-scale(MOCHnf$SpringSnow, center=TRUE, scale=TRUE) MOCHnf$StormsMayJune.s<-scale(MOCHnf$StormsMayJune, center=TRUE, scale = TRUE) MOCHnf$Templag2.s<-scale(MOCHnf$Templag2, center=TRUE, scale = TRUE) MOCHnf$Rainlag2.s<-scale(MOCHnf$Rainlag2, center=TRUE, scale = TRUE) MOCHnf$BeetleSA.s<-scale(MOCHnf$BeetleSA, center=TRUE, scale = TRUE) MOCHnf$LepTrees.s<-scale(MOCHnf$LepTrees, center=TRUE, scale = TRUE) MOCHnf$BeetlePHA.s<-scale(MOCHnf$BeetlePHA, center = TRUE, scale = TRUE) MOCHnf$TAHU.s<-scale(MOCHnf$TAHU, center=TRUE, scale = TRUE) MOCHnf$MOCH.s<-scale(MOCHnf$MOCH, center=TRUE, scale = TRUE) #Brood Size models nf.N<-glmmTMB(NoFledge ~ 1 + (1| Year), family = compois, data=MOCHnf) summary(nf.N) nf.I<-glmmTMB(NoFledge ~ CS.s + Ht.s + relevel(CavSize2, ref="S") + (1| Year), family = compois, data=MOCHnf) summary(nf.I) confint(nf.I) cov2cor(vcov(nf.I)$cond)#check correlations model.sel(nf.I, nf.N) nf.reg1<-glmmTMB(NoFledge ~ SpringTemp.s*StormsMayJune.s + SpringSnow.s + (1| Year), family = compois, data=MOCHnf) summary(nf.reg1) cov2cor(vcov(nf.reg1)$cond)#check correlations nf.reg2<-glmmTMB(NoFledge ~ Templag2.s + Rainlag2.s + BeetleSA.s + (1| Year), family = compois, data=MOCHnf) summary(nf.reg2) cov2cor(vcov(nf.reg2)$cond)#check correlations model.sel(nf.reg1, nf.reg2, nf.N) # nf.FP<-glmmTMB(NoFledge ~ LepTrees.s + MOCH.s + TAHU.s + BeetlePHA.s + (1| Year), family = compois, data=MOCHnf) summary(nf.FP) cov2cor(vcov(nf.FP)$cond)#check correlations model.sel(nf.FP, nf.N) ################################################################# ########### Figs 3, 4, 5 and 6 ###### ########### Showing data and best fit models for: ###### ########### Laying date, Clutch Size, Nest Survival, & Brood Size.###### ########### Fig 5 (Nest survival) is made in Excel. ###### ################################################################# dev.off() ## ggplot customization: theme_set(theme_bw()) scale_colour_discrete <- function(...,palette="Set1") { scale_colour_brewer(...,palette=palette) } scale_colour_orig <- ggplot2::scale_colour_discrete scale_fill_discrete <- function(...,palette="Set1") { scale_fill_brewer(...,palette=palette) } ## to squash facets together ... zmargin <- theme(panel.spacing=grid::unit(0,"lines")) #par(mfrow=c(2,1),mar=c(5,9,8,1)*.7,las=1,xpd=TRUE) show_col(viridis_pal()(4)) viridis_pal()(4) #################################################################### ##### Fig 3 - Laying Date (Date of First Egg - DFE) ####### #################################################################### #obtain scaled values of 0 and 3 stormsJanApr: mean(MOCHdfe$StormsJanApr) sd(MOCHdfe$StormsJanApr) s0<-(0-mean(MOCHdfe$StormsJanApr))/(sd(MOCHdfe$StormsJanApr)) #0 storms s3<-(3-mean(MOCHdfe$StormsJanApr))/(sd(MOCHdfe$StormsJanApr)) #3 storms p1<-data.frame (SpringTemp.s = seq(from=min(MOCHdfe$SpringTemp.s), to=max(MOCHdfe$SpringTemp.s), by = 0.02)) p1$SpringTemp.backscaled=p1$SpringTemp.s*sd(MOCHdfe$SpringTemp)+mean(MOCHdfe$SpringTemp) p1$StormsJanApr.s<-s0 #0 storms p1$SpringSnow.s<-median(MOCHdfe$SpringSnow.s) p1$Year<-"2011" p1$preds0<-predict(dfe.reg1, newdata=p1, type="response") p2<-data.frame (SpringTemp.s = seq(from=min(MOCHdfe$SpringTemp.s), to=max(MOCHdfe$SpringTemp.s), by = 0.02)) p2$SpringTemp.backscaled=p2$SpringTemp.s*sd(MOCHdfe$SpringTemp)+mean(MOCHdfe$SpringTemp) p2$StormsJanApr.s<-s3 #3 storms p2$SpringSnow.s<-median(MOCHdfe$SpringSnow.s) p2$Year<-"2011" p2$preds3<-predict(dfe.reg1, newdata=p2, type="response") g2 <- ggplot(MOCHdfe,aes(x=SpringTemp,y=DFE, colour = StormsJanApr)) + geom_jitter(width = 0.2, height = 0.2) + labs(x = "Spring temperature (celsius)", y = "Laying date") + geom_point() + theme(aspect.ratio=3/3,panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + scale_colour_viridis(direction=-1, begin =0, end=1, name = "Storms") g2 + geom_line(data=p1, aes(x=SpringTemp.backscaled, y = preds0), colour= "#FDE725FF", size = 1) + #0 storms geom_line(data=p2, aes(x=SpringTemp.backscaled, y = preds3), colour= "#44015477", size = 1) #3 storms ############################################################ ############# Fig 4 Clutch size vs laying date ###### ############# col = SpringTemp, symbols = Storms ###### ############# in Panels by Year ###### ############################################################ MOCHcs$Storms<-as.factor(MOCHcs$StormsJanApr) MOCHcs$Temperature<-MOCHcs$SpringTemp g0 <- ggplot(MOCHcs,aes(x=DFE,y=CS, colour = Temperature, shape = Storms)) + geom_jitter(width = 0.1, height = 0.1) + xlim(125,179) + labs(x = "Laying date", y = "Clutch size") + geom_point() + scale_shape_manual(values=c(17, 15,18,16))+ theme(aspect.ratio=3/3) + scale_colour_viridis(direction=1, begin =0, end=1) + facet_wrap(.~Year) + ## 1 row of panels by roost zmargin ## squash together csNew<-data.frame(DFE.s=rep(seq(from=-2, to=4, by =0.01), each = 12)) csNew$DFE.backscaled<-csNew$DFE.s*sd(MOCHcs$DFE)+mean(MOCHcs$DFE) csNew$Year<-rep(2000:2011, length(csNew$DFE.s)/12) csNew$Ht.s<-median(MOCHcs$Ht.s) csNew$CavSize2<-"S" csNew$predicted<-predict(cs.I, newdata=csNew, type = "response") g0 + geom_point(data=csNew, aes(x=DFE.backscaled, y = predicted), colour= "black", shape=1, size = 0.01) ################################## ### Fig 5: Nest survival ######### ################################## #Panel A. Top Regional scale model. #pull b estimates from model summ.ns.reg2<-summary(ns.reg2) #model summary b0<-summ.ns.reg2$coefficients[1, 1] #intercept b0 b1<-summ.ns.reg2$coefficients[2, 1] #b1 templag - the variable of interest b2<-summ.ns.reg2$coefficients[3, 1] #b2 rainlag b3<-summ.ns.reg2$coefficients[4, 1] #b3 beetle #templag (scaled) is going to vary, while rainlag and beetleSA are held constant med.rainlag.s<-median(MOCHns$Rainlag2.s) #median scaled-rainlag med.beetle.s<-median(MOCHns$BeetleSA.s) #median scaled-beetleSA mean(MOCHns$Templag2) sd(MOCHns$Templag2) newdat<-data.frame(Templag2=seq(min(MOCHns$Templag2), max(MOCHns$Templag2), by = 0.02)) #original values of Templag2 for graph newdat$Templag2.s<-(newdat$Templag2-mean(MOCHns$Templag2))/sd(MOCHns$Templag2) #scaled values of Templag2 for calculation newdat$exponent<-exp(b0 + b1*newdat$Templag2.s + b2*med.rainlag.s + b3*med.beetle.s) newdat$DSR<-newdat$exponent/(1+newdat$exponent) #plot below #Panel B. Nest-scale model. #pull b estimates from model summ.ns.I<-summary(ns.I) #model summary b0<-summ.ns.I$coefficients[1, 1] #intercept b0 b1<-summ.ns.I$coefficients[2, 1] #b1 Clutch size - of interest b2<-summ.ns.I$coefficients[3, 1] #b2 Height - fix at median b3<-summ.ns.I$coefficients[4, 1] #b3 Cavsize - of interest (2 levels: S = 0, ML = 1) #CS (scaled) is going to vary, while Height held constant at median and Cavsize either 0 or 1 med.ht.s<-median(MOCHns$Ht.s) #median scaled-cavity-height newdat2<-data.frame(CS=seq(min(MOCHns$CS), max(MOCHns$CS), by = 0.02)) #original values of CS for graph newdat2$CS.s<-(newdat2$CS-mean(MOCHns$CS))/sd(MOCHns$CS) #scaled values of CS for calculation newdat2$exponent.sm<-exp(b0 + b1*newdat2$CS.s + b2*med.ht.s + b3*0) #for small cavities newdat2$exponent.ml<-exp(b0 + b1*newdat2$CS.s + b2*med.ht.s + b3*1) #for medium & large cavities newdat2$DSR.sm<-newdat2$exponent.sm/(1+newdat2$exponent.sm) #DSR for small cavities newdat2$DSR.ml<-newdat2$exponent.ml/(1+newdat2$exponent.ml) #DSR for med & large cavities #plot in Base R (not used) #par(mfcol=c(1,2)) #plot(newdat$Templag2, newdat$DSR, col="white", las=1, ylim = c(0.94, 1.0), xlab = "2-year lagged temperature (°C)") #lines(newdat$Templag2, newdat$DSR, lwd=2) #mtext("Daily nest survival", side=2, line = 3, cex=1.2) #plot(newdat2$CS, newdat2$DSR.sm, col="white", las=1, xlab = "Clutch size", ylim=c(0.94,1.0)) #lines(newdat2$CS, newdat2$DSR.sm, col = "brown", lwd=2) #lines(newdat2$CS, newdat2$DSR.ml, lwd=2) # plot in ggplot show_col(viridis_pal()(3)) viridis_pal()(3) g7 <- ggplot(newdat,aes(x=Templag2,y=DSR)) + ylim (0.94, 1.0) + labs(x = "2-year lagged temperature (°C)", y = "Daily nest survival") + theme(aspect.ratio=3/3,panel.grid.major = element_blank(), panel.grid.minor = element_blank()) g7.1<-g7 + geom_line(data=newdat, aes(x=Templag2, y = DSR), colour= "#21908CFF", size = 1) + labs(tag = "A") + theme(plot.tag.position = c(0.18, 1.04)) g8 <- ggplot(newdat2,aes(x=CS,y=DSR.sm)) + ylim (0.94, 1.0) + xlim(3,9) + labs(x = "Clutch size", y = "") + theme(aspect.ratio=3/3,panel.grid.major = element_blank(), panel.grid.minor = element_blank()) g8.1<- g8 + geom_line(data=newdat2, aes(x=CS, y = DSR.sm), colour= "#440154FF", size = 1) + geom_line(data=newdat2, aes(x=CS, y = DSR.ml), colour= "#FDE725FF", size =1) + labs(tag = "B") + theme(plot.tag.position = c(0.18, 1.04)) grid.arrange(g7.1, g8.1, nrow = 1) ##################################################### ### Fig 6 - Brood Size ####### ##################################################### g1 <- ggplot(MOCHnf,aes(x=CS,y=NoFledge, colour = DFE)) + geom_jitter(width = 0.2, height = 0.2) + labs(x = "Clutch size", y = "Brood size") + geom_point() + theme(aspect.ratio=3/3,panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + scale_colour_viridis(direction=-1, begin =0, end=1, name = "Laying date") + xlim(2.8,9.2) g1 nfNew<-data.frame(CS.s=rep(seq(from=-3, to=3, by =0.01), each = 11)) nfNew$CS.backscaled<-mean(MOCHnf$CS)+nfNew$CS.s*sd(MOCHnf$CS) nfNew$Year<-2004 nfNew$Ht.s<-median(MOCHnf$Ht.s) nfNew$CavSize2<-"S" nfNew$CavSize2<-as.factor(nfNew$CavSize2) nfNew$predicted<-predict(nf.I, newdata=nfNew, type = "response") g1 + geom_point(data=nfNew, aes(x=CS.backscaled, y = predicted), colour= "black", size = 0.01)