#### R Script to produce the analyses and graphs for the following paper:

# Early-life begging effort reduces adult body mass but strengthens behavioural 
# defence of the rate of energy intake in European starlings (Sturnus vulgaris)

# Authors: Jonathon Dunn1, Clare Andrews1, Daniel Nettle1 and Melissa Bateson1

# Script written by Jonathon Dunn
# 13 October 2017: first version of script for submission

# Notes on the variable names:
# 1. The developmental treatment called Effort in the paper is referred to as 'Work' in this script
# 2. When the data frames have been manipulated to give summaries, the naming convention is: 
# <name data frame data based on> . <response variable(s)> . <variable(s) data split by>

packages<-c("ggplot2",
            "plyr",
            "dplyr",
            "reshape2",
            "lme4",
            "broom",
            "tidyr",
            "gridExtra",
            "lubridate",
            "MuMIn",
            "arm",
            "scales")

lapply(packages, require, character.only = TRUE)

# Set the ggplot2 theme to be theme_bw() globally and specify accessory functions
# needed for analysis and plotting
theme_set(theme_bw())

# Create length function that can handle NAs
length2 <- function (x, na.rm=FALSE) {
  if (na.rm) sum(!is.na(x))
  else       length(x)
}

# Create function to adjust number of decimal places in ggplot axes
scaleFUN <- function(x) sprintf("%.2f", x)

#### Read in the data ####

# Setwd
setwd("X:/batesonlab/Early life stress/2014 birds/CFL operant/JD1_manuscript/Scripts and data")

# Read in .csv files
a<-read.csv("weights.csv")
b<-read.csv("latencies.csv")
d<-read.csv("individual.csv")

# Set up the colour scheme for any graphs that show Amount
a$Amount<-factor(a$Amount,levels=c("Lean","Plenty"))
b$Amount<-factor(b$Amount,levels=c("Lean","Plenty"))
d$Amount<-factor(d$Amount,levels=c("Lean","Plenty"))

# Set up the colour scheme for any graphs that show Work
a$Work<-factor(a$Work,levels=c("Hard","Easy"))
b$Work<-factor(b$Work,levels=c("Hard","Easy"))
d$Work<-factor(d$Work,levels=c("Hard","Easy"))

################ Weights ##############

# Use dataframe a

# Get dates in the correct format 
a$date_weighed=as.Date(dmy(a$date_weighed))

# Calculate the study duration and min/max dates
max(a$date_weighed)-min(a$date_weighed)
max(a$date_weighed)
min(a$date_weighed)

# Calculate the approximate number of times each bird was weighed for 'general' reasons
a%>%
  filter(reason_weighed=="general")%>%
  group_by(Bird)%>%
  summarise(number_general=n())%>%
  summarise(mean=mean(number_general),sd=sd(number_general))

# Calculate the approximate number of times each bird was caught
# after D56 and for reasons other than 'general'
a%>%
  filter(reason_weighed!="general"&reason_weighed!="D56")%>%
  group_by(Bird)%>%
  summarise(number_general=n())%>%
  summarise(mean=mean(number_general),sd=sd(number_general))

# Exclude PURWHI bird on 02-2015 as was ill, giving uncharacteristically low weight
a1<-subset(a,weight!=59.4)

# Count of how many times each bird was weighed
a1%>%
    group_by(Bird) %>%
    summarise(count=n())%>%
    summarise(mean=mean(count),sd=sd(count))

# Calculate the total number of mass measurements after exclusion
nrow(a1)

# Model 1: Predictors of body mass
m1<-lmer(weight~sex+tarsusD56+time.elapsed*Amount*Work*come_from+(1|NatalNest/Bird),
                    data=a1,
                    REML=F,
                    na.action=na.fail)
m1.s<-standardize(m1,standardize.y=F)
summary(m1)
r.squaredGLMM(m1)
m1.d<-dredge(m1.s)
(m1.ds<-subset(m1.d, delta < 2))
nrow(m1.d)
nrow(m1.ds)
m1.a<-model.avg(m1.ds)
summary(m1.a)
confint(m1.a)

# Create Figure 1A
(a.bc.housing<-a1%>%
    group_by(come_from,Bird)%>%
    summarise(mean.bc=mean(body.condition,na.rm=T))%>%
    group_by(come_from)%>%
    summarise(mean=mean(mean.bc,na.rm=T), sd=sd(mean.bc,na.rm=T),se = sd(mean.bc,na.rm=T)/sqrt(length2(mean.bc,na.rm=T))))

levels(a.bc.housing$come_from)
levels(a.bc.housing$come_from)<-c("Aviary","Cage")

(g.bc.housing<-ggplot(a.bc.housing,aes(x=come_from,y=mean))+
    geom_bar(stat="identity")+
    geom_errorbar(aes(ymin=mean-se,ymax=mean+se),position=position_dodge(0.1),width=0.1,linetype="solid")+
    xlab("Housing")+
    ylab("Residual mass")+
    ggtitle("A"))

# Create Figure 1B
(a.bc.sex<-a1%>%
    group_by(sex,Bird)%>%
    summarise(mean.bc=mean(body.condition,na.rm=T))%>%
    group_by(sex)%>%
    summarise(mean=mean(mean.bc,na.rm=T), sd=sd(mean.bc,na.rm=T),se = sd(mean.bc,na.rm=T)/sqrt(length2(mean.bc,na.rm=T))))

levels(a.bc.sex$sex)<-c("Female","Male")

(g.bc.sex<-ggplot(a.bc.sex,aes(x=sex,y=mean))+
    geom_bar(stat="identity")+
    geom_errorbar(aes(ymin=mean-se,ymax=mean+se),position=position_dodge(0.1),width=0.1,linetype="solid")+
    xlab("Sex")+
    ylab("Residual mass")+
    ggtitle("B"))

# Create Figure 1C
(a.mass.tarsus<-a1%>%
    group_by(tarsusD56)%>%
    summarise(mean=mean(weight,na.rm=T),se = sd(weight,na.rm=T)/sqrt(length2(weight,na.rm=T))))

(g.mass.tarsus<-ggplot(a.mass.tarsus,aes(x=tarsusD56,y=mean))+
    geom_point()+
    xlab("Tarsus length (mm)")+
    ylab("Mass (g)")+
    stat_smooth(method="lm",se=F)+
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se),position=position_dodge(0.1),width=0.1,linetype="solid")+ 
    ggtitle("C"))

# Create Figure 1D
(g.age<-ggplot(a1,aes(x=time.elapsed+56,y=body.condition))+
    geom_point()+
    xlab("Age (days)")+
    ylab("Residual mass")+
    stat_smooth(method="lm")+
    ggtitle("D"))

# Create Figure 1E
housing_names <- c(
  flight = "Aviary",
  individual = "Cage"
)

(g.age.housing<-ggplot(a1,aes(x=time.elapsed+56,y=body.condition))+
    geom_point()+
    xlab("Age (days)")+
    ylab("Residual mass")+
    stat_smooth(method="lm")+
    facet_wrap(~come_from,labeller=labeller(come_from=housing_names),ncol=2)+
    ggtitle("E"))

### Draw Figure 1 for the paper ###
Fig1<-grid.arrange(arrangeGrob(g.bc.housing,g.bc.sex, ncol=2),
                    arrangeGrob(g.mass.tarsus, g.age, ncol=2),
                    g.age.housing, 
                    ncol=1)

ggsave("Fig1.jpeg",Fig1,dpi=300,height=10,width=6)

# Create Figure 2A
(a.bc.work<-a1%>%
    group_by(Bird,Work)%>%
    summarise(mean.bc=mean(body.condition,na.rm=T))%>%
    group_by(Work)%>%
    summarise(mean=mean(mean.bc),se=sd(mean.bc,na.rm=T)/sqrt(length2(mean.bc,na.rm=T))))

(g.bc.work<-ggplot(a.bc.work,aes(x=Work,y=mean,fill=Work))+
    geom_bar(stat="identity")+
    geom_errorbar(aes(ymin=mean-se,ymax=mean+se),position=position_dodge(0.1),width=0.1,linetype="solid")+
    xlab("Effort")+
    ylab("Residual mass")+
    scale_fill_manual(values=c("gray40","gray80"),guide=F)+ 
    theme(axis.text.x = element_text(angle = 0, hjust = 1))+
    ggtitle("A"))

# Create Figure 2B
(a.bc.work.housing<-a1%>%
  group_by(Work,come_from,Bird)%>%
  summarise(mean.bc=mean(body.condition,na.rm=T))%>%
  group_by(Work,come_from)%>%
  summarise(mean=mean(mean.bc,na.rm=T), se = sd(mean.bc,na.rm=T)/sqrt(length2(mean.bc,na.rm=T))))

(g.bc.work.housing<-ggplot(a.bc.work.housing,aes(x=Work,y=mean,fill=Work))+
    geom_bar(stat="identity")+
    geom_errorbar(aes(ymin=mean-se,ymax=mean+se),position=position_dodge(0.1),width=0.1,linetype="solid")+
    facet_wrap(~come_from,labeller=labeller(come_from=housing_names))+
    xlab("Effort")+
    ylab("Residual mass")+
    scale_fill_manual(values=c("gray40","gray80"),guide=F)+ 
    theme(axis.text.x = element_text(angle = 0, hjust = 1))+
    ggtitle("B"))

# Create Figure 2C
(a.bc.work.housing.age<-a1%>%
  group_by(come_from,Work,bin.4m,Bird)%>%
  summarise(mean.bc=mean(body.condition,na.rm=T))%>%
  group_by(come_from,Work,bin.4m)%>%
  summarise(mean=mean(mean.bc,na.rm=T), se = sd(mean.bc,na.rm=T)/sqrt(length2(mean.bc,na.rm=T))))

(g.bc.work.housing.age<-ggplot(a.bc.work.housing.age,aes(x=bin.4m,y=mean,col=Work))+
    geom_point(position=position_dodge(0.1))+
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se), position=position_dodge(0.1),width=0.1,linetype="solid") +
    facet_wrap(~come_from,labeller=labeller(come_from=housing_names),ncol=2)+
    xlab("Date")+ ylab("Residual mass")+
    geom_line(aes(col=Work,group=Work),position=position_dodge(0.1))+ 
    theme(axis.text.x = element_text(angle = 90, hjust = 1,size=8))+
    ggtitle("C")+
    scale_color_manual(values=c("gray40","gray80"),labels=c("Hard","Easy"),name="Effort")+
    theme(strip.text.x = element_blank(),
          strip.background = element_rect(colour="white", fill="white"),
          legend.position=c(.9,.75)))

### Draw Figure 2 for the paper ###
Fig2<-grid.arrange(arrangeGrob(g.bc.work,g.bc.work.housing,ncol=2),g.bc.work.housing.age,ncol=1)

ggsave("Fig2.jpeg",Fig2,dpi=300,height=6,width=6)

##########Latencies#######

# Use dataframe b

# Model 2: predictors of trial latency in the operant experiment
m2<-lmer(log(Latency)~Ratio*Work*Amount+sex+tarsusD56+(1|NatalNest/Bird),data=b,REML=F,na.action=na.fail)
summary(m2)
r.squaredGLMM(m2)
m2.s<-standardize(m2,standardize.y=F)
m2.d <- dredge(m2.s) 
nrow(m2.d) 
(m2.ds<-subset(m2.d, delta < 2))
nrow(m2.ds) 
m2.a<-model.avg(m2.ds)
summary(m2.a)  
confint(m2.a)  

# Plot Figure 3A

# Calculate mean trial latency for each ratio value with between-bird SEs
(b.latency.ratio.amount<-b%>%
    group_by(Ratio,Bird,Amount)%>%
    summarise(mean1=mean(log(Latency),na.rm=T))%>%
    group_by(Ratio,Amount)%>%
    summarise(mean=mean(mean1,na.rm=T), se = sd(mean1,na.rm=T)/sqrt(length2(mean1,na.rm=T))))

(g.latency.ratio.amount<-ggplot(b.latency.ratio.amount,aes(x=Ratio,y=mean,col=Amount))+
    geom_point()+
    stat_smooth(method="lm",se=F)+
    facet_wrap(~Amount)+
    scale_color_manual(values=c("hotpink2","royalblue2"),guide=F)+
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se),position=position_dodge(0.1),width=0.1,linetype="solid")+
    ylab("Trial latency log(s)")+
    xlab("Ratio requirement")+
    scale_x_continuous(breaks=pretty_breaks())+
    ggtitle("A")+
    scale_y_continuous(labels=scaleFUN,limits=c(1.5,3.75)))

# Plot Figure 3B

# Calculate mean trial latency for each ratio value with between-bird SEs
(b.latency.ratio.work<-b%>%
    group_by(Ratio,Bird,Work)%>%
    summarise(mean1=mean(log(Latency),na.rm=T))%>%
    group_by(Ratio,Work)%>%
    summarise(mean=mean(mean1,na.rm=T), se = sd(mean1,na.rm=T)/sqrt(length2(mean1,na.rm=T))))

b.latency.ratio.work$Work<-factor(b.latency.ratio.work$Work,levels=c("Hard","Easy"))

(g.latency.ratio.work<-ggplot(b.latency.ratio.work,aes(x=Ratio,y=mean,col=Work))+
    geom_point()+
    stat_smooth(method="lm",se=F)+
    facet_wrap(~Work)+
    scale_color_manual(values=c("gray40","gray80"),guide=F)+
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se),position=position_dodge(0.1),width=0.1,linetype="solid")+
    ylab("Trial latency log(s)")+
    xlab("Ratio requirement")+
    scale_x_continuous(breaks=pretty_breaks())+
    ggtitle("B")+
    scale_y_continuous(labels=scaleFUN,limits=c(1.5,3.75)))

# Plot Figure 3C

# Use dataframe d

# Calculate mean trial latency with between-bird SEs
(d.latency.amount<-d%>%
  group_by(Amount)%>%
  summarise(mean=mean(latency.mean,na.rm=T), se = sd(latency.mean,na.rm=T)/sqrt(length2(latency.mean,na.rm=T))))

(g.latency.amount<-ggplot(d.latency.amount, aes(x=Amount, y=mean, fill=Amount)) + 
    geom_bar(stat="identity",position=position_dodge(0.1)) +
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1, position=position_dodge(0.1),linetype="solid") +
    xlab("Amount")+ 
    ylab("Trial latency log(s)")+
    ggtitle("C")+
    theme(plot.title = element_text(hjust = 0))+ 
    scale_fill_manual(values=c("hotpink2","royalblue2"),guide=F)+
    scale_y_continuous(labels=scaleFUN,limits=c(0,3.2)))

# Plot Figure 3D

# Calculate mean trial latency with between-bird SEs
(d.latency.work<-d%>%
  group_by(Work)%>%
  summarise(mean=mean(latency.mean,na.rm=T), se = sd(latency.mean,na.rm=T)/sqrt(length2(latency.mean,na.rm=T))))

d.latency.work$Work<-factor(d.latency.work$Work,levels=c("Hard","Easy"))

(g.latency.work<-ggplot(d.latency.work, aes(x=Work, y=mean, fill=Work)) + 
    geom_bar(stat="identity",position=position_dodge(0.1)) +
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1, position=position_dodge(0.1),linetype="solid") +
    xlab("Effort")+ 
    ylab("Trial latency log(s)")+
    ggtitle("D")+
    theme(plot.title = element_text(hjust = 0))+ 
    scale_fill_manual(values=c("gray40","gray80"),guide=F)+
    scale_y_continuous(labels=scaleFUN,limits=c(0,3.2)))

# Plot Figure 3E

# Calculate mean latency slopes with between-bird SEs from bird-level summary data frame
(d.slope.amount<-d%>%
    group_by(Amount)%>%
    summarise(mean=mean(latency.slope.mean,na.rm=T), se = sd(latency.slope.mean,na.rm=T)/sqrt(length2(latency.slope.mean,na.rm=T))))

(g.slope.amount<-ggplot(d.slope.amount, aes(x=Amount, y=mean, fill=Amount)) + 
    geom_bar(stat="identity",position=position_dodge(0.1)) +
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1, position=position_dodge(0.1),linetype="solid") +
    xlab("Amount")+ 
    ylab("Latency slope")+
    ggtitle("E")+
    theme(plot.title = element_text(hjust = 0))+ 
    scale_fill_manual(values=c("hotpink2","royalblue2"),guide=F)+
    scale_y_reverse(lim=c(0.12,0),labels=scaleFUN))

# Plot Figure 3F

# Calculate mean latency slopes with between-bird SEs from bird-level summary data frame
(d.slope.work<-d%>%
    group_by(Work)%>%
    summarise(mean=mean(latency.slope.mean,na.rm=T), se = sd(latency.slope.mean,na.rm=T)/sqrt(length2(latency.slope.mean,na.rm=T))))

d.slope.work$Work<-factor(d.slope.work$Work,levels=c("Hard","Easy"))

(g.slope.work<-ggplot(d.slope.work, aes(x=Work, y=mean, fill=Work)) + 
    geom_bar(stat="identity",position=position_dodge(0.1)) +
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1, position=position_dodge(0.1),linetype="solid") +
    xlab("Effort")+
    ylab("Latency slope")+
    ggtitle("F")+
    theme(plot.title = element_text(hjust = 0))+
    scale_fill_manual(values=c("gray40","gray80"),guide=F)+
    scale_y_reverse(lim=c(0.12,0),labels=scaleFUN))

### Draw Figure 3 for paper ###
Fig3<-grid.arrange(g.latency.ratio.amount,
                   g.latency.ratio.work,
                   g.latency.amount,
                   g.latency.work,
                   g.slope.amount,
                   g.slope.work,
                   ncol=2)


ggsave("Fig3.jpeg",Fig3,dpi=300,width=6,height=8)

########### Relationship between latency slope and body mass

# Use dataframe d

# Model 3: relationship between latency slope and body mass
m3<-lmer(latency.slope.mean~mean.wt*Work+tarsusD56+(1|NatalNest),data=d)
m3.s<-standardize(m3,standardize.y=F)
summary(m3.s)
confint(m3.s)
m3a.s<-update(m3.s,.~.-mean.wt:Work)
summary(m3a.s)
confint(m3a.s)  

### Draw Figure 4 for the paper ###
(Fig4<-ggplot(d,aes(x=mean.bc,y=latency.slope.mean,col=Work,group=Work))+
    geom_point()+
    stat_smooth(method="lm",col="black",se=F)+
    xlab("Residual mass")+
    ylab("Latency slope")+
    ggtitle("")+
    scale_color_manual(values=c("gray40","gray80"),labels=c("Hard","Easy"),name="Effort")+
    scale_y_reverse()+ 
    theme(strip.text.x = element_blank(),
          strip.background = element_rect(colour="white", fill="white"),
          legend.position=c(.8,.25)))

ggsave("Fig4.jpeg",Fig4,dpi=300,height=3,width=3)

################# Online Supporting Material ############

# Comparing Staddon parameters to latency parameters

# Create subset from dataframe d where only -ve Staddon slopes included for comparison
d1<-subset(d,Staddon.slope<0)

# Create Figure S1A
(g.slopes<-ggplot(d1,aes(x=Staddon.slope,y=latency.slope.median))+
    geom_point()+
    stat_smooth(method="lm")+
    xlab("Staddon slope")+
    ylab("Defence of average trial latency log(s)")+
    ggtitle("A"))

# Create Figure S1B
(g.intercept.latency<-ggplot(d1,aes(x=Staddon.y.intercept,y=latency.median))+
    geom_point()+
    stat_smooth(method="lm")+
    xlab("Staddon y.intercept")+
    ylab("Average trial latency log(s)")+
    ggtitle("B")+ 
    scale_y_continuous(labels = function(x) format(x, nsmall = 2)))

### Draw Figure S1 for the online supporting materials ###
FigS1<-grid.arrange(g.slopes,g.intercept.latency,ncol=2)

ggsave("FigS1.jpeg",FigS1,dpi=300,width=6,height=3)

# Calculate correlation co-efficient for relationships shown in Figure S1
cor.test(d1$Staddon.slope,d1$latency.slope.median,alternative="g",method="pearson") #greater corresponds to +ve slope
cor.test(d1$Staddon.y.intercept,d1$latency.median,alternative="l",method="pearson") #greater corresponds to -ve slope