## Last run in R v 3.0.2 GUI 1.62 Snow Leopard Build 04/2014 ## Set up R environment before analysis ## Set working directory to folder with data file setwd("") ## Load the data. 'Taff_Data.csv' has all of the data organized with one observed breeding attempt ## in each row. Subsets of data are created here for first year and after second year males ## for separate analyses. New columns are made for some binomial variables that are saved ## as text rather than 0/1. Testosterone values are natural log transformed. Oxidative stress measures ## are arcsin sqrt transformed as in Freeman-Gallant 2011. d<-read.csv("Taff_Data.csv") d$Age2<-ifelse(d$Age=="FY",0,ifelse(d$Age=="SY",1,NA)) d$Surv<-ifelse(d$Survive=="Yes",1,ifelse(d$Survive=="No",0,NA)) d$Imp2<-ifelse(d$Implant=="T",1,ifelse(d$Implant=="C",0,NA)) d$Test_Band2<-log(d$Test_Band+1) d$Test_Recap2<-log(d$Test_Recap+1) d$OXS.a<-asin(sqrt(d$OXS/100))*100 d$OXSRecap.a<-asin(sqrt(d$OXSRecap/100))*100 d.fy<-subset(d,d$Age=="FY") d.sy<-subset(d,d$Age=="SY") d.sy<-subset(d.sy,abs(d.sy$syMask)>0) d.age<-subset(d,d$Age2!="NA") ## Load packages that contain functions that will be used for analysis later on. ## These packages must be downloaded from CRAN to make script work library(rethinking) ## Version 1.30 library(lme4) ## Version 1.0-5 library(plotrix) ## Version 3.5-1 library(Hmisc) ## Version 3.12-2 library(lmerTest) ## Version 2.0-6 library(car) ## Version 2.0-19 ## Should plots be sent to R device or saved to tiff files? plot.to.dev<-1 # set as 1 (plot on screen) or 0 (plot to tiff) ## Table 1. Relationship between plumage ornaments and testosterone at banding. ## Testosterone measure used is natural log tranformed. Results are presented for ## a single linear mixed model that includes ornaments, male age, and two way interactions ## between male age and ornaments. Male identity is included as a random effect. P-values ## are provided by the lmerTest package based on Satterthwaite approximation for df. Table 1 ## is just a summary of the output from this model. ## Full model ## Fit with raw testosterone concentrations m.int<-lmer(Test_Band~syUVB+syYelB+syCCar+syBib+syMask+Age2+ Age2*syUVB+Age2*syYelB+Age2*syCCar+Age2*syBib+Age2*syMask+(1|Male),data=d) ## Check residuals - not normal shapiro.test(resid(m.int)) ## p <0.0001, not normal at all par(mfrow=c(1,2)) plot(resid(m.int)~fitted(m.int)) qqPlot(resid(m.int)) ## Fit with natural log transformed testosterone instead m.int<-lmer(Test_Band2~syUVB+syYelB+syCCar+syBib+syMask+Age2+ Age2*syUVB+Age2*syYelB+Age2*syCCar+Age2*syBib+Age2*syMask+(1|Male),data=d) ## Check residuals - OK shapiro.test(resid(m.int)) ## p = 0.14, proceed with this response par(mfrow=c(1,2)) plot(resid(m.int)~fitted(m.int)) qqPlot(resid(m.int)) ## Step One: Remove Age2*syMask m.int.1<-lmer(Test_Band2~syUVB+syYelB+syCCar+syBib+syMask+Age2+ Age2*syUVB+Age2*syYelB+Age2*syCCar+Age2*syBib+(1|Male),data=d) ## Step Two: Remove Age2*syCCar m.int.2<-lmer(Test_Band2~syUVB+syYelB+syCCar+syBib+syMask+Age2+ Age2*syUVB+Age2*syYelB+Age2*syBib+(1|Male),data=d) ## Step Three: Remove syCCar m.int.3<-lmer(Test_Band2~syUVB+syYelB+syBib+syMask+Age2+ Age2*syUVB+Age2*syYelB+Age2*syBib+(1|Male),data=d) ## Step Four: Remove syMask = Final model reported in paper m.int.4<-lmer(Test_Band2~syUVB+syYelB+syBib+Age2+ Age2*syUVB+Age2*syYelB+Age2*syBib+(1|Male),data=d) ## Check residuals - OK (all were ok on intermediate models too) shapiro.test(resid(m.int.4)) ## p = 0.2, par(mfrow=c(1,2)) plot(resid(m.int.4)~fitted(m.int.4)) qqPlot(resid(m.int.4)) ## Draw samples from the final model for plotting and confidence intervals ## Set the number of samples to be taken from the fit model. post.size<-100000 ## Sample from the multivariate normal distribution of parameter estimates. post<-mvrnorm(n=post.size,mu=fixef(m.int.4),Sigma=vcov(m.int.4)) ## Table 2. Point estimates for testosterone-plumage relationships for each experience class. table.x<-matrix(nrow=3,ncol=6) colnames(table.x)<-c("fyUVB","fyYelB","fyBib","syUVB","syYelB","syBib") rownames(table.x)<-c("mean","low","high") table.x[1,]<-c(mean(post[,2]),mean(post[,3]),mean(post[,4]), mean(post[,2]+post[,6]),mean(post[,3]+post[,7]),mean(post[,4]+post[,8])) table.x[2:3,]<-c(PCI(post[,2]),PCI(post[,3]),PCI(post[,4]), PCI(post[,2]+post[,6]),PCI(post[,3]+post[,7]),PCI(post[,4]+post[,8])) ## Table 3 pearson correlations for plumage ornaments to each other. ## Each pair of correlations, pooling males d2<-as.data.frame(d[,32:36]) d.tab3<-subset(d2,d2$syMask!="NA") cor(d.tab3,use="complete.obs") d3<-as.matrix(d.tab3) rcorr(d3) ## Figure 1. Commands to plot model predictions from table 1 for the panels included in Figure 1. ## Note that all of these values derive from the fit model above in table 1. ## Link function to generate predictions for different paramenter combinations from fit model. mu.link<-function(U,Y,B,A){exp( post[,1] + post[,2]*U + post[,3]*Y + post[,4]*B + post[,5]*A + post[,6]*U*A + post[,7]*Y*A + post[,8]*B*A) } ## Sequence of values over which to generate predictions for plotting (note that predictors are ## standardized, so this range represents three standard deviations above and below the mean). r<-seq(from=-3.4,to=3.4,by=0.05) ## Mean values of each predictor. All plumage ornaments but one are held at these values for plotting. mean.uvb<-mean(d$syUVB) mean.yelb<-mean(d$syYelB) mean.ccar<-mean(d$syCCar) mean.bib<-mean(na.omit(d$syBib)) mean.mask<-mean(na.omit(d$syMask)) ## Set up plotting commands if(plot.to.dev==1){ tiff("Fig1New.tiff",width=10.2,height=3.8,units="in",res=300,compression="lzw") } par(mfrow=c(1,3)) ## Plot UV Brightness relationship by age. First panel of figure 1. ## Calculate mean and ci across sequence for first year males. mu.uvb.fy<-sapply(r,function(r)mean(mu.link(r,mean.yelb,mean.bib,0))) ci.uvb.fy<-sapply(r,function(r)PCI(mu.link(r,mean.yelb,mean.bib,0))) ## Calculate mean and ci across sequence for second year males. mu.uvb.sy<-sapply(r,function(r)mean(mu.link(r,mean.yelb,mean.bib,1))) ci.uvb.sy<-sapply(r,function(r)PCI(mu.link(r,mean.yelb,mean.bib,1))) ## Make the plot plot(r,mu.uvb.fy,xlim=c(-3,3),type="n",ylim=c(0,11),xlab="UV Brightness (SD from mean)", ylab="Testosterone (ng/mL)") lines(r,mu.uvb.fy,lwd=1,col="slateblue") lines(r,mu.uvb.sy,lwd=1) shade(ci.uvb.fy,r,col=col.alpha("slateblue",0.25)) shade(ci.uvb.sy,r,col=col.alpha("black",0.4)) corner.label("a") ## Plot Yellow Brightness relationship by age. Second panel of figure 1. ## Calculate mean and ci across sequence for first year males. mu.yelb.fy<-sapply(r,function(r)mean(mu.link(mean.uvb,r,mean.bib,0))) ci.yelb.fy<-sapply(r,function(r)PCI(mu.link(mean.uvb,r,mean.bib,0))) ## Calculate mean and ci across sequence for second year males. mu.yelb.sy<-sapply(r,function(r)mean(mu.link(mean.uvb,r,mean.bib,1))) ci.yelb.sy<-sapply(r,function(r)PCI(mu.link(mean.uvb,r,mean.bib,1))) ## Make the plot plot(r,mu.yelb.fy,xlim=c(-3,3),type="n",ylim=c(0,11),xlab="Yellow Brightness (SD from mean)", ylab="Testosterone (ng/mL)") lines(r,mu.yelb.fy,lwd=1,col="slateblue") lines(r,mu.yelb.sy,lwd=1) shade(ci.yelb.fy,r,col=col.alpha("slateblue",0.25)) shade(ci.yelb.sy,r,col=col.alpha("black",0.4)) corner.label("b") ## Plot Bib Size relationship by age. Third panel of figure 1. ## Calculate mean and ci across sequence for first year males. mu.bib.fy<-sapply(r,function(r)mean(mu.link(mean.uvb,mean.yelb,r,0))) ci.bib.fy<-sapply(r,function(r)PCI(mu.link(mean.uvb,mean.yelb,r,0))) ## Calculate mean and ci across sequence for second year males. mu.bib.sy<-sapply(r,function(r)mean(mu.link(mean.uvb,mean.yelb,r,1))) ci.bib.sy<-sapply(r,function(r)PCI(mu.link(mean.uvb,mean.yelb,r,1))) ## Make the plot plot(r,mu.bib.fy,xlim=c(-3,3),type="n",ylim=c(0,11),xlab="Bib Size (SD from mean)", ylab="Testosterone (ng/mL)") lines(r,mu.bib.fy,lwd=1,col="slateblue") lines(r,mu.bib.sy,lwd=1) shade(ci.bib.fy,r,col=col.alpha("slateblue",0.25)) shade(ci.bib.sy,r,col=col.alpha("black",0.4)) corner.label("c") if(plot.to.dev==1){dev.off()} ## Figure 2. GLMM with binomial fit of survival predicted by plasma testosterone ## concentration. Male ID is included as a random effect as there are many ## males with multiple observations. Plot shows the raw points indicating ## testosterone levels for males that did or did not survive and lines show ## the predicted likelihood of survival at each level of testosterone from ## the fit model. ## Conditionally plot to tiff file if(plot.to.dev==1){ tiff("Figure2new.tiff",width=5,height=6.2,units="in",res=300) } ## Run model and plot model predictions and raw data data.surv<-subset(d.age,d.age$Test_Band2!="NA") m<-glmm(Surv~Test_Band+Age2+(1|Male),data=data.surv,family="binomial") post.s<-mvrnorm(n=post.size,mu=fixef(m),Sigma=vcov(m)) r<-seq(from=-1,to=11,by=0.1) plot(d$Test_Band,d$Surv,pch=16,ylab="Likelihood of Surviving", xlab="Testosterone (ng/mL)",col=col.alpha("black",alpha=0.5),xlim=c(0,10)) mean.s<-sapply(r,function(z)logistic(mean(post.s[,1]+post.s[,2]*z+post.s[,3]*mean(na.omit(d$Age2))))) lines(r,mean.s) ci.s<-sapply(r,function(z)logistic(HPDI(post.s[,1]+post.s[,2]*z+post.s[,3]*mean(na.omit(d$Age2))))) shade(ci.s,r,col=col.alpha("black",0.2)) ## Conditionally turn off plotting device if(plot.to.dev==1){dev.off()} ## Figure 3. Effect of testosterone or control implants on change in TAC, change in OXS ## and survival percentage. ## Conditionally plot to tiff instead of r device if(plot.to.dev==1){ tiff("Figure3.tiff",width=9.3,height=5.4,units="in",res=300) } ## Plotting parameters par(mfrow=c(1,3)) ## Panel 1 ## Calculate means and SE t.m<-mean(na.omit(subset(d$OXSChange,d$Implant=="T"))) sd(na.omit(subset(d$OXSChange,d$Implant=="T"))) t.se<-sd(na.omit(subset(d$OXSChange,d$Implant=="T"))/sqrt(14)) c.m<-mean(na.omit(subset(d$OXSChange,d$Implant=="C"))) sd(na.omit(subset(d$OXSChange,d$Implant=="C"))) c.se<-sd(na.omit(subset(d$OXSChange,d$Implant=="C")))/sqrt(14) ## Make the plot barplot(c(c.m,t.m),names.arg=c("Control","Testosterone"), ylab="Change in % DNA Damage",ylim=c(-1,5)) abline(h=0,lwd=1.5) corner.label(label="a") ## Add error bars to plot error<-0.05 lines(rep(0.7,2),c(c.m+c.se,c.m-c.se)) lines(c(0.7-error,0.7+error),rep(c.m+c.se,2)) lines(c(0.7-error,0.7+error),rep(c.m-c.se,2)) lines(rep(1.9,2),c(t.m+t.se,t.m-t.se)) lines(c(1.9-error,1.9+error),rep(t.m+t.se,2)) lines(c(1.9-error,1.9+error),rep(t.m-t.se,2)) ## Test for difference using linear model with TAC at recapture as response ## and both TAC at banding and implant as predictors d.oxs<-subset(d,d$Implant!="") m<-lm(OXSRecap~OXS+Age2+Imp2+CaptoCapDays,data=d.oxs) summary(m) ## Check residuals of model - OK #shapiro.test(resid(m)) #par(mfrow=c(1,2)) #plot(resid(m),fitted(m)) #qqPlot(resid(m)) ## Put p-value on plot lines(c(0.7,1.9),rep(4.2,2)) text(1.3,4.4,"P = 0.82") ## Panel 2. Implants and change in TAC. Plot shows mean TAC change for T vs. C ## implanted males. P-value from t-test ## Calculate means and SE t.m<-mean(na.omit(subset(d$TACChange,d$Implant=="T"))) sd(na.omit(subset(d$TACChange,d$Implant=="T"))) t.se<-sd(na.omit(subset(d$TACChange,d$Implant=="T"))/sqrt(14)) c.m<-mean(na.omit(subset(d$TACChange,d$Implant=="C"))) sd(na.omit(subset(d$TACChange,d$Implant=="C"))) c.se<-sd(na.omit(subset(d$TACChange,d$Implant=="C")))/sqrt(14) ## Make the plot barplot(c(c.m,t.m),names.arg=c("Control","Testosterone"), ylab="Change in TAC (mM)",ylim=c(-.5,.6)) abline(h=0,lwd=1.5) corner.label(label="b") ## Add error bars to plot error<-0.05 lines(rep(0.7,2),c(c.m+c.se,c.m-c.se)) lines(c(0.7-error,0.7+error),rep(c.m+c.se,2)) lines(c(0.7-error,0.7+error),rep(c.m-c.se,2)) lines(rep(1.9,2),c(t.m+t.se,t.m-t.se)) lines(c(1.9-error,1.9+error),rep(t.m+t.se,2)) lines(c(1.9-error,1.9+error),rep(t.m-t.se,2)) ## Test for difference using linear model with TAC at recapture as response ## and both TAC at banding and implant as predictors d.tac<-subset(d,d$Implant!="") m<-lm(TACRecap~TACBand+Imp2+Age2+CaptoCapDays,data=d.tac) summary(m) ## Check residuals of model - OK #shapiro.test(resid(m)) #par(mfrow=c(1,2)) #plot(resid(m),fitted(m)) #qqPlot(resid(m)) ## Put p-value on plot lines(c(0.7,1.9),rep(0.5,2)) text(1.3,0.53,"P = 0.02") ## Panel 3. Implants and survival. Plot just shows percentage of each treatment ## group that survived to next year. Sample size for implant birds is very ## small. Can show fisher exact test p-value but note that effect size is ## large even though result is ns. barplot(c(54,58,33),names.arg=c("None","C","T"), ylim=c(0,100),ylab="Percent Surviving to Next Year") abline(h=0,lwd=2) corner.label(label="c") ## Calculate mean and SE of the mean for each treatment group survival mean(na.omit(subset(d$Surv,d$Implant=="T"))) sd(na.omit(subset(d$Surv,d$Implant=="T"))) sd(na.omit(subset(d$Surv,d$Implant=="T")))/sqrt(9) mean(na.omit(subset(d$Surv,d$Implant=="C"))) sd(na.omit(subset(d$Surv,d$Implant=="C"))) sd(na.omit(subset(d$Surv,d$Implant=="C")))/sqrt(12) mean(na.omit(subset(d$Surv,d$Implant==""))) sd(na.omit(subset(d$Surv,d$Implant==""))) sd(na.omit(subset(d$Surv,d$Implant=="")))/sqrt(72) ## Add error bars to the plot error<-0.05 lines(rep(0.7,2),c(54-5.9,54+5.9)) lines(c(0.7-error,0.7+error),rep(54+5.9,2)) lines(c(0.7-error,0.7+error),rep(54-5.9,2)) lines(rep(1.9,2),c(58-14.9,58+14.9)) lines(c(1.9-error,1.9+error),rep(58+14.9,2)) lines(c(1.9+error,1.9-error),rep(58-14.9,2)) lines(rep(3.1,2),c(33-16.7,33+16.7)) lines(c(3.1-error,3.1+error),rep(33+16.7,2)) lines(c(3.1-error,3.1+error),rep(33-16.7,2)) ## Add sample sizes height<-3.5 text(0.7,height,"n = 72") text(1.9,height,"n = 12") text(3.1,height,"n = 9") ## Fisher exact tests comparing survival for control vs. t implants x<-matrix(nrow=2,ncol=2) x[1,]<-c(7,5) x[2,]<-c(3,6) fisher.test(x,alternative="greater") ## Add test statistic lines(c(1.9,3.1),rep(80,2)) text(2.5,83.5,"P = 0.25") ## Conditionally turn off plotting device if(plot.to.dev==1){dev.off()} ## Other tests reported in text but not in figures or tables ## Observational dataset ## Does testosterone level differ between experienced and inexperienced males? m<-lmer(Test_Band2~Age2+(1|Male),data=d.age) summary(m) ## Get summary stats for each age level in normal scale mean(na.omit(subset(d$Test_Band,d$Age2==0))) sd(na.omit(subset(d$Test_Band,d$Age2==0))) mean(na.omit(subset(d$Test_Band,d$Age2==1))) sd(na.omit(subset(d$Test_Band,d$Age2==1))) ## Check residuals for normality each age level resi<-cbind(resid(m),subset(d.age$Age2,d.age$Test_Band>0)) ## inexperienced males - borderline in S-W, but not significant shapiro.test(subset(resi[,1],resi[,2]==0)) qqPlot(subset(resi[,1],resi[,2]==0)) ## experienced males shapiro.test(subset(resi[,1],resi[,2]==1)) qqPlot(subset(resi[,1],resi[,2]==1)) ## Is testosterone related to OXS? m<-lmer(OXS.a~Test_Band+Age2+JDateBand+(1|Male),data=d.age) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## By age separately - only referenced in paper, numbers not reported ## FY m<-lm(OXS.a~Test_Band+JDateBand,data=d.fy) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## SY m<-lmer(OXS.a~Test_Band+JDateBand+(1|Male),data=d.sy) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## Is testosterone related to TAC m<-lmer(TACBand~Test_Band2+Age+JDateBand+(1|Male),data=d.age) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## By age separately - only referenced in paper, numbers not reported ## FY m<-lm(TACBand~Test_Band2,data=d.fy) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## SY m<-lmer(TACBand~Test_Band2+(1|Male),data=d.sy) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## Is endogeneous testosterone related to apparent overwinter survival m<-glmm(Surv~Test_Band2+Age2+(1|Male),data=d.age,family="binomial") summary(m) ## Experimental dataset d$CaptoCapDays2<-(d$CaptoCapDays-mean(na.omit(d$CaptoCapDays)))/sd(na.omit(d$CaptoCapDays)) ## Effect of implants on testosterone at recapture m<-lm(Test_Recap2~Test_Band2+Imp2+Age2+CaptoCapDays,data=d) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## Effect of implants on weight, controlling for initial weight d$WeightBand2<-(d$WeightBand-mean(na.omit(d$WeightBand)))/sd(na.omit(d$WeightBand)) d$WeightRecap2<-(d$WeightRecap-mean(na.omit(d$WeightRecap)))/sd(na.omit(d$WeightRecap)) m<-lm(WeightRecap2~WeightBand2+Imp2+Age2+CaptoCapDays,data=d) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## Effect of implants on TAC, controlling for initial TAC d$TACBand2<-(d$TACBand-mean(na.omit(d$TACBand)))/sd(na.omit(d$TACBand)) d$TACRecap2<-(d$TACRecap-mean(na.omit(d$TACRecap)))/sd(na.omit(d$TACRecap)) m<-lm(TACRecap2~TACBand2+Imp2+Age2+CaptoCapDays2,data=d) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m)) ## Effect of implants on OXS, controlling for initial OXS d$OXSRecap2<-(d$OXSRecap-mean(na.omit(d$OXSRecap)))/sd(na.omit(d$OXSRecap)) d$OXS2<-(d$OXS-mean(na.omit(d$OXS)))/sd(na.omit(d$OXS)) d.oxs<-subset(d,d$Implant!="") m<-lm(OXSRecap2~OXS2+Age2+Imp2+CaptoCapDays,data=d.oxs) summary(m) ## Check residuals shapiro.test(resid(m)) par(mfrow=c(1,2)) plot(resid(m)~fitted(m)) qqPlot(resid(m))