library(gdata) library(rptR) library(survival) library(MuMIn) library(sem) library(plot3D) GWA<-read.xls("GWA-DataforEvolution.xlsx",sheet = "TrialScores", 1) GWAave<-read.xls("GWA-DataforEvolution.xlsx",sheet = "AverageScores", 1) #GWA stands for "The Great Walkabout", the name we somehow #ended up giving to this project. #repeatability of soft songs and wing waves rpt.aov(GWA$Ssrate,GWA$Bird,npermut=1000) rpt.aov(GWA$Wwrate,GWA$Bird,npermut=1000) #Confirmatory factor analysis on average scores #putative aggressive behaviors: rates of flights, closest approach, time spent within 5m #putative aggressive signals: rates of soft songs and wing waves GWAaveResp<-GWAave[,2:6] #take the response variables as above GWAavecov<-cov(GWAaveResp) #calculate the covariance matrix #specify model 1: single latent variable gwa.model.ave1<-specify.model() F1 -> AveWW, lam1, NA F1 -> AveSS, lam2, NA F1 -> AveFlight, lam3, NA F1 -> AveTime5, lam4, NA F1 -> AveClosest, lam5, NA AveWW <-> AveWW, e1, NA AveSS <-> AveSS, e2, NA AveFlight <-> AveFlight, e3, NA AveTime5 <-> AveTime5, e4, NA AveClosest <-> AveClosest, e5, NA F1 <-> F1, NA, 1 #specify model 1: two latent variables gwa.model.ave2<-specify.model() F1 -> AveWW, lam1, NA F1 -> AveSS, lam2, NA F2 -> AveFlight, lam3, NA F2 -> AveTime5, lam4, NA F2 -> AveClosest, lam5, NA AveWW <-> AveWW, e1, NA AveSS <-> AveSS, e2, NA AveFlight <-> AveFlight, e3, NA AveTime5 <-> AveTime5, e4, NA AveClosest <-> AveClosest, e5, NA F1 <-> F1, NA, 1 F2 <-> F2, NA, 1 F1 <-> F2, F1F2, NA #Confirmatory factor analysis aveGWA.sem<-sem(gwa.model.ave1,GWAavecov,nrow(GWAaveResp)) summary(aveGWA.sem) aveGWA.sem2<-sem(gwa.model.ave2,GWAavecov,nrow(GWAaveResp)) summary(aveGWA.sem2) # Delta AIC: 33.256 (Model 1) - 23.70413 (Model 2)= 9.55187 # Indicates that Model 2 (two latent variables) have a better fit #Lande-Arnold phenotypic selection analysis GWAaveS<-GWAave[!is.na(GWAave$YearsFrom2009),] GWAaveS$zsoft<-(GWAaveS$AveSS-mean(GWAaveS$AveSS))/sd(GWAaveS$AveSS) GWAaveS$zflight<-(GWAaveS$AveFlight-mean(GWAaveS$AveFlight))/sd(GWAaveS$AveFlight) GWAaveS$zclosest<-(GWAaveS$AveClosest-mean(GWAaveS$AveClosest))/sd(GWAaveS$AveClosest) GWAaveS$ztime5<-(GWAaveS$AveTime5-mean(GWAaveS$AveTime5))/sd(GWAaveS$AveTime5) GWAaveS$zAgg<-(GWAaveS$AveAggPCA-mean(GWAaveS$AveAggPCA))/sd(GWAaveS$AveAggPCA) GWAaveS$relfitness<-GWAaveS$YearsFrom2009/mean(GWAaveS$YearsFrom2009) #Lande Arnold model on soft song and aggression PCA scores only (reported in the main paper) LandeArnoldSSAggPCA<-lm(relfitness~zsoft+zAgg+I(zsoft^2)+I(zAgg^2)+zsoft*zAgg,data=GWAaveS) summary(LandeArnoldSSAggPCA) #selection surface for soft song and aggression scores x=seq(-1,3,length=50) y=seq(-2,2,length=50) f=function(x,y) {0.89643+(0.04279*x)-(0.13319*y)+(0.35217*I(x^2))+(0.17481*I(y^2))-(0.79487*x*y)} z=outer(x,y,f) ribbon3D(x,y,z,theta = 00, phi = 00, expand = 0.5, col=ramp.col(c("blue", "yellow", "red","black")), colvar=z, colkey=TRUE, ltheta = 150, along="xy", space= 0.5, shade = 0.1,ticktype="detailed", xlab = "soft songs (z-scores)", ylab = "aggression score (z-scores)", zlab = "survival") GWAaveS$zx<-1:50 for(i in 1:50) zx[i]<-mean(z[i,]) plot(type="l",x,zx,xlab="soft song scores (z-scores)", ylab="predicted survival") #3-D scatterplot of actual data. scatter3D(GWAaveS$zsoft,GWAaveS$zAgg,GWAaveS$YearsFrom2009, theta=30, phi= 30,bty="b" ,xlab = "soft songs (z-scores)", ylab = "aggression score (z-scores)", zlab = "survival" , col=ramp.col(c("blue", "yellow", "red","black")),ticktype = "detailed") #Additional Analyses on soft songs and individual aggressive behaviors #all three models yield significant negative correlational selection #between soft songs and aggressive behaviors #Lande-Arnold on soft songs and time w 5 only LandeArnoldSSTime<-lm(relfitness~zsoft+ztime5+I(zsoft^2)+I(ztime5^2)+zsoft*ztime5,data=GWAaveS) summary(LandeArnoldSSTime) #Lande-Arnold on soft song and flights only LandeArnoldSSflight<-lm(relfitness~zsoft+zflight+I(zsoft^2)+I(zflight^2)+zsoft*zflight,data=GWAaveS) summary(LandeArnoldSSflight) #Lande-Arnold on soft song and closest only LandeArnoldSSclose<-lm(relfitness~zsoft+zclosest+I(zsoft^2)+I(zclosest^2)+zsoft*zclosest,data=GWAaveS) summary(LandeArnoldSSclose) #survival analysis (Cox Regression) with age as a covariate GWAaveS$startyear= 0 S<-Surv(time=GWAaveS$startyear,time2=I(GWAaveS$YearsFrom2009+1),event=GWAaveS$survivedTo2014) SurvModelSSPCA0<-coxph(S~AgeIn2009+AveAggPCA*AveSS+I(AveSS^2)+I(AveAggPCA^2) ,data=GWAaveS, na.action="na.fail", subset= !is.na(AgeIn2009)) summary(SurvModelSSPCA0) SurvModelAve<-dredge(SurvModelSSPCA0) print(subset(SurvModelAve, delta<2), abbrev=FALSE) #best model AveragedSurvModelSS<-model.avg(SurvModelAve, delta<2) summary(AveragedSurvModelSS) #averaged model shows a significant interaction effect #is conserved when age is added as a covariate in the Cox-Regression model.