library(ggplot2) 
library(lme4)
library(lmerTest)
library(brms)
library(ppcor)
library(grid)
library(RColorBrewer)
library(jcolors)

######################## ANALYSES ################################
# Analysis 1 #
dat1 <- read.table("Analysis 1.txt", header=T)

# correlations #
cor.test(dat1$On.bout_duration,dat1$Off.bout_duration)
cor.test(dat1$On.bout_duration,dat1$Off.bout_number)
cor.test(dat1$Off.bout_duration,dat1$Off.bout_number)

# Analysis 2 on #
dat2on <- read.table("Analysis 2 on.txt", header=T)
twoon<-lmer(Model_egg_temperature~On.bout_duration+Ambient_temperature+Wind_speed+(On.bout_duration|NestID), data=dat2on)
summary(twoon)
anova(twoon)

# interaction #
twoon<-lmer(Model_egg_temperature~On.bout_duration*Ambient_temperature+Wind_speed+(On.bout_duration|NestID), data=dat2on)
summary(twoon)
anova(twoon)

# Analysis 2 off #
dat2off <- read.table("Analysis 2 off.txt", header=T)
twooff<-lmer(Model_egg_temperature~Off.bout_duration+Ambient_temperature+(Off.bout_duration|NestID), data=dat2off)
summary(twooff)
anova(twooff)

# interaction #
twooff<-lmer(Model_egg_temperature~Off.bout_duration*Ambient_temperature+(Off.bout_duration|NestID), data=dat2off)
summary(twooff)
anova(twooff)

## Analysis 3 on##
dat3on <- read.table("Analysis 3 on.txt", header=T)

threeon3<-lmer(log(On.bout_duration)~poly(Centered_ambient_temperature,3)+Centered_wind_speed+(Centered_ambient_temperature|NestID), data=dat3on)
threeon2<-lmer(log(On.bout_duration)~poly(Centered_ambient_temperature,2)+Centered_wind_speed+(Centered_ambient_temperature|NestID), data=dat3on)
threeon1<-lmer(log(On.bout_duration)~poly(Centered_ambient_temperature,1)+Centered_wind_speed+(Centered_ambient_temperature|NestID), data=dat3on)

# test of poly2 and poly3 effects
anova(threeon1,threeon3)
anova(threeon1,threeon2)
summary(threeon3)
summary(threeon2)

# no poly effects
summary(threeon1)
anova(threeon1)

# test of random effects
threeonfull<-lmer(log(On.bout_duration)~Centered_ambient_temperature+Centered_wind_speed+(Centered_ambient_temperature|NestID), data=dat3on)
threeonnocorr<-lmer(log(On.bout_duration)~Centered_ambient_temperature+Centered_wind_speed+(1|NestID)+(0+Centered_ambient_temperature|NestID), data=dat3on)
threeonnoslope<-lmer(log(On.bout_duration)~Centered_ambient_temperature+Centered_wind_speed+(1|NestID), data=dat3on)
threeonnone<-lm(log(On.bout_duration)~Centered_ambient_temperature+Centered_wind_speed, data=dat3on)

# tests of random slope/int correlation, random slope, and random intercept
anova(threeonfull,threeonnocorr)
anova(threeonfull,threeonnoslope)
anova(threeonnoslope,threeonnonone)


## Analysis 3 off ##
dat3off <- read.table("Analysis 3 off.txt", header=T)

# test of random effects
threeofffull<-lmer(log(Off.bout_duration)~Centered_ambient_temperature+(Centered_ambient_temperature|NestID), data=dat3off)
threeoffnocorr<-lmer(log(Off.bout_duration)~Centered_ambient_temperature+(1|NestID)+(0+Centered_ambient_temperature|NestID), data=dat3off)
threeoffnoslope<-lmer(log(Off.bout_duration)~Centered_ambient_temperature+(1|NestID), data=dat3off)
threeoffnone<-lm(log(Off.bout_duration)~Centered_ambient_temperature, data=dat3off)

# tests of random slope/int correlation, random slope, and random intercept
anova(threeofffull,threeoffnocorr)
anova(threeofffull,threeoffnoslope)
anova(threeoffnoslope,threeoffnonone)

# test of wind speed #
threeoff<-lmer(log(Off.bout_duration)~Centered_ambient_temperature+Centered_wind_speed+(1|NestID)+(0+Centered_ambient_temperature|NestID), data=dat3off)
summary(threeoff)
anova(threeoff)


# Analysis 4 #
dat4 <- read.table("Analysis 4.txt", header=T)

coevegg <- bf(Coev_model_egg_temperature|mi() ~ Coev_ambient_temperature + (1|p|NestID))
onduration <- bf(Log_on.bout_duration|mi() ~ Centered_ambient_temperature + Centered_wind_speed + (ambic|p|NestID))
offduration <- bf(Log_off.bout_duration|mi() ~ Centered_ambient_temperature + (1|p|NestID)+(0+ambic|p|NestID))
offnumber <- bf(Off.bout_number|mi() ~ (1|p|NestID))
attentive <- bf(Attentiveness|mi() ~ (1|p|NestID))
eggtemp <- bf(Model_egg_temperature|mi() ~ Ambient_temperature + (1|p|NestID))
four <- brm(coevegg+onduration+offduration+offnumber+attentive+eggtemp, data = dat4)
summary(four)

######################## FIGURES ################################

# Figure 1 #
fig1 <- read.table("Figure 1.txt", header=T)

p2<-ggplot(data=fig1, aes(x=Julian, y=Mean_temp))+labs(x="Breeding initiation date (from 1 June)", y="Mean ambient temp (°C)")+geom_point(size=5)+theme_bw()+geom_smooth(method=lm, color="black", se=FALSE, size=1, linetype=2)+theme(legend.position="none")+theme(axis.line = element_line(colour = 'black', size = 1))+theme(axis.title.y=element_text(margin=margin(0,15,0,0)))+theme(axis.title.x=element_text(margin=margin(15,0,0,0)))+theme(axis.title = element_text(size =28)) +theme(axis.text.x=element_text(angle=0, size=25, vjust=0.5, color="black"))+theme(axis.text.y=element_text(angle=0, size=25, vjust=0.5, color="black"))
p3<-p2+scale_x_continuous(breaks=c(0,25,50,75), limits=c(0,75))
p4<-p3+scale_y_continuous(breaks=c(5,10,15,20), limits=c(5,20))+theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line.x = element_line(colour = 'black', size = 1),axis.line.y = element_line(colour = 'black', size = 1))
p5<-p4 + theme(axis.ticks.length = unit(0.5, "cm"))
p5+geom_pointrange(aes(ymin=Lower, ymax=Upper))

# Climate Analysis #
climate<-lm(Mean_temp~Julian, data=fig1)
summary(climate)

# Figure 2 #
# Figure 2 was created online and made public in plotly chart studio #
# Figure 2A in plotly:https://plotly.com/~agco232/8/
# Figure 2A data in plotly: https://plotly.com/~agco232/7/
# Figure 2B in plotly:https://plotly.com/~agco232/10/
# Figure 2B data in plotly: https://plotly.com/~agco232/9/


# Figure 3 #
fig3ab <- read.table("Figure 3AB.txt", header=T)
fig3cd <- read.table("Figure 3CD.txt", header=T)

# custom color palette #
nb.cols <- 23
mycolors <- colorRampPalette(jcolors(palette="rainbow"))(23)

# Figure 3A #
p1<-ggplot(data=fig3ab, aes(x=Ambient_temperature, y=On.bout_duration))+geom_smooth(method="lm", se=F)+aes(group=NestID, color=NestID)+scale_color_manual(values=mycolors)
p2<-p1+labs(x="Centred ambient temp (°C)", y="On-bout duration (min)")+theme_bw()+theme(legend.position="none")+theme(axis.line = element_line(colour = 'black', size = 1))+theme(axis.title.y=element_text(margin=margin(0,15,0,0)))+theme(axis.title.x=element_text(margin=margin(15,0,0,0)))+theme(axis.title = element_text(size = 32)) +theme(axis.text.x=element_text(angle=0, size=25, vjust=0.5, color="black"))+theme(axis.text.y=element_text(angle=0, size=25, vjust=0.5, color="black"))
p3<-p2+scale_x_continuous(breaks=c(-20,-10,0,10,20), limits=c(-20,20))
p4<-p3+scale_y_continuous(breaks=c(0,15,30,45,60), limits=c(0,60))+theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line.x = element_line(colour = 'black', size = 1),axis.line.y = element_line(colour = 'black', size = 1))
p4+ theme(axis.ticks.length = unit(0.5, "cm"))+geom_vline(xintercept = 0, linetype="dashed", size=0.5)

# Figure 3B #
p1<-ggplot(data=fig3ab, aes(x=Ambient_temperature, y=Off.bout_duration))+geom_smooth(method="lm", se=F)+aes(group=NestID, color=NestID)+scale_color_manual(values=mycolors)
p2<-p1+labs(x="Centred ambient temp (°C)", y="Off-bout duration (min)")+theme_bw()+theme(legend.position="none")+theme(axis.line = element_line(colour = 'black', size = 1))+theme(axis.title.y=element_text(margin=margin(0,15,0,0)))+theme(axis.title.x=element_text(margin=margin(15,0,0,0)))+theme(axis.title = element_text(size = 32)) +theme(axis.text.x=element_text(angle=0, size=25, vjust=0.5, color="black"))+theme(axis.text.y=element_text(angle=0, size=25, vjust=0.5, color="black"))
p3<-p2+scale_x_continuous(breaks=c(-20,-10,0,10,20), limits=c(-20,20))
p4<-p3+scale_y_continuous(breaks=c(0,7,14,21,28), limits=c(0,28))+theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line.x = element_line(colour = 'black', size = 1),axis.line.y = element_line(colour = 'black', size = 1))
p4 + theme(axis.ticks.length = unit(0.5, "cm"))+ geom_vline(xintercept = 0, linetype="dashed", size=0.5)

# Figure 3C #
p2<-ggplot(data=fig3cd, aes(x=On.bout_intercept, y=On.bout_slope))+geom_point(size=5)+labs(x="On-bout intercept", y="On-bout slope")+theme_bw()+theme(legend.position="none")+theme(axis.line = element_line(colour = 'black', size = 1))+theme(axis.title.y=element_text(margin=margin(0,15,0,0)))+theme(axis.title.x=element_text(margin=margin(15,0,0,0)))+theme(axis.title = element_text(size = 32)) +theme(axis.text.x=element_text(angle=0, size=25, vjust=0.5, color="black"))+theme(axis.text.y=element_text(angle=0, size=25, vjust=0.5, color="black"))
p3<-p2+scale_x_continuous(breaks=c(10,15,20,25,30,35), limits=c(10,36))
p4<-p3+scale_y_continuous(breaks=c(-1,0,1,2,3), limits=c(-1,3))+theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line.x = element_line(colour = 'black', size = 1),axis.line.y = element_line(colour = 'black', size = 1))
p4 + theme(axis.ticks.length = unit(0.5, "cm"))+aes(color=NestID)+scale_color_manual(values=mycolors)

# Figure 3D #
p2<-ggplot(data=fig3cd, aes(x=Off.bout_intercept, y=Off.bout_slope))+geom_point(size=5)+labs(x="Off-bout intercept", y="Off-bout slope")+theme_bw()+theme(legend.position="none")+theme(axis.line = element_line(colour = 'black', size = 1))+theme(axis.title.y=element_text(margin=margin(0,15,0,0)))+theme(axis.title.x=element_text(margin=margin(15,0,0,0)))+theme(axis.title = element_text(size = 32)) +theme(axis.text.x=element_text(angle=0, size=25, vjust=0.5, color="black"))+theme(axis.text.y=element_text(angle=0, size=25, vjust=0.5, color="black"))
p3<-p2+scale_x_continuous(breaks=c(10,12,14,16,18,20), limits=c(10,20))
p4<-p3+scale_y_continuous(breaks=c(0,0.5,1,1.5,2), limits=c(0,2))+theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line.x = element_line(colour = 'black', size = 1),axis.line.y = element_line(colour = 'black', size = 1))
p4 + theme(axis.ticks.length = unit(0.5, "cm"))+aes(color=NestID)+scale_color_manual(values=mycolors)

# FIGURE 4 #
fig4 <- read.table("Figure 4.txt", header=T)

# Figure 4A #
p2<-ggplot(data=fig4, aes(x=On.bout_slope, y=Attentiveness))+geom_smooth(method="lm", se=FALSE, color="black")+geom_point(size=5)+labs(x="", y="Attentiveness intercept")+theme_bw()+theme(legend.position="none")+theme(axis.line = element_line(colour = 'black', size = 1))+theme(axis.title.y=element_text(margin=margin(0,15,0,0)))+theme(axis.title = element_text(size = 32)) +theme(axis.text.x=element_text(angle=0, size=25, vjust=0.5, color="black"))+theme(axis.text.y=element_text(angle=0, size=25, vjust=0.5, color="black"))
p3<-p2+scale_y_continuous(breaks=c(0.55,0.60,0.65,0.70,0.75), limits=c(0.55,0.75))
p4<-p3+scale_x_continuous(breaks=c(-1,0,1,2,3), limits=c(-1,3))+theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line.x = element_line(colour = 'black', size = 1),axis.line.y = element_line(colour = 'black', size = 1))
p4 + theme(axis.ticks.length = unit(0.5, "cm"))

# Figure 4B #
p2<-ggplot(data=fig4, aes(x=On.bout_slope, y=Model_egg_temperature))+geom_smooth(method="lm", se=FALSE, color="black")+geom_point(size=5)+labs(x="", y="Egg temp intercept (°C)")+theme_bw()+theme(legend.position="none")+theme(axis.line = element_line(colour = 'black', size = 1))+theme(axis.title.y=element_text(margin=margin(0,15,0,0)))+theme(axis.title = element_text(size = 32)) +theme(axis.text.x=element_text(angle=0, size=25, vjust=0.5, color="black"))+theme(axis.text.y=element_text(angle=0, size=25, vjust=0.5, color="black"))
p3<-p2+scale_y_continuous(breaks=c(24,26,28,30,32), limits=c(24,32))
p4<-p3+scale_x_continuous(breaks=c(-1,0,1,2,3), limits=c(-1,3))+theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line.x = element_line(colour = 'black', size = 1),axis.line.y = element_line(colour = 'black', size = 1))
p4 + theme(axis.ticks.length = unit(0.5, "cm"))

# Figure 4C #
p2<-ggplot(data=fig4, aes(x=On.bout_slope, y=Coev_egg_temperature))+geom_smooth(method="lm", se=FALSE, color="black")+geom_point(size=5)+labs(x="On-bout slope", y="CV egg temp intercept")+theme_bw()+theme(legend.position="none")+theme(axis.line = element_line(colour = 'black', size = 1))+theme(axis.title.y=element_text(margin=margin(0,15,0,0)))+theme(axis.title.x=element_text(margin=margin(15,0,0,0)))+theme(axis.title = element_text(size = 32)) +theme(axis.text.x=element_text(angle=0, size=25, vjust=0.5, color="black"))+theme(axis.text.y=element_text(angle=0, size=25, vjust=0.5, color="black"))
p3<-p2+scale_y_continuous(breaks=c(0.02,0.04,0.06,0.08,0.1), limits=c(0.02,0.1))
p4<-p3+scale_x_continuous(breaks=c(-1,0,1,2,3), limits=c(-1,3))+theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line.x = element_line(colour = 'black', size = 1),axis.line.y = element_line(colour = 'black', size = 1))
p4 + theme(axis.ticks.length = unit(0.5, "cm"))