######################################################### #### Summarizing thermal thresholds for activity #### ######################################################### # Script by Octavio Jiménez-Robles. octavio.jimenez.robles (at) gmail.com # Used in the analyses of Jiménez-Robles & De la Riva. 2019. # Lizards in the mist: Thermal niches constrained by habitat and microclimates in the Andes of southern Bolivia. # Journal of Biogeography. #--------------------------------------------------------------------------------------------------------------- # 1. Defining Thermal thresholds for activity #### # The maximum threshold whe chose is the average of the voluntary thermal maxima # The minimum threshold we chose is the 10th percentile of the body temperatures recorded in the field in active individuals physioData <- read.table("input_data/PhysiologyData.csv", header = T, sep=";") species.fisio <- data.frame(Species=c("puna", "ornatus", "sp", "orientalis"), N.Tb=NA, TEmerge=NA, N.Tp=NA, LowerTset=NA, UpperTset=NA, TVmax=NA, TVmax2=NA, N.TbP=NA, TEmergeP=NA, N.TpP=NA, LowerTsetP=NA, UpperTsetP=NA, TVmaxP=NA, TVmaxP2=NA) for(i in 1:length(species.fisio$Species)){ speciesdata <- subset(physioData, Species==species.fisio$Species[i]) species.fisio$N.Tb[i] <- nrow(subset(speciesdata, Tb>0)) species.fisio$TEmerge[i] <- quantile(subset(speciesdata)$Tb, 0.10, na.rm=T) species.fisio$N.Tp[i] <- nrow(subset(speciesdata, TpMax>0)) species.fisio$LowerTset[i] <- mean(speciesdata$SPRlow, na.rm=T) species.fisio$UpperTset[i] <- mean(speciesdata$SPRup, na.rm=T) species.fisio$TVmax[i] <- mean(speciesdata$TpMax, na.rm=T) species.fisio$TVmax2[i] <- paste(round(mean(speciesdata$TpMax, na.rm=TRUE),1), " \u00B1 ", round(sd(speciesdata$TpMax, na.rm=TRUE),1), " (",round(min(speciesdata$TpMax, na.rm=TRUE),1), "-", round(max(speciesdata$TpMax, na.rm=T),1), ")", sep="") groupdata <- subset(speciesdata, Sex=='P') species.fisio$N.TbP[i] <- nrow(subset(groupdata, Tb>0 )) species.fisio$TEmergeP[i] <- quantile(subset(groupdata,)$Tb, 0.10, na.rm=T) species.fisio$N.TpP[i] <- nrow(subset(groupdata, TpMax>0)) species.fisio$LowerTsetP[i] <- mean(groupdata$SPRlow, na.rm=T) species.fisio$UpperTsetP[i] <- mean(groupdata$SPRup, na.rm=T) species.fisio$TVmaxP[i] <- mean(groupdata$TpMax, na.rm=T) species.fisio$TVmaxP2[i] <- paste(round(mean(groupdata$TpMax, na.rm=TRUE),1), " \u00B1 ", round(sd(groupdata$TpMax, na.rm=TRUE),1), " (",round(min(groupdata$TpMax, na.rm=TRUE),1), "-", round(max(groupdata$TpMax, na.rm=T),1), ")", sep="") } rm(speciesdata, groupdata, i) View(species.fisio) write.table(species.fisio, "output_data/SpeciesPhysiology.csv", sep=";") # 2. Voluntary Thermal Maxima comparison #### summary(lm(TpMax~Species, data=physioData)) TukeyHSD(aov(TpMax~Species, data=physioData)) summary(lm(TpMax~Species, data=subset(physioData, Sex=="P"))) TukeyHSD(aov(TpMax~Species, data=subset(physioData, Sex=="P"))) # 3. Graphs of Temergence ###### pun<-subset(merged, Species=='puna') orn<-subset(merged, Species=='ornatus') pan<-subset(merged, Species=='pantherinus') ori<-subset(merged, Species=='orientalis') load("/Users/octavio/Dropbox/PhD-dropbox/PhD_Analyses/PhDCleanData/Sama.fisio.RData") par(mfrow=c(2,2), mar=c(.5,.5,.5,.5), oma=c(5,5,4,3)) hist(physioData$Tb[physioData$Species == 'puna'], col='dodgerblue', breaks=15, ylim=c(0,30), xlim=c(10,60), xlab ='', main="") axis(1, labels=F); text(50,25, labels="Tb", cex=1.2); axis(3, at=c(35), labels='Liolaemus puna', las=0, tck=0, lwd=0, font=3) hist(physioData$Tb[physioData$Species == 'puna' & physioData$Sex=="NP" | physioData$Species == 'puna' & physioData$Sex=="P" ], add=T, col=rgb(.3,.3,.3,.5), breaks=15) hist(physioData$Tb[physioData$Species == 'puna' & physioData$Sex=="P"], add=T, col=rgb(0,0,0,.5), breaks=10, lwd=5) abline(v=c(species.fisio$TEmerge[1], species.fisio$TVmax[1] ), col='red', lty=2) rect(species.fisio$LowerTset[1], par('usr')[3], species.fisio$UpperTset[1], par('usr')[4], col=rgb(1,0,0,.3), density=30, lwd=2, border=rgb(1,0,0,.2)) text(55,3, labels=paste0("N = ", species.fisio$N.Tb[1]), cex=.8) hist(physioData$Tb[physioData$Species == 'ornatus'], col='orange', border=NA, breaks=30, ylim=c(0,30), xlim=c(10,60), xlab='', main="", yaxt='n', ylab="") axis(4); axis(1, labels=F); text(50,25, labels="Tb", cex=1.2); axis(3, at=35, labels='Liolaemus ornatus', las=0, tck=0, lwd=0, font=3) hist(physioData$Tb[physioData$Species == 'ornatus' & physioData$Sex=="NP" | physioData$Species == 'ornatus' & physioData$Sex=="P" ], add=T, col=rgb(.3,.3,.3,.5), border=NA, breaks=30) hist(physioData$Tb[physioData$Species == 'ornatus' & physioData$Sex=="P"], add=T, col=rgb(0,0,0,0.5), border=NA,#density=20, angle=135, breaks=30, lwd=5) abline(v=c(species.fisio$TEmerge[2], species.fisio$TVmax[2] ), col='red', lty=2) rect(species.fisio$LowerTset[2], par('usr')[3], species.fisio$UpperTset[2], par('usr')[4], col=rgb(1,0,0,.3), density=30, lwd=2, border=rgb(1,0,0,.2)) text(55,3, labels=paste0("N = ", species.fisio$N.Tb[2]), cex=.8) hist(physioData$Tb[physioData$Species == 'sp'], col='yellow', breaks=10, ylim=c(0,30), xlim=c(10,60), xlab='', main="", ylab="") axis(1, labels=F); text(50,25, labels="Tb", cex=1.2); axis(1, at=35, labels='Liolaemus sp', las=0, tck=0, lwd=0, font=3, line=2) hist(physioData$Tb[physioData$Species == 'sp' & physioData$Sex=="NP" | physioData$Species == 'sp' & physioData$Sex=="P" ], add=T, col=rgb(.3,.3,.3,.5), breaks=10) hist(physioData$Tb[physioData$Species == 'sp' & physioData$Sex=="P"], add=T, col=rgb(0,0,0,0.5), #density=20, angle=135, breaks=5, lwd=5) abline(v=c(species.fisio$TEmerge[3], species.fisio$TVmax[3] ), col='red', lty=2) rect(species.fisio$LowerTset[3], par('usr')[3], species.fisio$UpperTset[3], par('usr')[4], col=rgb(1,0,0,.3), density=30, lwd=2, border=rgb(1,0,0,.2)) text(55,3, labels=paste0("N = ", species.fisio$N.Tb[3]), cex=.8) hist(physioData$Tb[physioData$Species == 'orientalis'], col='forestgreen', breaks=15, ylim=c(0,30), xlim=c(10,60), xlab='', main="", yaxt='n', ylab="") axis(4); axis(1, labels=F); text(50,25, labels="Tb", cex=1.2); axis(1, at=35, labels='Liolaemus orientalis', las=0, tck=0, lwd=0, font=3 line=2) hist(physioData$Tb[physioData$Species == 'orientalis' & physioData$Sex=="NP" | physioData$Species == 'orientalis' & physioData$Sex=="P" ], add=T, col=rgb(.3,.3,.3,.5), breaks=10) hist(physioData$Tb[physioData$Species == 'orientalis' & physioData$Sex=="P"], add=T, col=rgb(0,0,0,0.5), #density=20, angle=135, breaks=3, lwd=5) abline(v=c(species.fisio$TEmerge[4], species.fisio$TVmax[4] ), col='red', lty=2) rect(species.fisio$LowerTset[4], par('usr')[3], species.fisio$UpperTset[1], par('usr')[4], col=rgb(1,0,0,.3), density=30, lwd=2, border=rgb(1,0,0,.2)) text(55,3, labels=paste0("N = ", species.fisio$N.Tb[4]), cex=.8) mtext("Temperature (ºC)", side=1, outer = TRUE, cex = 1, line=3) mtext("Frequency of observations", side=2, outer = TRUE, cex = 1, line=3) # pdf 4 x 8 # png 900 x 400 pero ojo con subíndices # 3. Calibration of operative temperature models #### # Let's check whether every one unit of temperature increase in # body temperature measurement there is a corresponding one unit of temperature in # operative temperature measurement. If the 95% regression bands do not include the identity # expression y = x, then the equation should be used to adjust measurements. data <- read.table("input_data/TevsTbCalibrations.csv", sep=";", header=T) # 3.1. Individual level analyses #### # The calculations have to go in pair of columns, in which the first one is always the lizard and the second one is the model ResultsTable <- as.data.frame(t(data[1:4, seq(2,40, by=2)])) colnames(ResultsTable) <- data[1:4,1] for(i in 1:(length(data)/2)){ i <- (i*2) lizard <- data[1:11, i] lizardT <- as.numeric(as.character(data[5: 172, i])) tubeT <- as.numeric(as.character(data[5: 172, i+1])) ResultsTable$pearson[i/2] <- round(cor(lizardT, tubeT, use="pairwise.complete.obs", "pearson"), 3) model <- lm(lizardT~tubeT) ResultsTable$R2[i/2] <- round(summary(model)$r.squared, 2) ResultsTable$Fvalue[i/2] <- round(summary(model)$fstatistic[1], 1) ResultsTable$DF[i/2] <- paste(summary(model)$fstatistic[2], summary(model)$fstatistic[3], sep=",") ResultsTable$PValue[i/2] <- as.character(pf(summary(model)$fstatistic[1],summary(model)$fstatistic[2], summary(model)$fstatistic[3], lower.tail = FALSE)) ResultsTable$Equation[i/2] <- paste ("y = ", round(summary(model)$coefficients[1,1],2), " + ", round(summary(model)$coefficients[2,1],2), "x", sep="") ResultsTable$Intercept[i/2] <- summary(model)$coefficients[1,1] ResultsTable$InterceptSE[i/2] <- summary(model)$coefficients[1,2] ResultsTable$Slope[i/2] <- summary(model)$coefficients[2,1] ResultsTable$SlopeSE[i/2] <- summary(model)$coefficients[2,2] } rm(lizard, lizardT, tubeT, model, i) ResultsTable <- ResultsTable[order(ResultsTable$Species, ResultsTable$Sex),] write.table(ResultsTable, "output_data/TevsTbCalibrationResults.csv", sep=";") # 3.2. Species level analyses #### individuals <- as.data.frame(t(data[1:4, -1])) # 20 individuals plus their 20 models colnames(individuals) <- data[1:4,1] individuals$SVL <- as.numeric(as.character(individuals$SVL)) class(individuals$SVL) SpeciesResults <- data.frame(Species=levels(individuals$Species), SVL=c(mean(subset(individuals, Species=="orientalis")$SVL, na.rm=T), mean(subset(individuals, Species=="ornatus")$SVL, na.rm=T), mean(subset(individuals, Species=="puna")$SVL, na.rm=T), mean(subset(individuals, Species=="sp")$SVL, na.rm=T))) library(MASS); library(lme4); library(MuMIn); library(merTools) for (i in 1:length(levels(individuals$Species))){ SP <- levels(individuals$Species)[i] lizardColumns <- grep(SP, individuals$Species)+1 SpeciesT <- as.numeric(as.character(data[6:172, lizardColumns[1]])) SpeciesT[168:334] <- as.numeric(as.character(data[6:172, lizardColumns[2]])) SpeciesT[335:501] <- as.numeric(as.character(data[6:172, lizardColumns[3]])) tubeSPT <- as.numeric(as.character(data[6:172, lizardColumns[1]+1])) tubeSPT[168:334] <- as.numeric(as.character(data[6:172, lizardColumns[2]+1])) tubeSPT[335:501] <- as.numeric(as.character(data[6:172, lizardColumns[3]+1])) if(length(lizardColumns)>3){ SpeciesT[502:668] <- as.numeric(as.character(data[6:172, lizardColumns[4]])) tubeSPT[502:668] <- as.numeric(as.character(data[6:172, lizardColumns[4]+1])) } if(length(lizardColumns)>4){ SpeciesT[669:835] <- as.numeric(as.character(data[6:172, lizardColumns[5]])) tubeSPT[669:835] <- as.numeric(as.character(data[6:172, lizardColumns[5]+1])) } if(length(lizardColumns)>5){ SpeciesT[836:1002] <- as.numeric(as.character(data[6:172, lizardColumns[6]])) tubeSPT[836:1002] <- as.numeric(as.character(data[6:172, lizardColumns[6]+1])) } temp <- data.frame(lizardT=SpeciesT, tubeT=tubeSPT, lizardID=c(rep(1, 167), rep(2, 167), rep(3, 167), rep(4,167), rep(5, 167), rep(6,167))[1:length(SpeciesT)], SVL=c(rep(individuals$SVL[lizardColumns[1]-1], 167), rep(individuals$SVL[lizardColumns[2]-1], 167), rep(individuals$SVL[lizardColumns[3]-1], 167), rep(individuals$SVL[lizardColumns[4]-1],167), rep(individuals$SVL[lizardColumns[5]-1], 167), rep(individuals$SVL[lizardColumns[6]-1],167))[1:length(SpeciesT)]) ModelMixedS <- lmer(lizardT ~ tubeT + (1+tubeT|lizardID), data=temp) # both random intercept and random slope ModelNull <- lmer(lizardT ~ 1 + (1+tubeT|lizardID), data=temp) summary(ModelMixedS) SpeciesResults$EquationMixed[i] <- paste ("y = ", round(summary(ModelMixedS)$coefficients[1,1],2), " + ", round(summary(ModelMixedS)$coefficients[2,1],2), "x", sep="") SpeciesResults$InterceptMixed[i] <- summary(ModelMixedS)$coefficients[1,1] SpeciesResults$InterceptMixedSE[i] <- summary(ModelMixedS)$coefficients[1,2] SpeciesResults$SlopeMixed[i] <- summary(ModelMixedS)$coefficients[2,1] SpeciesResults$SlopeMixedSE[i] <- summary(ModelMixedS)$coefficients[2,2] # The proportion of variance explained by the random factor "Individual" is the variance explained by it, divided by the sum of variances explained by it, the other random factor, the residuals SpeciesResults$IndivVar[i] <- VarCorr(ModelMixedS)$lizardID[1,1]/((attr(VarCorr(ModelMixedS), "sc")^2)+VarCorr(ModelMixedS)$lizardID[1,1]) SpeciesResults$DeltaAICc[i] <- AICc(ModelNull)-AICc(ModelMixedS) } rm(tubeSPT, SpeciesT, SP, lizardColumns, temp, ModelNull, ModelMixedS, i) SpeciesResults <- SpeciesResults[order(SpeciesResults$Species),] write.table(SpeciesResults, "output_data/TevsTbCalibrationResults_perSpecies.csv", sep=";") TeCorrections <- data.frame(Species=SpeciesResults$Species, MixedInt= ifelse(SpeciesResults$InterceptMixed-2*SpeciesResults$InterceptMixedSE<0 & SpeciesResults$InterceptMixed+2*SpeciesResults$InterceptMixedSE>0, 0, SpeciesResults$InterceptMixed), MixedSlp= ifelse(SpeciesResults$SlopeMixed-2*SpeciesResults$SlopeMixedSE<1 & SpeciesResults$SlopeMixed+2*SpeciesResults$SlopeMixedSE>1, 1, SpeciesResults$SlopeMixed)) # Because the 95% confidence intervals of the intercept for L. ornatus do not include the 0, then # we have to use also the average coefficient for the slope, instead of 1. # In the case of L. puna and L. orientalis, we do not need to transform the Te data, # but in the case of Liolaemus sp and L. ornatus, we have to use the intercept and slope of their respective model TeCorrections$MixedSlp[2] <- SpeciesResults$SlopeMixed[2] TeCorrections <- TeCorrections[c(3,2, 4, 1),] rownames(TeCorrections) <- NULL write.table(TeCorrections, "output_data/TeCorrections.csv", sep=";") # 3.3. Plots ###### ResultsTable$SVL <- as.numeric(as.character(ResultsTable$SVL)) par(mfrow=c(1,2)) plot(Slope~SVL, data=ResultsTable, type="n", ylim=c(0.5,1.3), ylab="Slopes for heating and cooling", xlab="SVL (mm)") arrows(ResultsTable$SVL, ResultsTable$Slope-ResultsTable$SlopeSE*2, lwd=0.5, ResultsTable$SVL, ResultsTable$Slope+ResultsTable$SlopeSE*2, length=0, angle=90, code=3) points(Slope~SVL, data=ResultsTable, pch=c(21, 22)[ResultsTable$Sex], cex=1, bg=c(rgb(0,1,0,1), rgb(1,.65,0,1), rgb(0,0,1,1), rgb(1,1,0, 1))[ResultsTable$Species]) abline(h=1, lty=2) arrows(SpeciesResults$SVL, SpeciesResults$SlopeMixed-SpeciesResults$SlopeMixedSE*2, SpeciesResults$SVL, SpeciesResults$SlopeMixed+SpeciesResults$SlopeMixedSE*2, length=0.05, angle=90, code=3) points(SlopeMixed~SVL, data=SpeciesResults, pch=23, cex=2, bg=c("green", "orange", "blue", "yellow")[SpeciesResults$Species]) TeSam2012 <- read.table("input_data/TeSam2012.csv", sep=";", header=T) timeseries <- TeSam2012[which(TeSam2012$JULDAY==350), 39] time <- seq(0, 1438, by=2)/60 plot(timeseries[,1]~time, type="l", lwd=2, ylim=c(20, 50), xlim=c(8, 17), ylab="Operative Temperatures (ºC)", xlab="Time of the day") coolheatSeries <- apply(timeseries, 2, function(x) SpeciesResults$InterceptMixed[3] + x*SpeciesResults$SlopeMixed[3]) lines(coolheatSeries[,1]~time, col="blue", lwd=1) coolheatSeries <- apply(timeseries, 2, function(x) TeCorrections$MixedInt[2] + x*TeCorrections$MixedSlp[2]) lines(coolheatSeries[,1]~time, col="orange", lwd=2) coolheatSeries <- apply(timeseries, 2, function(x) TeCorrections$MixedInt[3] + x*TeCorrections$MixedSlp[3]) lines(coolheatSeries[,1]~time, col="yellow", lwd=2) coolheatSeries <- apply(timeseries, 2, function(x) SpeciesResults$InterceptMixed[1] + x*SpeciesResults$SlopeMixed[1]) lines(coolheatSeries[,1]~time, col="green", lwd=1) lines(timeseries[,1]~time) rm(coolheatSeries, time, timeseries, TeSam2012) # 3.4. Function to plot Mixed Models #### plotMixModel <- function(SP, individuals, data){ lizardColumns <- grep(SP, individuals$Species)+1 SpeciesT <- as.numeric(as.character(data[6:172, lizardColumns[1]])) SpeciesT[168:334] <- as.numeric(as.character(data[6:172, lizardColumns[2]])) SpeciesT[335:501] <- as.numeric(as.character(data[6:172, lizardColumns[3]])) tubeSPT <- as.numeric(as.character(data[6:172, lizardColumns[1]+1])) tubeSPT[168:334] <- as.numeric(as.character(data[6:172, lizardColumns[2]+1])) tubeSPT[335:501] <- as.numeric(as.character(data[6:172, lizardColumns[3]+1])) if(length(lizardColumns)>3){ SpeciesT[502:668] <- as.numeric(as.character(data[6:172, lizardColumns[4]])) tubeSPT[502:668] <- as.numeric(as.character(data[6:172, lizardColumns[4]+1])) } if(length(lizardColumns)>4){ SpeciesT[669:835] <- as.numeric(as.character(data[6:172, lizardColumns[5]])) tubeSPT[669:835] <- as.numeric(as.character(data[6:172, lizardColumns[5]+1])) } if(length(lizardColumns)>5){ SpeciesT[836:1002] <- as.numeric(as.character(data[6:172, lizardColumns[6]])) tubeSPT[836:1002] <- as.numeric(as.character(data[6:172, lizardColumns[6]+1])) } temp <- data.frame(lizardT=SpeciesT, tubeT=tubeSPT, lizardID=c(rep(1, 167), rep(2, 167), rep(3, 167), rep(4,167), rep(5, 167), rep(6,167))[1:length(SpeciesT)], SVL=c(rep(individuals$SVL[lizardColumns[1]-1], 167), rep(individuals$SVL[lizardColumns[2]-1], 167), rep(individuals$SVL[lizardColumns[3]-1], 167), rep(individuals$SVL[lizardColumns[4]-1],167), rep(individuals$SVL[lizardColumns[5]-1], 167), rep(individuals$SVL[lizardColumns[6]-1],167))[1:length(SpeciesT)]) ModelMixedS <- lmer(lizardT ~ tubeT + (1+tubeT|lizardID), data=temp) # both random intercept and random slope colores <- c(rgb(0,0,1,.5), rgb(0,1,0,.5), rgb(0,1,1,.5), rgb(0.6,.12,.8,.5), rgb(1,.65,0,.5), rgb(1,1,1,.5)) plot(lizardT ~tubeT, data=temp, col=colores[temp$lizardID], pch=19, cex=.5, xlim=c(10,40), ylim=c(10,40), xlab="Te from physical model (ºC)", ylab="T from inmobilized lizard (ºC)") abline(a=summary(ModelMixedS)$coefficients[1,1], b=summary(ModelMixedS)$coefficients[2,1]) #abline(a=summary(model)$coefficients[1,1], b=summary(model)$coefficients[2,1], col="red") newx <- data.frame(tubeT=seq(min(temp$tubeT, na.rm=T), max(temp$tubeT, na.rm=T), length.out=100)) preds <- predict(ModelMixedS, newdata = newx, re.form=~0) ci_line<-bootMer(ModelMixedS,FUN=function(.) predict(.,newdata=newx,re.form=~0),nsim=200) # https://www.r-bloggers.com/linear-mixed-effect-model-workflow/ ci_regT<-apply(ci_line$t,2,function(x) x[order(x)][c(5,195)]) lines(newx$tubeT,ci_regT[1,],lty=2) lines(newx$tubeT,ci_regT[2,],lty=2) abline(a=0, b=1, lty=2, col="red") VarF <- var(as.vector(fixef(ModelMixedS) %*% t(ModelMixedS@pp$X))) # R2 <- (VarF + VarCorr(ModelMixedS)$lizardID[1])/(VarF + VarCorr(ModelMixedS)$lizardID[1] + (attr(VarCorr(ModelMixedS), "sc")^2)) # text(15, 40, paste("R2GLMM =", round(R2, 3))) IndivVar <- VarCorr(ModelMixedS)$lizardID[1,1]/((attr(VarCorr(ModelMixedS), "sc")^2)+VarCorr(ModelMixedS)$lizardID[1,1]) text(15, 38, paste("Expl. var. random f. =", round(IndivVar, 3))) ModelNull <- lmer(lizardT ~ 1 + (1+tubeT|lizardID), data=temp) # both random intercept and random slope text(20, 36, paste("∆AICc", "=", round(AICc(ModelNull)-AICc(ModelMixedS), 1))) text(20, 34, paste ("y = ", round(summary(ModelMixedS)$coefficients[1,1],2), " + ", round(summary(ModelMixedS)$coefficients[2,1],2), "x", sep="")) } par(mfrow=c(2,2)) plotMixModel(SP="pun", individuals = individuals, data=data) plotMixModel(SP="orn", individuals = individuals, data=data) plotMixModel(SP="sp", individuals = individuals, data=data) plotMixModel(SP="ori", individuals = individuals, data=data)