####PIPFIL Telomere Analysis Graphs ##March 27th 2019 #Ben Vernasco library(dplyr) library(ggplot2) library(lmerTest) library(lubridate) library(visreg) library(lubridate) library(gridExtra) DF<-read.csv("PIPFIL.Telo.Data.csv") levels(DF$TH.FL)<-c("Floater","Territory-holder") DF$Count<-rep(1,187) DF$Field.Season<-as.factor(DF$Field.Season) ##Age-based changes in Telomere Lengths summary(lm1<-lmer(rTL~(Avg.Age+Delta.Age)*TH.FL+(1|Band.Num), data = DF)) #Graphing age-related changes a<-visreg(lm1, "Avg.Age", by = "TH.FL") aaa<-a$fit DF$part.res<-a$res[,6] aaa<-aaa[-c(49:101),] ab<-ggplot(filter(aaa, TH.FL != "Floater"))+ aes(x = Avg.Age, y = visregFit)+ geom_line(size = 1)+ geom_point(aes(x = Avg.Age, y = visregRes),shape = 1,size = 2, data = filter(a$res, TH.FL != "Floater"))+ theme_classic(base_size = 14)+ #geom_smooth(aes(x = Avg.Age, y = rTL),method = "lm", data = DF, color = "black")+ xlab("Average Age (years)")+ ylab("")+ theme(strip.background = element_blank(),strip.text.x = element_text(size = 14))+ ggtitle("Territory-holders")+ theme(plot.title = element_text(hjust = 0.5))+ scale_x_continuous(limits = c(0,16), breaks = c(0,4,8,12,16))+ scale_y_continuous(limits = c(0.8,1.5), breaks = c(0.8,1,1.2,1.4)) ac<-ggplot(filter(aaa, TH.FL == "Floater"))+ aes(x = Avg.Age, y = visregFit)+ geom_line(size = 1, linetype = "dashed")+ geom_point(aes(x = Avg.Age, y = visregRes),shape = 1,size = 2, data = filter(a$res, TH.FL == "Floater"))+ theme_classic(base_size = 14)+ #geom_smooth(aes(x = Avg.Age, y = rTL),method = "lm", data = DF, color = "black")+ xlab("Average Age (years)")+ ylab("Relative Telomere Lengths")+ theme(strip.background = element_blank(),strip.text.x = element_text(size = 14))+ ggtitle("Floaters")+ theme(plot.title = element_text(hjust = 0.5))+ scale_x_continuous(limits = c(0,8), breaks = c(0,2,4,6,8))+ scale_y_continuous(limits = c(0.8,1.5), breaks = c(0.8,1,1.2,1.4)) grid.arrange(ac,ab, ncol = 2) ###Comparing behavior models using AICc aictable<-function(X,m){ # # This function will create a full model selection table based on AICc. # # Inputs: # # X is the AIC table produced by R by using AIC(model1,model2, ...) # m is the sample size # # rnames<-row.names(X) AICc<-X$AIC+2*X$df*(X$df+1)/(m-X$df-1) #small-sample correction logL<-X$df-X$AIC/2 #Log-likelihood tab<-data.frame(X[,1],logL,AICc) #Remove AIC column; add logL and AICc colnames(tab)[1]<-c("Params") #Rename "df" column row.names(tab)<-rnames tab<-tab[order(tab$AICc),] #Sort by ascending AICc value deltaAICc<-tab$AICc-min(tab$AICc) #Delta AICc weight<-exp(-deltaAICc/2)/sum(exp(-deltaAICc/2)) #Weights cumwt<-weight #Column for cumulative weight for(i in 2:dim(X)[1]){ cumwt[i]<-cumwt[i-1]+cumwt[i] #Accumulate weight from the top } tab<-data.frame(tab,deltaAICc,weight,cumwt) tab<-round(tab,4) tab } #degree DF2<-arrange(DF,degree) DF2<-DF2[-c(87:187),] DF2$degree<-log(DF2$degree+1) summary(lm8<-lmer(degree~rTL+Age.at.Cap2+Field.Season+(1|Band.Num), data = DF2)) summary(lm9<-lmer(degree~rTL+Field.Season+(1|Band.Num), data = DF2)) summary(lm10<-lmer(degree~(rTL+Age.at.Cap2)*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm11<-lmer(degree~rTL*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm12<-lmer(degree~Age.at.Cap2+Field.Season+(1|Band.Num), data = DF2)) summary(lm13<-lmer(degree~Age.at.Cap2*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm14<-lmer(degree~Field.Season+(1|Band.Num), data = DF2)) summary(lm15<-lmer(degree~rTL+Age.at.Cap2+Field.Season+TH.FL+(1|Band.Num), data = DF2)) summary(lm16<-lmer(degree~rTL+Field.Season+TH.FL+(1|Band.Num), data = DF2)) summary(lm17<-lmer(degree~Age.at.Cap2+Field.Season+TH.FL+(1|Band.Num), data = DF2)) X<-AIC(lm8,lm9,lm10,lm11,lm12,lm13,lm14,lm15,lm16,lm17) m<-86 aictable(X,m) b<-visreg(lm9, "rTL", scale = "response") bbb<-b$fit DF2$part.res<-b$res[,5] ab<-ggplot(bbb)+ aes(x = rTL, y = visregFit)+ geom_point(aes(x = rTL, y = part.res), shape = 1, size = 2, data = DF2)+ #geom_line(size = 1)+ theme_classic(base_size = 12)+ xlab("Relative Telomere Length")+ ylab("Degree")+ theme(plot.margin = margin(0.5,0.5,0.5,2, "cm")) #strength DF2<-arrange(DF,strength) DF2<-DF2[-c(87:187),] DF2$strength<-log(DF2$strength+1) summary(lm8<-lmer(strength~rTL+Age.at.Cap2+Field.Season+(1|Band.Num), data = DF2)) summary(lm9<-lmer(strength~rTL+Field.Season+(1|Band.Num), data = DF2)) summary(lm10<-lmer(strength~(rTL+Age.at.Cap2)*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm11<-lmer(strength~rTL*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm12<-lmer(strength~Age.at.Cap2+Field.Season+(1|Band.Num), data = DF2)) summary(lm13<-lmer(strength~Age.at.Cap2*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm14<-lmer(strength~Field.Season+(1|Band.Num), data = DF2)) summary(lm15<-lmer(strength~rTL+Age.at.Cap2+Field.Season+TH.FL+(1|Band.Num), data = DF2)) summary(lm16<-lmer(strength~rTL+Field.Season+TH.FL+(1|Band.Num), data = DF2)) summary(lm17<-lmer(strength~Age.at.Cap2+Field.Season+TH.FL+(1|Band.Num), data = DF2)) X<-AIC(lm8,lm9,lm10,lm11,lm12,lm13,lm14,lm15,lm16,lm17) m<-86 aictable(X,m) b<-visreg(lm9, "rTL", scale = "response") bbb<-b$fit DF2$part.res<-b$res[,5] ac<-ggplot(bbb)+ aes(x = rTL, y = visregFit)+ geom_point(aes(x = rTL, y = part.res), shape = 1, size = 2, data = DF2)+ geom_line(size = 1)+ theme_classic(base_size = 12)+ xlab("Relative Telomere Length")+ ylab("Strength")+ theme(plot.margin = margin(0.5,0.5,0.5,0.5, "cm")) #effort DF2<-arrange(DF,effort) DF2<-DF2[-c(87:187),] DF2$effort<-log(DF2$effort+1) summary(lm8<-lmer(effort~rTL+Age.at.Cap2+Field.Season+(1|Band.Num), data = DF2)) summary(lm9<-lmer(effort~rTL+Field.Season+(1|Band.Num), data = DF2)) summary(lm10<-lmer(effort~(rTL+Age.at.Cap2)*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm11<-lmer(effort~rTL*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm12<-lmer(effort~Age.at.Cap2+Field.Season+(1|Band.Num), data = DF2)) summary(lm13<-lmer(effort~Age.at.Cap2*TH.FL+Field.Season+(1|Band.Num), data = DF2)) summary(lm14<-lmer(effort~Field.Season+(1|Band.Num), data = DF2)) summary(lm15<-lmer(effort~rTL+Age.at.Cap2+Field.Season+TH.FL+(1|Band.Num), data = DF2)) summary(lm16<-lmer(effort~rTL+Field.Season+TH.FL+(1|Band.Num), data = DF2)) summary(lm17<-lmer(effort~Age.at.Cap2+Field.Season+TH.FL+(1|Band.Num), data = DF2)) X<-AIC(lm8,lm9,lm10,lm11,lm12,lm13,lm14,lm15,lm16,lm17) m<-86 aictable(X ,m) b<-visreg(lm12, "Age.at.Cap2", scale = "response") bbb<-b$fit DF2$part.res<-b$res[,5] ad<-ggplot(bbb)+ aes(x = Age.at.Cap2, y = visregFit)+ geom_point(aes(x = Age.at.Cap2, y = part.res), shape = 1, size = 2, data = DF2)+ geom_line(size = 1)+ theme_classic(base_size = 12)+ xlab("Age at Capture")+ ylab("Effort")+ theme(plot.margin = margin(0.5,0.5,0.5,0.5, "cm")) #importance DF3<-arrange(DF,import) DF3<-DF3[-c(87:210),] summary(lm8<-lmer(import~rTL+Age.at.Cap2+Field.Season+(1|Band.Num), data = DF3)) summary(lm9<-lmer(import~rTL+Field.Season+(1|Band.Num), data = DF3)) summary(lm10<-lmer(import~(rTL+Age.at.Cap2)*TH.FL+Field.Season+(1|Band.Num), data = DF3)) summary(lm11<-lmer(import~rTL*TH.FL+Field.Season+(1|Band.Num), data = DF3)) summary(lm12<-lmer(import~Age.at.Cap2+Field.Season+(1|Band.Num), data = DF3)) summary(lm13<-lmer(import~Age.at.Cap2*TH.FL+Field.Season+(1|Band.Num), data = DF3)) summary(lm14<-lmer(import~Field.Season+(1|Band.Num), data = DF3)) summary(lm15<-lmer(import~rTL+Age.at.Cap2+Field.Season+TH.FL+(1|Band.Num), data = DF3)) summary(lm16<-lmer(import~rTL+Field.Season+TH.FL+(1|Band.Num), data = DF3)) summary(lm17<-lmer(import~Age.at.Cap2+Field.Season+TH.FL+(1|Band.Num), data = DF3)) X<-AIC(lm8,lm9,lm10,lm11,lm12,lm13,lm14,lm15,lm16,lm17) m<-86 aictable(X,m) b<-visreg(lm9, "rTL", scale = "response") bbb<-b$fit DF3$part.res<-b$res[,5] ae<-ggplot(bbb)+ aes(x = rTL, y = visregFit)+ geom_point(aes(x = rTL, y = part.res), shape = 1, size = 2, data = DF3)+ geom_line(size = 1)+ theme_classic(base_size = 12)+ xlab("Relative Telomere Length")+ ylab("Importance")+ theme(plot.margin = margin(0.5,0.5,0.5,2, "cm")) ###Graphing behavior results, network depictions of behaviors found in paper not included grid.arrange(ad,ab,ac,ae) ###Repeatability mod <- lmer(rTL ~ 1 + (1|Band.Num), data=DF) myfun.lmer <- function(mod, group){ tab <- data.frame(VarCorr(mod))[,c('grp','vcov')] return(tab$vcov[tab$grp=='Band.Num']/sum(tab$vcov)) } myfun.lmer(mod) rTL.rep.sim <- bootMer(mod, myfun.lmer, nsim=1000)$t quantile(rTL.rep.sim, c(0.025,0.975))