library(dplyr)
library(lme4)
library(lmerTest)
library(glmmTMB)
library(MuMIn)
library(arm)
library(ggplot2)
library(GGally)
library(grid)
library(gridExtra)
library(cowplot)
library(raster)
library(patchwork)

#Pull and format data
setwd("C:/Users/swi/Desktop/Revision_Repro/Final")
data <- read.csv("data.csv")
data <- data %>% select(-X) 

data$bird_id <- as.factor(data$bird_id)
data$nest_id <- as.factor(data$nest_id)
data$year <- as.factor(data$year)
data$migration <- as.factor(data$migration)
data$day_return <- as.numeric(data$day_return)
data$nr_fledglings <- as.numeric(data$nr_fledglings)
data$day_incubate <- as.numeric(data$day_incubate)
data$day_incubate <- as.numeric(data$day_incubate)
data$experience <- as.factor(data$experience)
data$couple <- as.factor(data$couple)

data$size.z <- scale(data$size)
data$elevation.z <- scale(data$elevation)
data$day_incubate.z <- scale(data$day_incubate)
data$day_return.z <- scale(data$day_return)


###Nest-site elevation###


##Mean values 

#Residents
mean(data[data$migration == "0",]$elevation)
sd(data[data$migration == "0",]$elevation)
range(data[data$migration == "0",]$elevation)

#Migrants
mean(data[data$migration == "1",]$elevation)
sd(data[data$migration == "1",]$elevation)
range(data[data$migration == "1",]$elevation)


##Model

#data$year <- relevel(data$year, ref = "2018") #change reference level of factors to extract contrasts. Default = 2017

modelev <- lmer(elevation ~ migration + size.z + year + experience + migration*year + (1|couple), data)


##Results 

#Effects table
nsim <- 10000
samp <- sim(modelev, n.sim = nsim)

effecte <- as.data.frame(fixef(modelev))
credIe <- as.data.frame(apply(samp@fixef, 2, quantile, prob = c(0.025,0.975)))
credIte <- as.data.frame(t(credIe))
colnames(credIte) <- rownames(credIe)
rownames(credIte) <- colnames(credIe)
outpute <- cbind(effecte,credIte)
outpute

#Conditional R^2
performance::r2(modelev)

#Random effect
attributes(VarCorr(modelev)$couple)$stddev #couple


###Timing of breeding (All birds)###


##Mean values 

#Residents
mean(data[data$migration == "0",]$day_incubate)
sd(data[data$migration == "0",]$day_incubate)
range(data[data$migration == "0",]$day_incubate)

#Migrants
mean(data[data$migration == "1",]$day_incubate)
sd(data[data$migration == "1",]$day_incubate)
range(data[data$migration == "1",]$day_incubate)


##Model

#data$year <- relevel(data$year, ref = "2018") #change reference level of factors to extract contrasts. Default = 2017

modtb <- lmer(day_incubate ~ migration + size.z + year + experience + elevation.z + migration*year + (1|couple/nest_id), data = data) 


##Results

#Effects table
nsim <- 10000
samp <- sim(modtb, n.sim = nsim)

effecttb <- as.data.frame(fixef(modtb))
credItb <- as.data.frame(apply(samp@fixef, 2, quantile, prob = c(0.025,0.975)))
credIttb <- as.data.frame(t(credItb))
colnames(credIttb) <- rownames(credItb)
rownames(credIttb) <- colnames(credItb)
outputtb <- cbind(effecttb,credIttb)
outputtb

#Conditional R^2
performance::r2(modtb)

#Random effect
attributes(VarCorr(modtb)$couple)$stddev #couple
attributes(VarCorr(modtb)$nest_id)$stddev #couple/nest_id

#Estimate difference in timing of breeding per 100 m increase in elevation

m <- mean(data$elevation)
s <- sd(data$elevation)

newdat <- expand.grid(migration = c("0","1"), size = mean(data$size.z), year = c("2017","2018","2019","2020"), experience = c("0","1"), elevation = seq(min(data$elevation.z), scale((min(data$elevation) + 100), center = m, scale = s), length=2))

elevation_unz <- (newdat$elevation * s) + m 
newdat <- cbind(newdat, elevation_unz)

Xmat <- model.matrix(~ migration + size + year + experience + elevation + migration*year, data = newdat)

b <- fixef(modtb)
newdat$fit <- Xmat %*% b
fitmat <- matrix(ncol = nsim, nrow = nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- Xmat %*% samp@fixef[i,]

fitmat <- as.data.frame(fitmat[c(1,17),])
fitmatdiff <- as.data.frame(t(apply(fitmat, 2, diff, 1)))

newdat$fit[17]-newdat$fit[1] #difference in timing of breeding (Julian days)
apply(fitmatdiff, 1, quantile, prob=c(0.025, 0.975)) #credible intervals of difference in timing of breeding


###Timing of breeding (Migrants only)###


##Model

#data$year <- relevel(data$year, ref = "2018") #change reference level of factors to extract contrasts. Default = 2017

datam <- filter(data, migration == "1") #migrants only

modtbm <- lmer(day_incubate ~ day_return.z + size.z + year + experience + elevation.z + (1|couple/nest_id), data = datam)


##Results

#Effects table
nsim <- 10000
samp <- sim(modtbm, n.sim = nsim)

effecttbm <- as.data.frame(fixef(modtbm))
credItbm <- as.data.frame(apply(samp@fixef, 2, quantile, prob = c(0.025,0.975)))
credIttbm <- as.data.frame(t(credItbm))
colnames(credIttbm) <- rownames(credItbm)
rownames(credIttbm) <- colnames(credItbm)
outputtbm <- cbind(effecttbm,credIttbm)
outputtbm

#Conditional R^2
performance::r2(modtbm)

#Random effect
attributes(VarCorr(modtbm)$couple)$stddev #couple
attributes(VarCorr(modtbm)$nest_id)$stddev #couple/nest_id

#Estimate difference in timing of breeding per 100 m increase in elevation

m <- mean(datam$day_return)
s <- sd(datam$day_return)

newdat <- expand.grid(day_return = seq(min(datam$day_return.z), scale((min(datam$day_return) + 10), center = m, scale = s), length=2), size = mean(data$size.z), year = c("2017","2018","2019","2020"), experience = c("0","1"), elevation = mean(data$elevation.z))

day_return_unz <- (newdat$day_return * s) + m 
newdat <- cbind(newdat, day_return_unz)

Xmat <- model.matrix(~ day_return + size + year + experience + elevation, data = newdat)

b <- fixef(modtbm)
newdat$fit <- Xmat %*% b
fitmat <- matrix(ncol = nsim, nrow = nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- Xmat %*% samp@fixef[i,]

fitmat <- as.data.frame(fitmat[c(1,2),])
fitmatdiff <- as.data.frame(t(apply(fitmat, 2, diff, 1)))

newdat$fit[2]-newdat$fit[1] #difference in timing of breeding (Julian days)
apply(fitmatdiff, 1, quantile, prob=c(0.025, 0.975)) #credible intervals of difference in timing of breeding


###Number of fledglings###


##Model

#data$year <- relevel(data$year, ref = "2018") #change reference level of factors to extract contrasts. Default = 2017

modnf <-glmmTMB(nr_fledglings ~ migration + size.z + year + experience + elevation.z + day_incubate.z + migration*year + migration*elevation.z + (1|couple), family = compois(link = "log"), data = data)


##Results

#Effects table
nsim <- 10000
pp <- fixef(modnf)$cond
vv <- vcov(modnf)$cond
samp <- mvrnorm(nsim, mu=pp, Sigma=vv)

effectnf <- fixef(modnf)
effectnf <- as.data.frame(effectnf$cond)
credInf <- as.data.frame(apply(samp, 2, quantile, prob = c(0.025,0.975)))
credItnf <- as.data.frame(t(credInf))
colnames(credItnf) <- rownames(credInf)
rownames(credItnf) <- colnames(credInf)
outputnf <- cbind(effectnf,credItnf)
outputnf

#Conditional (pseudo) R^2
performance::r2(modnf)

#Random effect
summary(modnf)[9] #couple

#Dispersion parameter
summary(modnf)[7]

#Estimate difference in number of fledglings per 100 m increase in elevation

m <- mean(data$elevation)
s <- sd(data$elevation)

newdat <- expand.grid(migration = c("0","1"), size = mean(data$size.z), year = c("2017","2018","2019", "2020"), experience = c("0","1"), elevation =  seq(min(data$elevation.z), scale((min(data$elevation) + 100), center = m, scale = s), length = 2), day_incubate = mean(data$day_incubate.z))

elevation_unz <- (newdat$elevation * s) + m 
newdat <- cbind(newdat, elevation_unz)

Xmat <- model.matrix(~ migration + size + year + experience + elevation + day_incubate + migration*year + migration*elevation, data = newdat)

b <- fixef(modnf)$cond
newdat$fit <- exp(Xmat %*% b)

fitmat <- matrix(ncol = nsim, nrow = nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- exp(Xmat %*% samp[i,])

fitmat <- as.data.frame(fitmat[c(4,20),])
fitmatdiff <- as.data.frame(t(apply(fitmat, 2, diff, 1)))

newdat$fit[18]-newdat$fit[2] #difference in fledglings
apply(fitmatdiff, 1, quantile, prob=c(0.025,0.975)) #credible intervals of difference in fledglings

#Estimate difference in number of fledglings between the earliest and latest breeders

m <- mean(data$day_incubate)
s <- sd(data$day_incubate)

newdat <- expand.grid(migration = c("0","1"), size = mean(data$size.z), year = c("2017","2018","2019", "2020"), experience = c("0","1"), elevation = mean(data$elevation.z), day_incubate = seq(min(data$day_incubate.z), max(data$day_incubate.z), length = 2))
 
day_incubate_unz <- (newdat$day_incubate * s) + m 
newdat <- cbind(newdat, day_incubate_unz)

Xmat <- model.matrix(~ migration + size + year + experience + elevation + day_incubate + migration*year + migration*elevation, data = newdat)

b <- fixef(modnf)$cond
newdat$fit <- exp(Xmat %*% b)

fitmat <- matrix(ncol = nsim, nrow = nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- exp(Xmat %*% samp[i,])

fitmat <- as.data.frame(fitmat[c(1,17),])
fitmatdiff <- as.data.frame(t(apply(fitmat, 2, diff,1)))
 
newdat$fit[17]-newdat$fit[1] #difference in fledglings
apply(fitmatdiff, 1, quantile, prob=c(0.025, 0.975)) #credible intervals of difference in fledglings


##Graphics

path <- "C:/Users/swi/Desktop/Revision_Repro/Final"
my_colours <- c("Orange", "Sky Blue")

#Migration x elevation

m <- mean(data$elevation)
s <- sd(data$elevation)

newdat <- expand.grid(migration = c("0","1"), size = mean(data$size.z), year = c("2017","2018","2019", "2020"), experience = c("0","1"), elevation =  seq(min(data$elevation.z), max(data$elevation.z), length = 600), day_incubate = mean(data$day_incubate.z))

elevation_unz <- (newdat$elevation * s) + m 
newdat <- cbind(newdat, elevation_unz)

Xmat <- model.matrix(~ migration + size + year + experience + elevation + day_incubate + migration*year + migration*elevation, data = newdat)

b <- fixef(modnf)$cond
newdat$fit <- exp(Xmat %*% b)
fitmat <- matrix(ncol = nsim, nrow = nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- exp(Xmat %*% samp[i,])

newdat$lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdat$upr <- apply(fitmat, 1, quantile, prob=0.975)
newdat$fit <- as.vector(newdat$fit)

#2017
newdat17 <- as.data.frame(newdat[newdat$year=="2017" & newdat$experience=="0",])
data17 <- as.data.frame(data[data$year=="2017" & data$experience=="0",])
names(data17)[names(data17) == 'nr_fledglings'] <- 'fit'
names(data17)[names(data17) == 'elevation'] <- 'elevation_unz'

p1 <- ggplot(newdat17, aes(x = elevation_unz, y = fit, group = migration)) + scale_colour_manual(values = my_colours, labels = c('Resident', 'Migrant')) + scale_fill_manual(values = my_colours, labels = c('Resident', 'Migrant')) + geom_jitter(data17, mapping = aes(color = migration)) + geom_line(stat = "identity", aes(color = migration)) + geom_ribbon(aes(ymin = lwr, ymax = upr, group = migration, color = migration, fill = migration), alpha = .5) + labs(title = "2017") + ylim(0,5) + theme_classic() + theme(strip.text.x = element_text(size = 16), strip.background = element_blank(), text = element_text(size = 13), axis.text.y.right = element_blank(), axis.ticks.y.right = element_blank(), axis.line.y.right = element_blank(), legend.title=element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), plot.margin = unit(c(0.5,.25,.5,.75), 'cm')) + theme(legend.position = "bottom")

#2018
newdat18 <- as.data.frame(newdat[newdat$year=="2018" & newdat$experience=="0",])
data18 <- as.data.frame(data[data$year=="2018" & data$experience=="0",])
names(data18)[names(data18) == 'nr_fledglings'] <- 'fit'
names(data18)[names(data18) == 'elevation'] <- 'elevation_unz'

p2 <- ggplot(newdat18, aes(x = elevation_unz, y = fit, group = migration)) + scale_colour_manual(values = my_colours, labels = c('Resident', 'Migrant')) + scale_fill_manual(values = my_colours, labels = c('Resident', 'Migrant')) + geom_jitter(data18, mapping = aes(color = migration)) + geom_line(stat = "identity", aes(color = migration)) + geom_ribbon(aes(ymin = lwr, ymax = upr, group = migration, color = migration, fill = migration), alpha = .5) + labs(title = "2018") + ylim(0,5) + theme_classic() + theme(strip.text.x = element_text(size = 16), strip.background = element_blank(), text = element_text(size = 14), axis.text.y.right = element_blank(), axis.ticks.y.right = element_blank(), axis.line.y.right = element_blank(), legend.title=element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), plot.margin = unit(c(0.5,.25,.5,.75), 'cm')) + theme(legend.position = "none")

#2019
newdat19 <- as.data.frame(newdat[newdat$year=="2019" & newdat$experience=="0",])
data19 <- as.data.frame(data[data$year=="2019" & data$experience=="0",])
names(data19)[names(data19) == 'nr_fledglings'] <- 'fit'
names(data19)[names(data19) == 'elevation'] <- 'elevation_unz'

p3 <- ggplot(newdat19, aes(x = elevation_unz, y = fit, group = migration)) + scale_colour_manual(values = my_colours, labels = c('Resident', 'Migrant')) + scale_fill_manual(values = my_colours, labels = c('Resident', 'Migrant')) + geom_jitter(data19, mapping = aes(color = migration)) + geom_line(stat = "identity", aes(color = migration)) + geom_ribbon(aes(ymin = lwr, ymax = upr, group = migration, color = migration, fill = migration), alpha = .5) + labs(title = "2019") + ylim(0,5) + theme_classic() + theme(strip.text.x = element_text(size = 16), strip.background = element_blank(), text = element_text(size = 14), axis.text.y.right = element_blank(), axis.ticks.y.right = element_blank(), axis.line.y.right = element_blank(), legend.title=element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), plot.margin = unit(c(0.5,.25,.5,.75), 'cm')) + theme(legend.position = "none")

#2020
newdat20 <- as.data.frame(newdat[newdat$year=="2020" & newdat$experience=="0",])
data20 <- as.data.frame(data[data$year=="2020" & data$experience=="0",])
names(data20)[names(data20) == 'nr_fledglings'] <- 'fit'
names(data20)[names(data20) == 'elevation'] <- 'elevation_unz'

p4 <- ggplot(newdat20, aes(x = elevation_unz, y = fit, group = migration)) + scale_colour_manual(values = my_colours, labels = c('Resident', 'Migrant')) + scale_fill_manual(values = my_colours, labels = c('Resident', 'Migrant')) + geom_jitter(data20, mapping = aes(color = migration)) + geom_line(stat = "identity", aes(color = migration)) + geom_ribbon(aes(ymin = lwr, ymax = upr, group = migration, color = migration, fill = migration), alpha = .5) + labs(title = "2020") + ylim(0,5) + theme_classic() + theme(strip.text.x = element_text(size = 16), strip.background = element_blank(), text = element_text(size = 14), axis.text.y.right = element_blank(), axis.ticks.y.right = element_blank(), axis.line.y.right = element_blank(), legend.title=element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), plot.margin = unit(c(0.5,.25,1.5,.75), 'cm')) + theme(legend.position = "none")

legend1 <- get_legend(p1)
p1 <- p1 + theme(legend.position = "none")

#Landscape projection

newdat1 <- newdat %>% group_by(year, experience, elevation_unz) %>% mutate(RRO = diff(desc(fit)))
newdat2 <- newdat1 %>% distinct(size, year, experience, day_incubate, RRO, .keep_all = TRUE)

#2017
newdat3 <- newdat2[newdat2$year == "2017",]
newdat4 <- newdat3[newdat3$experience == "0",]
newdat5 <- newdat4[,c(7,11)]
newdat5$elevation_unz <- round(newdat5$elevation_unz)
newdat6 <- newdat5 %>% distinct(elevation_unz, .keep_all = TRUE)

DEM <- raster("Ch2_studyarea_elevation.gpkg", RAT = TRUE) #layer available from swisstopo (https://www.swisstopo.admin.ch/en/height-model-swissalti3d)

DEMagg <- aggregate(DEM, mean, fact = 500)
DEMaggmask <- mask(DEMagg, mask = min(newdat$elevation_unz) < DEMagg$Height & DEMagg$Height < max(newdat$elevation_unz), maskvalue = NA)
DEMaggmask <- round(DEMaggmask)
DEMdf <- as.data.frame(DEMaggmask, xy = TRUE)
DEMdf$order <- seq(1, nrow(DEMdf))
DEMdfattr <- merge(DEMdf, newdat6, by.x = "layer", by.y = "elevation_unz", all.x = TRUE)
DEMdfattr17 <- arrange(DEMdfattr, order)
DEMdfattr1 <- DEMdfattr17[,c(2:3,5)]
DEMrastRRO17 <- rasterFromXYZ(DEMdfattr1)

#2018
newdat3 <- newdat2[newdat2$year == "2018",]
newdat4 <- newdat3[newdat3$experience == "0",]
newdat5 <- newdat4[,c(7,11)]
newdat5$elevation_unz <- round(newdat5$elevation_unz)
newdat6 <- newdat5 %>% distinct(elevation_unz, .keep_all = TRUE)

DEMdfattr <- merge(DEMdf, newdat6, by.x = "layer", by.y = "elevation_unz", all.x = TRUE)
DEMdfattr18 <- arrange(DEMdfattr, order)
DEMdfattr1 <- DEMdfattr18[,c(2:3,5)]
DEMrastRRO18 <- rasterFromXYZ(DEMdfattr1)

#2019
newdat3 <- newdat2[newdat2$year == "2019",]
newdat4 <- newdat3[newdat3$experience == "0",]
newdat5 <- newdat4[,c(7,11)]
newdat5$elevation_unz <- round(newdat5$elevation_unz)
newdat6 <- newdat5 %>% distinct(elevation_unz, .keep_all = TRUE)
DEMdfattr <- merge(DEMdf, newdat6, by.x = "layer", by.y = "elevation_unz", all.x = TRUE)
DEMdfattr19 <- arrange(DEMdfattr, order)
DEMdfattr1 <- DEMdfattr19[,c(2:3,5)]
DEMrastRRO19 <- rasterFromXYZ(DEMdfattr1)

#2020
newdat3 <- newdat2[newdat2$year == "2020",]
newdat4 <- newdat3[newdat3$experience == "0",]
newdat5 <- newdat4[,c(7,11)]
newdat5$elevation_unz <- round(newdat5$elevation_unz)
newdat6 <- newdat5 %>% distinct(elevation_unz, .keep_all = TRUE)

DEMdfattr <- merge(DEMdf, newdat6, by.x = "layer", by.y = "elevation_unz", all.x = TRUE)
DEMdfattr20 <- arrange(DEMdfattr, order)
DEMdfattr1 <- DEMdfattr20[,c(2:3,5)]
DEMrastRRO20 <- rasterFromXYZ(DEMdfattr1)

DEMRROstack <- stack(DEMrastRRO17, DEMrastRRO18, DEMrastRRO19, DEMrastRRO20)

#2017
p17 <- ggplot(DEMdfattr17, aes(x, y, fill = RRO)) + geom_raster()  + scale_fill_gradientn(name = "RRO", colours = c("Sky Blue", "white", "Orange"), values = scales::rescale(c(min(minValue(DEMRROstack)), -0.2, 0, 0.2, max(maxValue(DEMRROstack)))), na.value = "grey35", limits=c(min(minValue(DEMRROstack)), max(maxValue(DEMRROstack))), breaks = c(0.5, 0, -0.5, -1.0, -1.5, -2.0), labels = c("0.5", 0, -0.5, -1.0, -1.5, "-2.0"), guide = guide_colorbar(title.vjust = .9, barwidth = 6)) + theme_classic() + theme(axis.line = element_blank(), axis.text = element_blank(), axis.ticks=element_blank(), axis.title=element_blank(), panel.border = element_rect(color = "black", fill = NA), legend.text = element_text(size = 10), legend.title = element_text(size = 16)) + coord_equal() + theme(legend.position = "bottom")

#2018
p18 <- ggplot(DEMdfattr18, aes(x, y, fill = RRO)) + geom_raster()  + scale_fill_gradientn(name = "RRO", colours = c("Sky Blue", "white", "Orange"), values = scales::rescale(c(min(minValue(DEMRROstack)), -0.2, 0, 0.2, max(maxValue(DEMRROstack)))), na.value = "grey35", limits=c(min(minValue(DEMRROstack)), max(maxValue(DEMRROstack))), breaks = c(0.5, 0, -0.5, -1.0, -1.5, -2.0), labels = c("0.5", 0, -0.5, -1.0, -1.5, "-2.0"), guide = guide_colorbar(title.vjust = .9, barwidth = 6)) + theme_classic() + theme(axis.line = element_blank(), axis.text = element_blank(), axis.ticks=element_blank(), axis.title=element_blank(), panel.border = element_rect(color = "black", fill = NA), legend.text = element_text(size = 10), legend.title = element_text(size = 16)) + coord_equal() + theme(legend.position = "none")

#2019
p19 <- ggplot(DEMdfattr19, aes(x, y, fill = RRO)) + geom_raster()  + scale_fill_gradientn(name = "RRO", colours = c("Sky Blue", "white", "Orange"), values = scales::rescale(c(min(minValue(DEMRROstack)), -0.2, 0, 0.2, max(maxValue(DEMRROstack)))), na.value = "grey35", limits=c(min(minValue(DEMRROstack)), max(maxValue(DEMRROstack))), breaks = c(0.5, 0, -0.5, -1.0, -1.5, -2.0), labels = c("0.5", 0, -0.5, -1.0, -1.5, "-2.0"), guide = guide_colorbar(title.vjust = .9, barwidth = 6)) + theme_classic() + theme(axis.line = element_blank(), axis.text = element_blank(), axis.ticks=element_blank(), axis.title=element_blank(), panel.border = element_rect(color = "black", fill = NA), legend.text = element_text(size = 10), legend.title = element_text(size = 16)) + coord_equal() + theme(legend.position = "none")

#2020
p20 <- ggplot(DEMdfattr20, aes(x, y, fill = RRO)) + geom_raster()  + scale_fill_gradientn(name = "RRO", colours = c("Sky Blue", "white", "Orange"), values = scales::rescale(c(min(minValue(DEMRROstack)), -0.2, 0, 0.2, max(maxValue(DEMRROstack)))), na.value = "grey35", limits=c(min(minValue(DEMRROstack)), max(maxValue(DEMRROstack))), breaks = c(0.5, 0, -0.5, -1.0, -1.5, -2.0), labels = c("0.5", 0, -0.5, -1.0, -1.5, "-2.0"), guide = guide_colorbar(title.vjust = .9, barwidth = 6))  + theme_classic() + theme(axis.line = element_blank(), axis.text = element_blank(), axis.ticks=element_blank(), axis.title=element_blank(), panel.border = element_rect(color = "black", fill = NA), legend.text = element_text(size = 10), legend.title = element_text(size = 16), plot.margin = unit(c(.2,.15,1.25,.15), 'cm')) + coord_equal() + theme(legend.position = "none")

legend2 <- get_legend(p17)
p17 <- p17 + theme(legend.position = "none")

fig3 <- grid.arrange(p1,p17,p2,p18,p3,p19,p4,p20,legend1,legend2, ncol = 2, nrow = 5, widths = c(2,1), heights = c(1,1,1,1.2,0.2), bottom = textGrob("Elevation (m.a.s.l.)", gp=gpar(fontsize=15, face = "bold"), x = .34, y = 4), left = textGrob("Fledgling Number", gp=gpar(fontsize=15, face = "bold"), x = 1, y = 0.55, rot = 90))

ggsave(filename = "Figure_3.jpeg", plot = fig3, path = path, width = 18, height = 24, units = "cm", device='jpeg', dpi=600)

#Timing of breeding

m <- mean(data$day_incubate)
s <- sd(data$day_incubate)

newdat <- expand.grid(migration = c("0","1"), size = mean(data$size.z), year = c("2017","2018","2019", "2020"), experience = c("0","1"), elevation = mean(data$elevation.z), day_incubate = seq(min(data$day_incubate.z), max(data$day_incubate.z), length = 100))

day_incubate_unz <- (newdat$day_incubate * s) + m 
newdat <- cbind(newdat, day_incubate_unz)

Xmat <- model.matrix(~ migration + size + year + experience + elevation + day_incubate + migration*year + migration*elevation, data = newdat)

b <- fixef(modnf)$cond
newdat$fit <- exp(Xmat %*% b)
fitmat <- matrix(ncol = nsim, nrow = nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- exp(Xmat %*% samp[i,])

newdat$lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdat$upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatdoy <- as.data.frame(newdat[newdat$migration =="0" & newdat$year=="2017" & newdat$experience=="0",])
datadoy <- as.data.frame(data)
names(datadoy)[names(datadoy) == 'nr_fledglings'] <- 'fit'
names(datadoy)[names(datadoy) == 'day_incubate'] <- 'day_incubate_unz'

fig4 <- ggplot(newdatdoy, aes(x = day_incubate_unz, y = fit)) + geom_jitter(data = datadoy, aes(x = day_incubate_unz, y = fit, group = migration, colour = migration), size = .75) + scale_colour_manual(values = c("Orange", "Sky Blue"), labels = c("Resident", "Migrant")) + geom_line(stat = "identity", aes()) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .5) + labs(x = "Breeding Start (Julian date)", y = "Fledgling Number", colour = "Migration\nStrategy") + theme_classic() + theme(text = element_text(size = 10), legend.title=element_blank(), plot.margin = unit(c(0.25,0.25,.25,0.25), "cm"), legend.position = c(.85, .9), legend.spacing.x = unit(.1, 'cm'), legend.key.size = unit(.3, "cm")) + coord_cartesian(clip = "off") + scale_x_continuous(breaks=c(80,90,100,110,120,130))

ggsave(filename = "Figure_4.jpeg", plot = fig4, path = path, width = 8.5, height = 6, units = "cm", device='jpeg', dpi=600)