######################################################################## ######################################################################## ##Midwest Ecotraps Analysis: fungal load change versus temp ##Copyright 2020 Skylar Hopkins #This analysis determines whether changes in fungal loads from November to March #were affected by November roosting temperatures and/or fungal loads. We use simple #correlation tests and a Logan-10 growth curve, fit in a Bayesian framework, to explore #these relationships. The Bayesian analysis uses weakly informative priors which allow #for any possible curve shape, as shown by the simulation included at the end of the script. #This script also makes Figure 2 from the paper. #Permission is hereby granted, free of charge, to any person obtaining a copy of this #software and associated documentation files (the "Software"), to deal in the Software #without restriction, including without limitation the rights to use, copy, modify, merge, #publish, distribute, sublicense, and/or sell copies of the Software, and to permit #persons to whom the Software is furnished to do so, subject to the following conditions: #The above copyright notice and this permission notice shall be included in all #copies or substantial portions of the Software. #THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR #PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR #ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, #ARISING FROM, OUT OF OR IN #CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ######################################################################## ######################################################################## #clear workspace rm(list=ls()) ######################################################################## #################load packages and data########################## ######################################################################## library(lme4) library(ggplot2) library(reshape2) library(brglm) library(MASS) library(tidyverse) library(viridis) library(naniar) library(ciTools) library(dplyr) library(nlme) library(predictmeans) library(mgcv) library(brms) library(pROC) library(cowplot) library(R2jags) library(lattice) library(mcmcplots) #Cleaned MYLU swab data dm <- read.csv("~/Dropbox/Hopkins_VT_postdoc/HopkinsGit/Midwest EcoTraps Paper/Cleaned MYLU Swab Data for Ecotraps Analysis.csv") #Cleaned MYLU Band data MYLUBands <- read.csv("~/Dropbox/Hopkins_VT_postdoc/HopkinsGit/Midwest EcoTraps Paper/Cleaned MYLU Band Data for Ecotraps Analysis.csv") #Cleaned Inf MYLU Band/Recap data MYLUBands_RecapAnalysis <- read.csv("~/Dropbox/Hopkins_VT_postdoc/HopkinsGit/Midwest EcoTraps Paper/Cleaned Inf MYLU Band Recap Data for Ecotraps Analysis.csv") table(dm$site) length(dm$swab_id[dm$ysw < 1 & dm$season=="hiber_earl"]) #573 length(dm$swab_id[dm$ysw %in% c(1, 2) & dm$season=="hiber_earl"]) #521 length(dm$swab_id[dm$ysw > 2 & dm$season=="hiber_earl"]) #178 ##################################################################################### ############Calculate the fungus growth rate on each bat################################ ################################################################################### #Calculate number of days between fall and spring samples MYLUBands_RecapAnalysis$DiffTime<-as.numeric(difftime(as.Date(MYLUBands_RecapAnalysis$date.y,"%m/%d/%y"), as.Date(MYLUBands_RecapAnalysis$date.x, "%m/%d/%y"), units = c("days"))) #Log lambda calculated "daily" - not great distributions MYLUBands_RecapAnalysis$LogLambda<-log((MYLUBands_RecapAnalysis$lategdL+0.0001)/(MYLUBands_RecapAnalysis$earlygdL+0.0001))^(1/MYLUBands_RecapAnalysis$DiffTime) c=ggplot(MYLUBands_RecapAnalysis, aes(x=earlytemp)) + geom_jitter(width = NULL, size=5, alpha=0.9, aes(y=LogLambda, col=as.factor(wyear.x))) + ylab("Fungus log lambda")+ xlab(expression("Early roosting temperature (degrees C)"))+ theme_bw()+ theme(panel.grid = element_blank(),axis.title=element_text(size=20), axis.text=element_text(size=15), axis.line=element_line(), legend.position="right",legend.text = element_text(size=20,face="italic"),strip.text = element_text(size=25,face="italic")) c #Log lambda calculately "monthly" (28 days) - better MYLUBands_RecapAnalysis$LogLambda<-log(((MYLUBands_RecapAnalysis$lategdL+0.0001)/(MYLUBands_RecapAnalysis$earlygdL+0.0001))^(1/(MYLUBands_RecapAnalysis$DiffTime/28))) c=ggplot(MYLUBands_RecapAnalysis, aes(x=earlytemp)) + geom_jitter(width = NULL, size=5, alpha=0.9, aes(y=LogLambda, col=as.factor(wyear.x))) + #scale_colour_viridis(option="inferno")+ ylab("Fungus log lambda")+ xlab(expression("Early roosting temperature (degrees C)"))+ theme_bw()+ theme(panel.grid = element_blank(),axis.title=element_text(size=20), axis.text=element_text(size=15), axis.line=element_line(), legend.position="right",legend.text = element_text(size=20,face="italic"),strip.text = element_text(size=25,face="italic")) c #Log lambda calculately "weekly" - also better. This is what I went with MYLUBands_RecapAnalysis$LogLambda<-log(((MYLUBands_RecapAnalysis$lategdL+0.0001)/(MYLUBands_RecapAnalysis$earlygdL+0.0001))^(1/(MYLUBands_RecapAnalysis$DiffTime/7))) c=ggplot(MYLUBands_RecapAnalysis, aes(x=earlytemp)) + geom_jitter(width = NULL, size=5, alpha=0.9, aes(y=LogLambda, col=as.factor(wyear.x))) + #scale_colour_viridis(option="inferno")+ ylab("Fungus log lambda")+ xlab(expression("Early roosting temperature (degrees C)"))+ theme_bw()+ theme(panel.grid = element_blank(),axis.title=element_text(size=20), axis.text=element_text(size=15), axis.line=element_line(), legend.position="right",legend.text = element_text(size=20,face="italic"),strip.text = element_text(size=25,face="italic")) c ##################################################################################### ########Run sinple analyses to confirm positive relationship######################### ################################################################################### #We did not to include a random site effect. Recaps probs do vary by site, #but mostly because ave temps vary by site table(MYLUBands_RecapAnalysis$site.x, MYLUBands_RecapAnalysis$ysw.x) SiteRecaps <- MYLUBands_RecapAnalysis %>% group_by(site.x) %>% summarise(NBats=n(), NRecaps=sum(EverRecapturedYN.x), PropRecaptured=NRecaps/NBats, meantemp=mean(earlytemp, na.rm=TRUE)) hist(SiteRecaps$PropRecaptured) plot(SiteRecaps$PropRecaptured~SiteRecaps$meantemp) #early and late temps are correlated and can't go in same model plot(MYLUBands_RecapAnalysis$latetemp~MYLUBands_RecapAnalysis$earlytemp) #Fit a simple linear model and GAM - including logearlyloads as a predictor #GAM confirms that relationship is nonlinear w/ the "peak" near the max observed temps test2<-lm(LogLambda~earlytemp + logearlyloads, data=MYLUBands_RecapAnalysis); summary(test2) test3<-gam(LogLambda~s(earlytemp) + logearlyloads, data=MYLUBands_RecapAnalysis, method = "REML"); summary(test3) cor.test(MYLUBands_RecapAnalysis$earlytemp, MYLUBands_RecapAnalysis$LogLambda, method="spearman") cor.test(MYLUBands_RecapAnalysis$logearlyloads, MYLUBands_RecapAnalysis$LogLambda, method="spearman") nd <- expand.grid(LogLambda=NA,earlytemp=seq(5,12,0.1), logearlyloads=seq(-5, 0, 1)) nd$phat <- predict(test2, nd, type="response", re.form=~0) c=ggplot(MYLUBands_RecapAnalysis, aes(x=earlytemp, col=logearlyloads)) + geom_jitter(width = NULL, size=5, alpha=0.9, aes(y=LogLambda)) + geom_line(data=nd[nd$logearlyloads==-1,], aes(x=earlytemp,y=phat),size=1)+ geom_line(data=nd[nd$logearlyloads==-2,], aes(x=earlytemp,y=phat),size=2)+ geom_line(data=nd[nd$logearlyloads==-3,], aes(x=earlytemp,y=phat),size=3)+ geom_line(data=nd[nd$logearlyloads==-4,], aes(x=earlytemp,y=phat),size=4)+ scale_x_continuous(limits =c(5,12), breaks=seq(5,12,1)) + scale_colour_viridis(option="inferno")+ ylab("Fungus log lambda")+ xlab(expression("Early roosting temperature (degrees C)"))+ theme_bw()+ theme(panel.grid = element_blank(),axis.title=element_text(size=20), axis.text=element_text(size=15), axis.line=element_line(), legend.position="right",legend.text = element_text(size=20,face="italic"),strip.text = element_text(size=25,face="italic")) c c=ggplot(MYLUBands_RecapAnalysis, aes(x=latetemp, col=logearlyloads)) + geom_jitter(width = NULL, size=5, alpha=0.9, aes(y=LogLambda)) + scale_colour_viridis(option="inferno")+ ylab("Fungus log lambda")+ xlab(expression("Late roosting temperature (degrees C)"))+ theme_bw()+ theme(panel.grid = element_blank(),axis.title=element_text(size=20), axis.text=element_text(size=15), axis.line=element_line(), legend.position="right",legend.text = element_text(size=20,face="italic"),strip.text = element_text(size=25,face="italic")) c ##################################################################################### ########Fit the Logan-10 in a Bayesian framework using JAGS############################ ################################################################################### ##Define model model <- " model { #Priors for fixed effects - a ~ dunif(0.419, 1.211) ##a shape parameter g ~ dunif(1.995, 19.459) ##a shape parameter p ~ dunif(0.156, 0.819) ##a shape parameter int ~ dunif(-0.5, 0.5) ##intercept parameter Topt ~ dunif(12.965, 15.162) ##Thermal optimum - based on 95% CI from Verant sigma ~ dunif(0, 100) # standard deviation tau <- 1 / (sigma * sigma) # sigma^2 doesn't work in JAGS #Likelihoods for (i in 1:NBats) { y[i] ~ dnorm(mu[i], tau) #tau is precision (1 / variance) mu[i] <- int*Load[i] + a*(1/(1+g*exp(-p*Temp[i]))-exp(-(21.5-Temp[i])/(21.5-Topt))) #mu[i] <- a*(1/(1+g*exp(-p*Temp[i]))-exp(-(21.5-Temp[i])/(21.5-Topt))) } }" #Bundle data to pass to Jags MYLUBands_RecapAnalysis<-MYLUBands_RecapAnalysis[!is.na(MYLUBands_RecapAnalysis$earlytemp),] NBats<-length(MYLUBands_RecapAnalysis$LogLambda) Temp<-MYLUBands_RecapAnalysis$earlytemp Load<-MYLUBands_RecapAnalysis$logearlyloads y<-MYLUBands_RecapAnalysis$LogLambda data<-list(y=y, Temp=Temp, Load=Load, NBats=NBats) ##Inits function inits<-function(){list(a=runif(n = 1, min = 0.45, max = 1.2), g=runif(n = 1, min = 2.1, max = 19), p=runif(n = 1, min = 0.16, max = 0.79), int=runif(n=1, min=-0.1, max=0.1), Topt=runif(n = 1, min = 12.95, max=15.21))} ##Parameters to estimate params<-c("a", "g", "p", "int", "Topt", "sigma") #Run the model and pull posterior samples jags.m <- jags(data = data, inits = inits, parameters.to.save = params, n.chains = 3, n.iter = 160000, #niter should be much higher n.burnin = 60000, model.file = textConnection(model)) jags.m #Check output using a plot Tmax<-21.5 Topt_out<-as.vector(jags.m$BUGSoutput$mean$Topt) a_out<-as.vector(jags.m$BUGSoutput$mean$a) g_out<-as.vector(jags.m$BUGSoutput$mean$g) p_out<-as.vector(jags.m$BUGSoutput$mean$p) int_out<-as.vector(jags.m$BUGSoutput$mean$int) Curve<-as.data.frame(expand.grid(earlytemp=seq(4.5, 22, 0.01), logearlyload=seq(-4, -0.5, by=0.5))) table(round(MYLUBands_RecapAnalysis$logearlyloads, 0)) Curve$NumBats<-121 Curve$NumBats[round(Curve$logearlyload, 0)==-3]<-92 Curve$NumBats[round(Curve$logearlyload, 0)==-2]<-33 Curve$NumBats[round(Curve$logearlyload, 0)==-1]<-10 #Can use these if I want to test how params change curve shape #a_out<-1 #g_out<-40 #p_out<-0.5 #Topt_out<-14 #int_out<-0 Curve$Pred<-int_out*Curve$logearlyload+a_out*(((1+g_out*exp(-p_out*Curve$earlytemp))^-1)-exp(-(Tmax-Curve$earlytemp)/(Tmax-Topt_out))) #Curve$Pred<-a_out*(((1+g_out*exp(-p_out*Curve$earlytemp))^-1)-exp(-(Tmax-Curve$earlytemp)/(Tmax-Topt_out))) a=ggplot(MYLUBands_RecapAnalysis, aes(x=earlytemp, col=logearlyloads)) + geom_line(data=Curve[Curve$logearlyload==-4,], aes(y=Pred, col=logearlyload, size=NumBats), alpha=0.4)+ geom_line(data=Curve[Curve$logearlyload==-3,], aes(y=Pred, col=logearlyload, size=NumBats), alpha=0.4)+ geom_line(data=Curve[Curve$logearlyload==-2,], aes(y=Pred, col=logearlyload, size=NumBats), alpha=0.4)+ geom_line(data=Curve[Curve$logearlyload==-1,], aes(y=Pred, col=logearlyload, size=NumBats), alpha=0.4)+ geom_line(data=nd[nd$logearlyloads==-4,], aes(x=earlytemp,y=phat),size=2)+ geom_line(data=nd[nd$logearlyloads==-3,], aes(x=earlytemp,y=phat),size=2)+ geom_line(data=nd[nd$logearlyloads==-2,], aes(x=earlytemp,y=phat),size=2)+ geom_line(data=nd[nd$logearlyloads==-1,], aes(x=earlytemp,y=phat),size=2)+ #data geom_point(size=3,aes(y=LogLambda, col=logearlyloads), alpha=0.9) + geom_hline(yintercept = 0, linetype="dashed", col="gray")+ scale_colour_viridis(breaks=seq(-4, 0, 1), begin=0.3, end=1, option="inferno")+ ylab(expression(Fungal~growth~rate~(lambda)))+ xlab(expression(Early~roosting~temperature~(degree~C))) + theme_bw()+ labs(color=expression(italic(Pd)~load~(log[10]))) a #Calculate an R2 for fit MYLUBands_RecapAnalysis$PredictedFungusLambda<-int_out*MYLUBands_RecapAnalysis$logearlyloads + a_out*(((1+g_out*exp(-p_out*MYLUBands_RecapAnalysis$earlytemp))^-1)-exp(-(Tmax-MYLUBands_RecapAnalysis$earlytemp)/(Tmax-Topt_out))) MYLUBands_RecapAnalysis$Resid<-MYLUBands_RecapAnalysis$LogLambda-MYLUBands_RecapAnalysis$PredictedFungusLambda plot(Resid~as.factor(site.x), data=MYLUBands_RecapAnalysis) points(Resid~as.factor(site.x), data=MYLUBands_RecapAnalysis) plot(MYLUBands_RecapAnalysis$LogLambda~MYLUBands_RecapAnalysis$PredictedFungusLambda) R2test<-lm(MYLUBands_RecapAnalysis$LogLambda~MYLUBands_RecapAnalysis$PredictedFungusLambda); summary(R2test) #Is any residual variation explained by late temp - no hist(MYLUBands_RecapAnalysis$Resid) #BEAUTIFUL abline(v=0, col="red") plot(MYLUBands_RecapAnalysis$Resid~MYLUBands_RecapAnalysis$latetemp) abline(h=0, lty=2) ##################################################################################### ##################################PLOTS############################################ ################################################################################### MYLUBands_RecapAnalysis$NumBats<-round(MYLUBands_RecapAnalysis$logearlyloads, 0) MYLUBands_RecapAnalysis$NumBats[MYLUBands_RecapAnalysis$NumBats==-4]<-121 MYLUBands_RecapAnalysis$NumBats[MYLUBands_RecapAnalysis$NumBats==-3]<-92 MYLUBands_RecapAnalysis$NumBats[MYLUBands_RecapAnalysis$NumBats==-2]<-33 MYLUBands_RecapAnalysis$NumBats[MYLUBands_RecapAnalysis$NumBats==-1]<-10 MYLUBands_RecapAnalysis$NumBats[MYLUBands_RecapAnalysis$NumBats==0]<-10 titlesize<-12 textsize<-12 legendtitle<-guide_legend(bquote(log[10]~italic(Pd)~"load"~"(n=num. banded bats)")) a=ggplot(MYLUBands_RecapAnalysis, aes(x=earlytemp)) + geom_rect(xmin=12.9, xmax=15.146, ymin=-0.2, ymax=0.62, fill="grey80", linetype=0) + geom_text(x=(12.9+(15.146-12.9)/2), y=0.6, label="Topt", col="black")+ geom_rect(xmin=20.53, xmax=22.52, ymin=-0.2, ymax=0.62, fill="grey80", linetype=0) + geom_text(x=(20.53+(22.52-20.53)/2), y=0.6, label="Tmax", col="black")+ geom_hline(yintercept = 0, linetype="dashed", col="gray")+ geom_line(data=Curve[Curve$logearlyload==-4,], aes(y=Pred, col=NumBats, size=NumBats), alpha=0.4)+ geom_line(data=Curve[Curve$logearlyload==-3,], aes(y=Pred, col=NumBats, size=NumBats), alpha=0.4)+ geom_line(data=Curve[Curve$logearlyload==-2,], aes(y=Pred, col=NumBats, size=NumBats), alpha=0.4)+ geom_line(data=Curve[Curve$logearlyload==-1,], aes(y=Pred, col=NumBats, size=NumBats), alpha=0.4)+ geom_point(size=3,aes(y=LogLambda, col=NumBats), alpha=0.9) + #geom_line(data=nd, aes(x=earlytemp,y=phat),size=2)+ scale_y_continuous(limits = c(-0.1, 0.6)) + scale_x_continuous(limits = c(4, 22)) + scale_colour_viridis(breaks=c(30, 60, 90, 120), option="B", direction = -1, begin=0.2, end=0.8, labels=c("-1 (n=11)", "-2 (n=33)", "-3 (n=92)", "-4 (n=121)"))+ scale_size_continuous(breaks=c(30, 60, 90, 120), labels=c("-1 (n=11)", "-2 (n=33)", "-3 (n=92)", "-4 (n=121)")) + #geom_errorbar(data=JohnsonLewin, aes(x=Temp, ymin=(Coef-2*SE), ymax=(Coef+2*SE)), width=0.2, col="grey20") + #geom_point(data=JohnsonLewin, aes(x=Temp, y=Coef), pch=17, size=2, col="grey20") + ylab(expression(Over~winter~fungal~load~increase~(log10~lambda))) + xlab(expression(Early~roosting~temperature~(degree~C))) + theme_classic()+ theme(panel.grid = element_blank(), axis.title=element_text(size=titlesize), axis.text=element_text(size=textsize), legend.position=c(0.21, 0.86), legend.title = element_text(size=10), legend.key.height = unit(0.1, units="in"), legend.box = "horizontal", legend.background = element_rect(fill="transparent"), legend.text = element_text(size=10)) + guides(color=legendtitle, size=legendtitle) a ggsave(file=paste("Fig2.pdf", sep = ""),width=6.5,height=5,units="in",limitsize=FALSE) ##################################################################################### #########Simulate how much flexibility there is in priors########################### ################################################################################### #This visualizes how much flexibility in curve shapes our priors allow - basically any curve shape #is possible. This plot has tons of lines on it, so only run if you have some time and computing power. Sim<-as.data.frame(expand.grid(a=seq(0.4, 1.2, 0.2), g=seq(2, 20, 2), p=seq(0.2, 0.8, 0.2), Topt=seq(13, 15, 1), int=c(-0.2, 0, -0.2))) Sim<-as.data.frame(expand.grid(a=c(0.42, 0.815, 1.21), g=c(2, 10.75, 19.5), p=c(0.16, 0.49, 0.82), Topt=seq(13, 15, 1), int=c(-0.5, 0, -0.5))) dim(Sim) Curve<-as.data.frame(expand.grid(earlytemp=seq(0, 22, 0.1), logearlyloads=seq(-4, -0.5, 0.5), pred=NA)) Curve$pred<-Sim$int[1]*Curve$logearlyloads + Sim$a[1]*(((1+Sim$g[1]*exp(-Sim$p[1]*Curve$earlytemp))^-1)-exp(-(21.5-Curve$earlytemp)/(21.5-Sim$Topt[1]))) c=ggplot(MYLUBands_RecapAnalysis, aes(x=earlytemp)) + geom_point(size=3,aes(y=LogLambda, col=logearlyloads), alpha=0.9) + geom_line(data=Curve[Curve$logearlyload==-4,], aes(y=pred, col=logearlyloads), alpha=0.2)+ geom_line(data=Curve[Curve$logearlyload==-3,], aes(y=pred, col=logearlyloads), alpha=0.2)+ geom_line(data=Curve[Curve$logearlyload==-2,], aes(y=pred, col=logearlyloads), alpha=0.2)+ geom_line(data=Curve[Curve$logearlyload==-1,], aes(y=pred, col=logearlyloads), alpha=0.2)+ scale_y_continuous(limits = c(-0.1, 0.6)) + scale_x_continuous(limits = c(4, 22)) + scale_colour_viridis(option="B", begin=0.2, end=0.8)+ ylab(expression(Fungal~growth~rate~(Delta~log10~ng~DNA)))+ xlab(expression(Early~roosting~temperature~(degree~C))) + theme_classic()+ theme(panel.grid = element_blank(), axis.title=element_text(size=titlesize), axis.text=element_text(size=textsize), legend.position="none", legend.title = element_text(size=10), legend.key.height = unit(0.1, units="in"), legend.box = "horizontal", legend.background = element_rect(fill="transparent"), legend.text = element_text(size=10)) c for(i in 1:length(Sim$a)) { Curve$pred<-Sim$int[i]*Curve$logearlyloads + Sim$a[i]*(((1+Sim$g[i]*exp(-Sim$p[i]*Curve$earlytemp))^-1)-exp(-(21.5-Curve$earlytemp)/(21.5-Sim$Topt[i]))) c <- c + geom_line(data=Curve[Curve$logearlyload==-4,], aes(y=pred, col=logearlyloads), alpha=0.2)+ geom_line(data=Curve[Curve$logearlyload==-3,], aes(y=pred, col=logearlyloads), alpha=0.2)+ geom_line(data=Curve[Curve$logearlyload==-2,], aes(y=pred, col=logearlyloads), alpha=0.2)+ geom_line(data=Curve[Curve$logearlyload==-1,], aes(y=pred, col=logearlyloads), alpha=0.2) } ggsave(plot=c, file=paste(getwd(), "/Midwest EcoTraps Paper/PriorSims.png", sep = ""),width=6.5,height=5,units="in",limitsize=FALSE) ## Counts<-c(9419, 172, 409, 0) Criteria<-c(1, 2, 3, 4) Data<-as.data.frame(cbind(Criteria, Counts)) X<-ggplot(data=Data, aes(x=Criteria, y=Counts)) + geom_col() + xlab("Number of criteria met") + ylab("Number of simulations") + theme_classic() ggsave(plot=X, file=paste(getwd(), "/Midwest EcoTraps Paper/RandomLoss_SigSims.png", sep = ""),width=6.5,height=5,units="in",limitsize=FALSE)