# A marker of biological age explains individual variation in the strength of the adult stress response # Clare Andrews, Daniel Nettle, Maria Larriva, Robert Gillespie, Sophie Reichert, Ben O Brilot, Thomas Bedford, Pat Monaghan, Karen A Spencer, Melissa Bateson* # R script accompanying the above paper. # Script written by Clare Andrews, Daniel Nettle and Melissa Bateson # Version for publication prepared by Melissa Bateson on 22nd August 2017 # This script produces the statistical analyses and figures in the paper library(AICcmodavg) library(gridExtra) library(dplyr) library(ggplot2) library(lme4) library(reshape2) library(psych) library(metafor) library(grid) # Read in the data for the two cohorts d12=read.csv("Data_file_1_Andrews_2012_cohort_CORT_data.csv") d13=read.csv("Data_file_2_Andrews_2013_cohort_CORT_data.csv") ##### Basic descriptive stats on CORT variables for Table 1 ##### # reverse the calculation of deltaCORT to make the direction the same as deltaTL (i.e.timepoint2-timepoint1) d12$deltacort=-d12$deltacort d13$DeltaCort=-d13$DeltaCort # 2012 describe(d12$Cort0mins) describe(d12$Cort15mins) describe(d12$Cort30mins) describe(d12$peakcort) describe(d12$deltacort) d12$increase=NA d12$decrease=NA d12$increase[d12$deltacort<=0] = 1 d12$decrease[d12$deltacort>0] = 1 inc12 = sum(d12$increase,na.rm = TRUE) dec12 = sum(d12$decrease,na.rm = TRUE) prop12 = dec12/(inc12+dec12) # 2013 describe(d13$Cort0) describe(d13$Cort15) describe(d13$Cort30) describe(d13$PeakCort) describe(d13$DeltaCort) d13$increase=NA d13$decrease=NA d13$increase[d13$DeltaCort<=0] = 1 d13$decrease[d13$DeltaCort>0] = 1 inc13 = sum(d13$increase,na.rm = TRUE) dec13 = sum(d13$decrease,na.rm = TRUE) prop13 = dec13/(inc13+dec13) # t tests t.test(d12$Cort0mins,d13$Cort0) t.test(d12$Cort15mins,d13$Cort15) t.test(d12$Cort30mins,d13$Cort30) t.test(d12$peakcort,d13$PeakCort) t.test(d12$deltacort,d13$DeltaCort) pat = data.frame(cohort=c("2012","2012","2013","2013"), change=c("dec","inc","dec","inc"),freq=c(dec12,inc12,dec13,inc13)) my.table = xtabs(freq~cohort+change, data=pat) ftable(my.table) summary(my.table) ##### Test for effects of potential nuisance variables ###### # 2012 # Body condition (table 1) m1=lmer(Cort0mins~body.condition + (1|Natal_nest), REML=FALSE, data=d12) summary(m1) drop1(m1, test="Chisq") m2=lmer(peakcort~body.condition + Cort0mins + (1|Natal_nest), REML=FALSE, data=d12) summary(m2) drop1(m2, test="Chisq") m3=lmer(deltacort~body.condition + Cort15mins + (1|Natal_nest), REML=FALSE, data=d12) summary(m3) drop1(m3, test="Chisq") # Age in days (table 1) m1=lmer(Cort0mins~Agedays + (1|Natal_nest), REML=FALSE, data=d12) summary(m1) drop1(m1, test="Chisq") m2=lmer(peakcort~Agedays + Cort0mins + (1|Natal_nest), REML=FALSE, data=d12) summary(m2) drop1(m2, test="Chisq") m3=lmer(deltacort~Agedays + Cort15mins + (1|Natal_nest), REML=FALSE, data=d12) summary(m3) drop1(m3, test="Chisq") # Days in cages (table 1) m1=lmer(Cort0mins~Days_in_cage + (1|Natal_nest), REML=FALSE, data=d12) summary(m1) drop1(m1, test="Chisq") m2=lmer(peakcort~Days_in_cage + Cort0mins + (1|Natal_nest), REML=FALSE, data=d12) summary(m2) drop1(m2, test="Chisq") m3=lmer(deltacort~Days_in_cage + Cort15mins + (1|Natal_nest), REML=FALSE, data=d12) summary(m3) drop1(m3, test="Chisq") # Time to get baseline sample (model S10) m1=lmer(Cort0mins~Delay_to_first_sample_sec + (1|Natal_nest), REML=FALSE, data=d12) summary(m1) drop1(m1, test="Chisq") # 2013 cohort of birds # Body condition (table 1) m1=lmer(Cort0~BodyCondition + (1|NatalNest), REML=FALSE, data=d13) summary(m1) drop1(m1, test="Chisq") m2=lmer(PeakCort~BodyCondition + Cort0 + (1|NatalNest), REML=FALSE, data=d13) summary(m2) drop1(m2, test="Chisq") m3=lmer(DeltaCort~BodyCondition + Cort15+ (1|NatalNest), REML=FALSE, data=d13) summary(m3) drop1(m3, test="Chisq") # Time to obtain sample (table 1) m1=lmer(Cort0~TimeToSampleS + (1|NatalNest), REML=FALSE, data=d13) summary(m1) drop1(m1, test="Chisq") #### Familial effects ##### # 2012 birds m1=lmer(Cort0mins~1 + (1|Natal_nest/Host_nest), REML=T, data=d12) summary(m1) # Natal nest component 23.17/(23.17+72.21) # Host nest component 0 m2=lmer(peakcort~1 + (1|Natal_nest/Host_nest), REML=T, data=d12) summary(m2) # Natal nest component for peak 219.1/(219.1+655.1) # Host nest 0 m3=lmer(deltacort~1 + (1|Natal_nest/Host_nest), REML=T, data=d12) summary(m3) # Natal nest 0 and host nest 0 col1=NULL col1[1:9]=rep("2012 cohort", times=9) col1[10:18]=rep("2013 cohort", times=9) col2=rep(c("Natal nest", "Host nest", "Residual"), times=6) col3=NULL col3[1:3]="Baseline" col3[4:6]="Peak" col3[7:9]="Delta" col3[10:12]="Baseline" col3[13:15]="Peak" col3[16:18]="Delta" fig1.data=data.frame(col1, col2, col3) View(fig1.data) fig1.data$value=NA fig1.data$value[1]=23.17/(23.17+72.21) fig1.data$value[2]=0/(23.17+72.21) fig1.data$value[3]=72.21/(23.17+72.21) fig1.data$value[4]=219.1/(219.1+655.1) fig1.data$value[5]=0/(219.1+655.1) fig1.data$value[6]=655.1/(219.1+655.1) fig1.data$value[7]=0 fig1.data$value[8]=0 fig1.data$value[9]=1 # 2013 birds # Base cort colnames(d13) m1=lmer(Cort0~1 + (1|NatalNest/Host.nest), REML=T, data=d13) summary(m1) fig1.data$value[10]=0.4927/(0.517+0.4927+1.1707) fig1.data$value[11]=0.517/(0.517+0.4927+1.1707) fig1.data$value[12]=1.1707/(0.517+0.4927+1.1707) # Peak cort m2=lmer(PeakCort~1 + (1|NatalNest/Host.nest), REML=T, data=d13) summary(m2) fig1.data$value[13]=0.2423/(0.2423+25.734) fig1.data$value[14]=0 fig1.data$value[15]=25.734/(0.2423+25.734) # Delta cort m3=lmer(DeltaCort~1 + (1|NatalNest/Host.nest), REML=T, data=d13) summary(m3) fig1.data$value[16]=0 fig1.data$value[17]=0 fig1.data$value[18]=1 # Make figure 2 - variance components analysis colnames(fig1.data)=c("Cohort", "Component", "Measure", "Value") fig2=ggplot(fig1.data, aes(x=Measure, y=Value, fill=Component)) + geom_bar(stat="identity", colour="black") + scale_x_discrete(limits=c("Baseline", "Peak", "Delta")) + scale_fill_manual(values=c("red", "black", "white")) + facet_wrap(~Cohort) + theme_bw() + ylab("Proportion of variation") + xlab("CORT variable") fig2 png("Figure2-variance components.png", res=300, units="in", width=6, height=3) print(fig2) dev.off() ##### Effects of developmental treatment on CORT variables ##### # 2012 birds (Models 1-3 in Table 3) m1=lmer(Cort0mins~as.factor(Treatment) + (1|Natal_nest), REML=FALSE, data=d12) summary(m1) drop1(m1, test="Chisq") m2=lmer(peakcort~as.factor(Treatment) + Cort0mins + (1|Natal_nest), REML=FALSE, data=d12) summary(m2) drop1(m2, test="Chisq") m3=lmer(deltacort~as.factor(Treatment) + Cort15mins + (1|Natal_nest), REML=FALSE, data=d12) summary(m3) drop1(m3, test="Chisq") # 2013 birds (Models 4-6 in Table 3) colnames(d13) m1=lmer(Cort0~Treatment + (1|NatalNest), REML=FALSE, data=d13) summary(m1) drop1(m1, test="Chisq") m2=lmer(PeakCort~Treatment + Cort0 + (1|NatalNest), REML=FALSE, data=d13) summary(m2) drop1(m2, test="Chisq") m3=lmer(DeltaCort~Treatment + Cort15+ (1|NatalNest), REML=FALSE, data=d13) summary(m3) drop1(m3, test="Chisq") ##### Effects of developmental telomere attrition on CORT variables ############# # 2012 birds (Models 7-9 in Table 3) m0=lmer(DN.verhulst4.55~TLday4 + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$TLday4))) summary(m0) drop1(m0, test="Chisq") # TLday4 does not predict attrition between d4 and d55 m7=lmer(Cort0mins~DN.verhulst4.55 + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$DN.verhulst4.55))) summary(m7) drop1(m7, test="Chisq") m8=lmer(peakcort~DN.verhulst4.55 + Cort0mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$DN.verhulst4.55))) summary(m8) drop1(m8, test="Chisq") m9=lmer(deltacort~DN.verhulst4.55 + Cort15mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(DN.verhulst4.55))) summary(m9) drop1(m9, test="Chisq") # 2013 birds (Models 10-12 in Table 3) m10=lmer(Cort0~TSTWO3 + VerhulstD3to12 + (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) summary(m10) drop1(m10, test="Chisq") m11=lmer(PeakCort~VerhulstD3to12 + TSTWO3 + Cort0 + (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) summary(m11) drop1(m11, test="Chisq") m12=lmer(DeltaCort~VerhulstD3to12 + TSTWO3 + Cort15+ (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) summary(m12) drop1(m12, test="Chisq") #### Meta-analysis of the two cohorts ###### # First rerun the above models scaling the variables so that effect sizes are comparable # baseline m7s=lmer(scale(Cort0mins)~scale(DN.verhulst4.55) + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$DN.verhulst4.55))) summary(m7s) m10s=lmer(scale(Cort0)~TSTWO3 + scale(VerhulstD3to12) + (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) summary(m10s) #peak m8s=lmer(scale(peakcort)~scale(DN.verhulst4.55) + Cort0mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$DN.verhulst4.55))) summary(m8s) m11s=lmer(scale(PeakCort)~scale(VerhulstD3to12) + TSTWO3 + Cort0 + (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) summary(m11s) # delta m9s=lmer(scale(deltacort)~scale(DN.verhulst4.55) + Cort15mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(DN.verhulst4.55))) summary(m9s) m12s=lmer(scale(DeltaCort)~scale(VerhulstD3to12) + TSTWO3 + Cort15+ (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) summary(m12s) # make a dataframe with the effect sizes and standard errors from the above models # (I can't for the life of me figure out how to extract the standard errors from the model objects, so need to type them in...) cohort = c("2012 cohort", "2013 cohort") base.effect = c(-0.0010855,0.08077) base.se = c(0.2139518,0.19004) peak.effect = c(0.47586,0.08453) peak.se = c(0.17249,0.19131) delta.effect = c(0.42041,0.43480) delta.se = c(0.19150,0.17109) meta = data.frame(cohort, base.effect, base.se, peak.effect,peak.se,delta.effect,delta.se) # Use a fixed effects model: asks how large is the average true effect in the two studies included. # questions: should I use weighted or unweighted estimation to fit the model? I can make an argument for weighting the two cohorts equally. # baseline m.base = rma(yi = base.effect, sei = base.se, measure = "GEN", method = "FE", data = meta, weighted = FALSE) summary(m.base) fitstats.rma(m.base) forest(m.base, slab = meta$cohort, xlab="Standardised beta", mlab="Summary estimate", annotate = T) grid.text("Effect of D on baseline cort", .5, .75, gp=gpar(cex=2)) # peaks m.peak = rma(yi = peak.effect, sei = peak.se, measure = "GEN", method = "FE", data = meta, weighted = F) summary(m.peak) fitstats(m.peak) forest(m.peak, slab = meta$cohort, xlab="Standardised beta", mlab="Summary estimate", annotate = F) grid.text("Effect of D on peak cort", .5, .75, gp=gpar(cex=2)) # deltacort m.delta = rma(yi = delta.effect, sei = delta.se, measure = "GEN", method = "FE", data = meta, weighted = F) summary(m.delta) fitstats(m.delta) forest(m.delta, slab = meta$cohort, xlab="Standardised beta", mlab="Summary estimate", annotate = F) grid.text("Effect of D on delta cort", .5, .75, gp=gpar(cex=2)) # output the plots for making a combined forest plot (Figure 4) win.metafile("forest_base.wmf") forest(m.base, slab = meta$cohort, xlab="Standardised beta", mlab="Summary estimate", annotate = T) dev.off() win.metafile("forest_peak.wmf") forest(m.peak, slab = meta$cohort, xlab="Standardised beta", mlab="Summary estimate", annotate = T) dev.off() win.metafile("forest_delta.wmf") forest(m.delta, slab = meta$cohort, xlab="Standardised beta", mlab="Summary estimate", annotate = T) dev.off() #### Model comparison for alternative telomere measures (Table 4) ##### # 2012 birds cor.test(d12$TLOneYear, d12$DN.verhulst4.55, use="complete.obs") cor.test(d12$TLOneYear, d12$DN.verhulst4.365, use="complete.obs") cor.test(d12$DN.verhulst4.55, d12$DN.verhulst4.365, use="complete.obs") cor.test(d12$TLOneYear, d12$TLday55, use="complete.obs") # Peak a1=lmer(peakcort~DN.verhulst4.55 + Cort0mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$DN.verhulst4.55))) summary(a1) AICc(a1) a1b=lmer(peakcort~TLOneYear + Cort0mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$DN.verhulst4.55))) summary(a1b) AICc(a1b)-AICc(a1) a1c=lmer(peakcort~DN.verhulst4.365 + Cort0mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$DN.verhulst4.55))) summary(a1c) AICc(a1c)-AICc(a1) a1d=lmer(peakcort~as.factor(Treatment) + Cort0mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(d12$DN.verhulst4.55))) summary(a1d) AICc(a1d)-AICc(a1) AICc(a1);AICc(a1b);AICc(a1c);AICc(a1d) evidence(aictab(c(a1, a1b, a1c, a1d)), model.low="Mod2") evidence(aictab(c(a1, a1b, a1c, a1d)), model.low="Mod3") evidence(aictab(c(a1, a1b, a1c, a1d)), model.low="Mod4") # Delta am1a=lmer(deltacort~DN.verhulst4.55 + Cort15mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(DN.verhulst4.55))) summary(am1a) AICc(am1a) am1b=lmer(deltacort~TLOneYear + Cort15mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(DN.verhulst4.55))) summary(am1b) AICc(am1b)-AICc(am1a) am1c=lmer(deltacort~DN.verhulst4.365 + Cort15mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(DN.verhulst4.55))) summary(am1c) AICc(am1c)-AICc(am1a) am1d=lmer(deltacort~as.factor(Treatment) + Cort15mins + (1|Natal_nest), REML=FALSE, data=subset(d12, !is.na(DN.verhulst4.55))) summary(am1c) AICc(am1d)-AICc(am1a) AICc(am1a);AICc(am1b);AICc(am1c);AICc(am1d) evidence(aictab(c(am1a, am1b, am1c, am1d)), model.low="Mod2") evidence(aictab(c(am1a, am1b, am1c, am1d)), model.low="Mod3") evidence(aictab(c(am1a, am1b, am1c, am1d)), model.low="Mod4") # 2013 birds # Delta cor.test(d13$VerhulstD3to12, d13$VerhulstD3to1yr, use="complete.obs") cor.test(d13$VerhulstD3to1yr, d13$TS1yr, use="complete.obs") cor.test(d13$VerhulstD3to12, d13$TS1yr, use="complete.obs") cor.test(d13$TSTWO12, d13$TS1yr, use="complete.obs") m1a=lmer(DeltaCort~VerhulstD3to12 + TSTWO3 + Cort15+ (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) AICc(m1a) m1b=lmer(DeltaCort~TS1yr + Cort15+ (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) summary(m1b) AICc(m1b)-AICc(m1a) m1c=lmer(DeltaCort~VerhulstD3to1yr + TSTWO3 + Cort15+ (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) AICc(m1c)-AICc(m1a) m1d=lmer(DeltaCort~Treatment + Cort15+ (1|NatalNest), REML=FALSE, data=subset(d13, !is.na(VerhulstD3to12))) AICc(m1d)-AICc(m1a) AICc(m1a);AICc(m1b);AICc(m1c);AICc(m1d) evidence(aictab(c(m1a, m1b, m1c, m1d)), model.low="Mod2") evidence(aictab(c(m1a, m1b, m1c, m1d)), model.low="Mod3") evidence(aictab(c(m1a, m1b, m1c, m1d)), model.low="Mod4") plot(d13$VerhulstD3to1yr~d13$TS1yr) ### Make figure 3 - main results ##### # 2012 birds error bar d12$split=(d12$DN.verhulst4.55>median(d12$DN.verhulst4.55, na.rm=TRUE)) d12$split l12=melt(d12[ ,c("BirdID", "split", "Cort0mins", "Cort15mins", "Cort30mins")], id=c("BirdID", "split")) l12$Time[l12$variable=="Cort0mins"] = 0 l12$Time[l12$variable=="Cort15mins"] = 15 l12$Time[l12$variable=="Cort30mins"] = 30 l12.2=summarize(group_by(l12, split, Time), mean=mean(value, na.rm=TRUE), sd=sd(value, na.rm=TRUE), length=length(value)) l12.2=l12.2[!is.na(l12.2$split), ] l12.2$se=l12.2$sd/sqrt(l12.2$length) l12.2$Group[l12.2$split==TRUE]="Less attrition" l12.2$Group[l12.2$split==FALSE]="More attrition" Fig3e=ggplot(l12.2, aes(x=as.factor(Time), y=mean, colour=Group, group=Group)) + scale_colour_manual(values=c("black", "darkred")) + scale_fill_manual(values=c("black", "white")) + geom_errorbar(aes(ymax=mean+se, ymin=mean-se), width=0.05, colour="darkgrey") + geom_line(aes(linetype=Group)) + geom_point(shape=21, aes(fill=Group), size=3) + theme_bw() + xlab("Time (min)") + ylab("CORT (ng/ml)") + annotate("text", x=0.7, y=60, label="(e)", size=6) + theme(legend.position="none") #theme(legend.justification=c(0,1), legend.position=c(0,1)) Fig3e # 2012 scatter plot for peak Fig3a = ggplot(data=d12, aes(x=DN.verhulst4.55, y=peakcort)) + geom_point(stat="identity", position="dodge") + geom_smooth(method="lm", colour = "darkred") + xlab(expression(paste(Delta, "TL"))) + ylab("Peak CORT (ng/ml)") + annotate("text", x=-0.8, y=148, label="(a)", size=6) + theme_bw() Fig3a # 2012 scatterplot for delta cort Fig3c = ggplot(data=d12, aes(x=DN.verhulst4.55, y=deltacort)) + geom_point(stat="identity", position="dodge") + geom_smooth(method="lm", colour = "darkred") + xlab(expression(paste(Delta, "TL"))) + ylab(expression(paste(Delta, "CORT (ng/ml)"))) + annotate("text", x=-0.8, y=110, label="(c)", size=6) + theme_bw() Fig3c # 2013 error bar plot d13$split=(d13$VerhulstD3to12>median(d13$VerhulstD3to12, na.rm=TRUE)) l13=melt(d13[ ,c("BirdID", "split", "Cort0", "Cort15", "Cort30")], id=c("BirdID", "split")) l13$Time[l13$variable=="Cort0"] = 0 l13$Time[l13$variable=="Cort15"] = 15 l13$Time[l13$variable=="Cort30"] = 30 l13.2=summarize(group_by(l13, split, Time), mean=mean(value, na.rm=TRUE), sd=sd(value, na.rm=TRUE), length=length(value)) l13.2=l13.2[!is.na(l13.2$split), ] l13.2$se=l13.2$sd/sqrt(l13.2$length) l13.2$Group[l13.2$split==TRUE]="Less attrition" l13.2$Group[l13.2$split==FALSE]="More attrition" Fig3f=ggplot(l13.2, aes(x=as.factor(Time), y=mean, colour=Group, group=Group)) + scale_colour_manual(values=c("black", "darkred")) + scale_fill_manual(values=c("black", "white")) + geom_errorbar(aes(ymax=mean+se, ymin=mean-se), colour="darkgrey", width=0.05) + geom_line(aes(linetype=Group)) + geom_point(shape=21, aes(fill=Group), size=3) + theme_bw() + xlab("Time (min)") + ylab("CORT (ng/ml)") + annotate("text", x=0.7, y=18, label="(f)", size=6) + theme(legend.justification=c(1,0), legend.position=c(1,0)) Fig3f # 2013 scatterplot for peak cort Fig3b = ggplot(data=d13, aes(x=VerhulstD3to12, y=PeakCort)) + geom_point(stat="identity", position="dodge")+ geom_smooth(method="lm", colour = "darkred")+ xlab(expression(paste(Delta, "TL"))) + ylab("Peak CORT (ng/ml)") + annotate("text", x=-0.8, y=31, label="(b)", size=6) + theme_bw() Fig3b # 2013 scatterplot for delta cort Fig3d = ggplot(data=d13, aes(x=VerhulstD3to12, y=DeltaCort)) + geom_point(stat="identity", position="dodge")+ geom_smooth(method="lm", colour = "darkred")+ xlab(expression(paste(Delta, "TL"))) + ylab(expression(paste(Delta, "CORT (ng/ml)"))) + annotate("text", x=-0.8, y=11, label="(d)", size=6) + theme_bw() Fig3d # Put figure 3 together png("Figure3-main results.png", res=300, units="in", width=12, height=8) grid.arrange(Fig3a, Fig3c, Fig3e, Fig3b, Fig3d, Fig3f, ncol=3, nrow=2) dev.off()