#################################################################################### # Code associated with paper by De Frenne et al. # # Global buffering of temperatures under forest canopies # #################################################################################### setwd("D:/Users/pdfrenne/Documents/R/Microclimate-meta") #set working directory library(nlme) library(lme4) library(boot) library(MCMCglmm) library(dplyr) library(CoordinateCleaner) library(MuMIn) library(ggplot2) mean.na <- function(x) mean(x, na.rm=T) max.na <- function(x) max(x, na.rm=T) se.na <- function(x) sd(x, na.rm=T)/sqrt(length(x)) mean.w <- function(x, w) sum(x*w) ################################################################## ################## Reading the data ########################### ################################################################## Data <- read.csv("Microclimate-database.csv",header=T) glimpse(Data) head(Data) #check file Data$Latitude<-as.numeric(as.character(Data$Latitude)) #make numeric Data$Longitude<-as.numeric(as.character(Data$Longitude)) Data$Elevation<-as.numeric(as.character(Data$Elevation)) Data$Buffering<-as.numeric(as.character(Data$Buffering)) Data$Vegetation.type<-as.factor(Data$Vegetation.type) #make factor Data$Height.sensor<-as.numeric(as.character(Data$Height.sensor)) Data$Height.sensor<-as.numeric(as.character(Data$Height.sensor)) Data$Temp.forest<-as.numeric(as.character(Data$Temp.forest)) Data$Temp.open<-as.numeric(as.character(Data$Temp.open)) Data$Season<-factor(Data$Season) Data$Number.days<-as.numeric(as.character(Data$Number.days)) #Add variable which is the absolute value of the latitude absLatitude<-abs(Data$Latitude) Data$absLatitude<-absLatitude #BIOME: Add variable as category of the latitude: tropical: 0-23.5, temperate: 23.5-55, boreal: >55 Lat.category<-numeric(length=nrow(Data)) for (i in 1:length(Lat.category)){ if(absLatitude[i] <23.5) Lat.category[i] <- 'Tropical' if((absLatitude[i] >23.5) & (absLatitude[i] <55)) Lat.category[i] <- 'Temperate' if(absLatitude[i] >55) Lat.category[i] <- 'Boreal' } Data$Lat.category<-as.factor(Lat.category) #Split in different datasets for Tmean, Tmin and Tmax dta.Mean<-subset(Data,Variable.statistic=="Tmean") dta.Max<-subset(Data,Variable.statistic=="Tmax") dta.Min<-subset(Data,Variable.statistic=="Tmin") ####################################################################################### ################## Descriptive stats of the data ########################## ####################################################################################### #[Supplementary Table 1] summary(dta.Max$Buffering) summary(dta.Mean$Buffering) summary(dta.Min$Buffering) summary(Data$Buffering) length(unique(Data$absLatitude)) length(unique(Data$Longitude)) unique(Data$First.author) quantile(dta.Max$Buffering, probs = c(0.05,0.95), na.rm = TRUE) quantile(dta.Mean$Buffering, probs = c(0.05,0.95), na.rm = TRUE) quantile(dta.Min$Buffering, probs = c(0.05,0.95), na.rm = TRUE) quantile(Data$Buffering, probs = c(0.05,0.95), na.rm = TRUE) ###########Linear mixed-effect models using lmer and study as random factor summary(lmer(Buffering~1+(1|First.author),data=dta.Max,na.action=na.omit,REML=T)) summary(lmer(Buffering~1+(1|First.author),data=dta.Mean,na.action=na.omit,REML=T)) summary(lmer(Buffering~1+(1|First.author),data=dta.Min,na.action=na.omit,REML=T)) summary(lmer(Buffering~1+(1|First.author),data=Data,na.action=na.omit,REML=T)) #Check significance of intercept M0<-lmer(Buffering~1+(1|First.author),data=Data,na.action=na.omit,REML=F) M1<-lmer(Buffering~-1+(1|First.author),data=Data,na.action=na.omit,REML=F) anova(M0,M1) M0<-lmer(Buffering~1+(1|First.author),data=dta.Max,na.action=na.omit,REML=F) M1<-lmer(Buffering~-1+(1|First.author),data=dta.Max,na.action=na.omit,REML=F) anova(M0,M1) M0<-lmer(Buffering~1+(1|First.author),data=dta.Mean,na.action=na.omit,REML=F) M1<-lmer(Buffering~-1+(1|First.author),data=dta.Mean,na.action=na.omit,REML=F) anova(M0,M1) M0<-lmer(Buffering~1+(1|First.author),data=dta.Min,na.action=na.omit,REML=F) M1<-lmer(Buffering~-1+(1|First.author),data=dta.Min,na.action=na.omit,REML=F) anova(M0,M1) ####################################################################################### ############ Buffering ~ predictors ########################## ####################################################################################### #[Supplementary Table 5 and 2] vars<-c(24,26,27,9,10,19,20,21,17,15,16) #List of which variables to include in the models results<-matrix(nrow = length(vars), ncol = 6) colnames(results)<-c('var','estimate','Chisq','p','ntotal','nstudies') for (i in 1:length(vars)){ #Which data to use? #change on this line here to dta.Min, dta.Mean, dta.Max or Data intermediateData<-Data intermediateData<-intermediateData[!is.na(intermediateData[,vars[i]]),] predictor <- intermediateData[,vars[i]] results[i,1]<-colnames(intermediateData)[vars[i]] M2<-lmer(Buffering~predictor + (1|First.author),data=intermediateData,na.action=na.omit,REML=T) M1<-lmer(Buffering~predictor + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) M0<-lmer(Buffering~1 + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) results[i,2]<-fixef(M2)[2] results[i,3]<-anova(M0,M1)$Chisq[2] results[i,4]<-anova(M0,M1)$Pr[2] results[i,5]<-sum(!is.na(intermediateData$Buffering)) results[i,6]<-nrow(coef(M1)$First.author) } results<-data.frame(results) write.table(results,"clipboard",sep="\t", col.names=NA) ####################################################################################### ################## Next analyses of the data ########################## ####### Temp.forest~Temp.open ########################## ####################################################################################### #[Supplementary Table 3] results<-matrix(nrow = 4, ncol = 9) colnames(results)<-c('var','intercept','estimate','Chisq','p','ntotal','nstudies','R2m','R2c') results[,1]<-c('Max','Mean','Min','All') #Then repeat this 4 times, once for each dataset #change on this line here to dta.Min, dta.Mean, dta.Max or Data intermediateData<-Data i=4 #1 = dta.Max, 2=dta.Mean, 3=dta.Min, 4=Data intermediateData<-intermediateData[!is.na(intermediateData[,24]),] #omit all NAs in Temp.open H2<-lmer(Temp.forest~Temp.open + (1|First.author),data=intermediateData,na.action=na.omit,REML=T) H1<-lmer(Temp.forest~Temp.open + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) H0<-lmer(Temp.forest~1 + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) # summary(H2) # anova(H1,H0) results[i,2]<-fixef(H2)[1] results[i,3]<-fixef(H2)[2] results[i,4]<-anova(H0,H1)$Chisq[2] results[i,5]<-anova(H0,H1)$Pr[2] results[i,6]<-sum(!is.na(intermediateData$Temp.open)) results[i,7]<-nrow(coef(H1)$First.author) results[i,8]<-r.squaredGLMM(H2)[1] results[i,9]<-r.squaredGLMM(H2)[2] #And then write away results<-data.frame(results) write.table(results,"clipboard",sep="\t", col.names=NA) ####################################################################################### ################### Interactions testing ################### ####################################################################################### #[Supplementary Table 6] vars<-c(26,27,9,10,19,20,21,17,15,16) #List of which variables to include in the models results<-matrix(nrow = length(vars), ncol = 9) colnames(results)<-c('var','estimateTemp.open','estimate-predictor','Chisq-predictor','p-predictor','Chisq-interact','p-interact','ntotal','nstudies') for (i in 1:length(vars)){ #Which data to use? #change on this line here to dta.Min, dta.Mean, dta.Max or Data intermediateData<-dta.Min intermediateData<-intermediateData[!is.na(intermediateData[,vars[i]]),] predictor <- intermediateData[,vars[i]] results[i,1]<-colnames(intermediateData)[vars[i]] M3<-lmer(Buffering~Temp.open * predictor + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) M2<-lmer(Buffering~Temp.open + predictor + (1|First.author),data=intermediateData,na.action=na.omit,REML=T) M1<-lmer(Buffering~Temp.open + predictor + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) M0<-lmer(Buffering~Temp.open + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) results[i,2]<-fixef(M3)[2] results[i,3]<-fixef(M3)[3] results[i,4]<-anova(M0,M1)$Chisq[2] results[i,5]<-anova(M0,M1)$Pr[2] results[i,6]<-anova(M3,M1)$Chisq[2] results[i,7]<-anova(M3,M1)$Pr[2] results[i,8]<-sum(!is.na(intermediateData$Temp.open)) results[i,9]<-nrow(coef(M1)$First.author) } results<-data.frame(results) write.table(results,"clipboard",sep="\t", col.names=NA) ####################################################################################### ################## Sensitivity analysis of the data ########################## ####### Excluding buffering of <-20°C ########################## ####################################################################################### #[Supplementary Table 8] results<-matrix(nrow = 4, ncol = 9) colnames(results)<-c('var','intercept','estimate','Chisq','p','ntotal','nstudies','R2m','R2c') results[,1]<-c('Max','Mean','Min','All') #Then repeat this 4 times, once for each dataset #change on this line here to dta.Min, dta.Mean, dta.Max or Data intermediateData<-dta.Max i=1 #1 = dta.Max, 2=dta.Mean, 3=dta.Min, 4=Data #omit buffering values smaller than -20°C #This is the 0.25% quantile quantile(Data$Buffering, probs = c(0.0025), na.rm = TRUE) intermediateData<-subset(intermediateData,Buffering>-20) intermediateData<-intermediateData[!is.na(intermediateData[,24]),] #omit all NAs in Temp.open H2<-lmer(Buffering~Temp.open + (1|First.author),data=intermediateData,na.action=na.omit,REML=T) H1<-lmer(Buffering~Temp.open + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) H0<-lmer(Buffering~1 + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) # summary(H2) # anova(H1,H0) results[i,2]<-fixef(H2)[1] results[i,3]<-fixef(H2)[2] results[i,4]<-anova(H0,H1)$Chisq[2] results[i,5]<-anova(H0,H1)$Pr[2] results[i,6]<-sum(!is.na(intermediateData$Temp.open)) results[i,7]<-nrow(coef(H1)$First.author) results[i,8]<-r.squaredGLMM(H2)[1] results[i,9]<-r.squaredGLMM(H2)[2] #And then write away results<-data.frame(results) write.table(results,"clipboard",sep="\t", col.names=NA) ####################################################################################### ############ Random slopes ########################## ####################################################################################### #[Supplementary Table 7] intData<-dta.Min M0<-lmer(Buffering~Temp.open +(1|First.author),data=intData,na.action=na.omit,REML=F) M1<-lmer(Buffering~Temp.open +(1+Temp.open|First.author),data=intData,na.action=na.omit,REML=F) M2<-lmer(Buffering~Temp.open +(1+Temp.open|First.author),data=intData,na.action=na.omit,REML=T) anova(M0,M1) fixef(M2)[2] r.squaredGLMM(M2) # summary(M2) ####################################################################################### ################## Supplementary: GAMMs ########################## ################# Testing for non-linearity ########################## ####################################################################################### library(mgcv) # [Supplementary Fig. 4] par(mfrow=c(2,2)) intData<-Data temprange<-seq(min(intData$Temp.open,na.rm=T),max(intData$Temp.open,na.rm=T),length.out=500) GM <- gamm(Buffering ~ s(Temp.open), random=list(First.author=~1),data=intData,na.action=na.omit) predvals <- predict(GM, list(Temp.open=temprange), type="response",se.fit=T) predvals<-as.data.frame(predvals) plot(intData$Buffering ~ intData$Temp.open,xlab="Temperature outside forest (°C)",ylab="Temperature offset (°C)",main="All data") MM<-lmer(Buffering ~ Temp.open +(1|First.author),data=intData,na.action=na.omit,REML=T) abline(a=fixef(MM)[1],b=fixef(MM)[2],col='blue', lwd=2) lines(temprange, predvals$fit, col="black",lwd=2) lines(temprange, predvals$fit-predvals$se.fit, lwd=2,col="grey",lty=2) lines(temprange, predvals$fit+predvals$se.fit, lwd=2,col="grey",lty=2) polygon(c(rev(temprange),temprange), c(sort(as.numeric(predvals$fit+predvals$se.fit)),rev(sort(as.numeric(predvals$fit-predvals$se.fit)))), col = "#8B000020", border = NA) abline(h=0,col='red', lwd=2, lty="dashed") mtext("a", side=3, line=1, adj=0, font=2, cex=2) intData<-dta.Max temprange<-seq(min(intData$Temp.open,na.rm=T),max(intData$Temp.open,na.rm=T),length.out=500) GM <- gamm(Buffering ~ s(Temp.open), random=list(First.author=~1),data=intData,na.action=na.omit) predvals <- predict(GM, list(Temp.open=temprange), type="response",se.fit=T) predvals<-as.data.frame(predvals) plot(intData$Buffering ~ intData$Temp.open,xlab="Temperature outside forest (°C)",ylab="Temperature offset (°C)",main="Max. temperatures") MM<-lmer(Buffering ~ Temp.open +(1|First.author),data=intData,na.action=na.omit,REML=T) abline(a=fixef(MM)[1],b=fixef(MM)[2],col='blue', lwd=2) lines(temprange, predvals$fit, col="black",lwd=2) lines(temprange, predvals$fit-predvals$se.fit, lwd=2,col="grey",lty=2) lines(temprange, predvals$fit+predvals$se.fit, lwd=2,col="grey",lty=2) polygon(c(rev(temprange),temprange), c(sort(as.numeric(predvals$fit+predvals$se.fit)),rev(sort(as.numeric(predvals$fit-predvals$se.fit)))), col = "#8B000020", border = NA) abline(h=0,col='red', lwd=2, lty="dashed") mtext("b", side=3, line=1, adj=0, font=2, cex=2) intData<-dta.Mean temprange<-seq(min(intData$Temp.open,na.rm=T),max(intData$Temp.open,na.rm=T),length.out=500) GM <- gamm(Buffering ~ s(Temp.open), random=list(First.author=~1),data=intData,na.action=na.omit) predvals <- predict(GM, list(Temp.open=temprange), type="response",se.fit=T) predvals<-as.data.frame(predvals) plot(intData$Buffering ~ intData$Temp.open,xlab="Temperature outside forest (°C)",ylab="Temperature offset (°C)",main="Mean temperatures") MM<-lmer(Buffering ~ Temp.open +(1|First.author),data=intData,na.action=na.omit,REML=T) abline(a=fixef(MM)[1],b=fixef(MM)[2],col='blue', lwd=2) lines(temprange, predvals$fit, col="black",lwd=2) lines(temprange, predvals$fit-predvals$se.fit, lwd=2,col="grey",lty=2) lines(temprange, predvals$fit+predvals$se.fit, lwd=2,col="grey",lty=2) polygon(c(rev(temprange),temprange), c(sort(as.numeric(predvals$fit+predvals$se.fit)),rev(sort(as.numeric(predvals$fit-predvals$se.fit)))), col = "#8B000020", border = NA) abline(h=0,col='red', lwd=2, lty="dashed") mtext("c", side=3, line=1, adj=0, font=2, cex=2) intData<-dta.Min temprange<-seq(min(intData$Temp.open,na.rm=T),max(intData$Temp.open,na.rm=T),length.out=500) GM <- gamm(Buffering ~ s(Temp.open), random=list(First.author=~1),data=intData,na.action=na.omit) predvals <- predict(GM, list(Temp.open=temprange), type="response",se.fit=T) predvals<-as.data.frame(predvals) plot(intData$Buffering ~ intData$Temp.open,xlab="Temperature outside forest (°C)",ylab="Temperature offset (°C)",main="Min. temperatures") MM<-lmer(Buffering ~ Temp.open +(1|First.author),data=intData,na.action=na.omit,REML=T) abline(a=fixef(MM)[1],b=fixef(MM)[2],col='blue', lwd=2) lines(temprange, predvals$fit, col="black",lwd=2) lines(temprange, predvals$fit-predvals$se.fit, lwd=2,col="grey",lty=2) lines(temprange, predvals$fit+predvals$se.fit, lwd=2,col="grey",lty=2) polygon(c(rev(temprange),temprange), c(sort(as.numeric(predvals$fit+predvals$se.fit)),rev(sort(as.numeric(predvals$fit-predvals$se.fit)))), col = "#8B000020", border = NA) abline(h=0,col='red', lwd=2, lty="dashed") mtext("d", side=3, line=1, adj=0, font=2, cex=2) ####################################################################################### ################## Supplementary Info PLOTS ########################## ####################################################################################### #graph with histograms [Supplementary Fig. 1] par(mfrow=c(2,2)) hist(Data$Elevation, main="",xlab="Elevation (m a.s.l.)") mtext("a", side=3, line=1, adj=0, font=2, cex=2) hist(Data$Latitude, main="",xlab="Latitude (°)") mtext("b", side=3, line=1, adj=0, font=2, cex=2) plot(Data$Vegetation.type, main="",xlab="Vegetation type",xaxt="n",ylab="Frequency") axis(1, at = c(0.8,2,3.2), labels = c("Deciduous", "Evergreen","Mixed"),cex.axis=1, tick=F) mtext("c", side=3, line=1, adj=0, font=2, cex=2) plot(Data$Season,xlab="Season",ylab="Frequency") mtext("d", side=3, line=1, adj=0, font=2, cex=2) #graph per study and country [Supplementary Fig. 2] par(mfrow=c(2,1)) plot(Data$Buffering~as.factor((Data$First.author)),ylab="Temperature offset (°C)",xlab="Study number") mtext("a", side=3, line=1, adj=0, font=2, cex=2) abline(h=0,col='red', lwd=2, lty="dashed") plot(Data$Buffering~Data$Country.code,ylab="Temperature offset (°C)",xlab="Country (ISO code)") abline(h=0,col='red', lwd=2, lty="dashed") mtext("b", side=3, line=1, adj=0, font=2, cex=2) #plots with Vegetation type [Supplementary Fig. 6] par(mfrow=c(2,2)) plot(Data$Buffering~Data$Vegetation.type,xlab="",ylab="Temperature offset (°C)", main="All data",xaxt="n") axis(1, at = c(1,2,3), labels = c("Deciduous", "Evergreen","Mixed"),las=1,cex.axis=1) abline(h=0,col='red', lwd=2, lty="dashed") mtext("a", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Max$Buffering~dta.Max$Vegetation.type,xlab="",ylab="Temperature offset (°C)", main="Max. temperature difference",xaxt="n") axis(1, at = c(1,2,3), labels = c("Deciduous", "Evergreen","Mixed"),las=1,cex.axis=1) abline(h=0,col='red', lwd=2, lty="dashed") mtext("b", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Mean$Buffering~dta.Mean$Vegetation.type,xlab="",ylab="Temperature offset (°C)", main="Mean temperature difference",xaxt="n") axis(1, at = c(1,2,3), labels = c("Deciduous", "Evergreen","Mixed"),las=1,cex.axis=1) abline(h=0,col='red', lwd=2, lty="dashed") mtext("c", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Min$Buffering~dta.Min$Vegetation.type,xlab="",ylab="Temperature offset (°C)", main="Min temperature difference",xaxt="n") axis(1, at = c(1,2,3), labels = c("Deciduous", "Evergreen","Mixed"),las=1,cex.axis=1) abline(h=0,col='red', lwd=2, lty="dashed") mtext("d", side=3, line=1, adj=0, font=2, cex=2) #relationship with sensor height [Supplementary Fig. 7] par(mfrow=c(2,2)) plot(Data$Buffering~(Data$Height.sensor),xlab="Height sensor (m)",ylab="Temperature offset (°C)", main="All data") abline(h=0,col='red', lwd=2, lty="dashed") abline(a=-1.8796,b=0.3965,col='blue', lwd=2) abline(a=-2.0729,b=0.7620,col='grey', lwd=2) mtext("a", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Max$Buffering~(dta.Max$Height.sensor),xlab="Height sensor (m)",ylab="Temperature offset (°C)", main="Max temperature difference") abline(h=0,col='red', lwd=2, lty="dashed") abline(a=-4.5775,b=0.6637,col='blue', lwd=2) abline(a=-5.6885,b=2.2725,col='grey', lwd=2) mtext("b", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Mean$Buffering~(dta.Mean$Height.sensor),xlab="Height sensor (m)",ylab="Temperature offset (°C)", main="Mean temperature difference") abline(h=0,col='red', lwd=2, lty="dashed") abline(a=-1.9054,b=0.5374,col='blue', lwd=2) mtext("c", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Min$Buffering~(dta.Min$Height.sensor),xlab="Height sensor (m)",ylab="Temperature offset (°C)", main="Min temperature difference") abline(h=0,col='red', lwd=2, lty="dashed") abline(a=1.4803,b=-0.4613,col='blue', lwd=2) mtext("d", side=3, line=1, adj=0, font=2, cex=2) #intercepts and slopes from mixed model temp<-dta.Max[!is.na(dta.Max[,16]),] M1<-lmer(Buffering~Height.sensor + (1|First.author),data=temp,na.action=na.omit,REML=T) summary(M1) #And without outliers with sensor height > 3 m temp<-dta.Max[!is.na(dta.Max[,16]),] temp<-subset(temp,Height.sensor<3) M1<-lmer(Buffering~Height.sensor + (1|First.author),data=temp,na.action=na.omit,REML=T) summary(M1) #Plots with Temp.open vs Temp.forest: [Supplementary Fig. 3] par(mfrow=c(2,2)) plot(Data$Temp.forest~Data$Temp.open,xlab="Temperature outside forest (°C)",ylab="Forest temperature (°C)", main="All data", xlim=c(-10,50),ylim=c(-10,50)) abline(a=0,b=1,col='red', lwd=2, lty="dashed") abline(a=2.823,b=0.753,col='blue', lwd=2) mtext("a", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Max$Temp.forest~dta.Max$Temp.open,xlab="Temperature outside forest (°C)",ylab="Forest temperature (°C)", main="Max temperature difference", xlim=c(0,50),ylim=c(0,50)) abline(a=0,b=1,col='red', lwd=2, lty="dashed") abline(a=4.190,b=0.681,col='blue', lwd=2) mtext("b", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Mean$Temp.forest~dta.Mean$Temp.open,xlab="Temperature outside forest (°C)",ylab="Forest temperature (°C)", main="Mean temperature difference", xlim=c(-5,35),ylim=c(-5,35)) abline(a=0,b=1,col='red', lwd=2, lty="dashed") abline(a=1.537,b=0.816,col='blue', lwd=2) mtext("c", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Min$Temp.forest~dta.Min$Temp.open,xlab="Temperature outside forest (°C)",ylab="Forest temperature (°C)", main="Min temperature difference", xlim=c(-10,30),ylim=c(-10,30)) abline(a=0,b=1,col='red', lwd=2, lty="dashed") abline(a=1.772,b=0.927,col='blue', lwd=2) mtext("d", side=3, line=1, adj=0, font=2, cex=2) #intercepts and slopes from mixed model temp<-dta.Max[!is.na(dta.Max[,24]),] M1<-lmer(Temp.forest~Temp.open + (1|First.author),data=temp,na.action=na.omit,REML=T) summary(M1) ####################################################################################### #### Effects of temperature anomalies ######## #### Analysis with the Worldclim V2 dataset ######## ####################################################################################### #[Supplementary Fig. 5 and Table 4] ##### Create coordinates crs.geo <- CRS("+proj=longlat +ellps=WGS84") # geographical, datum WGS84 Coord <- subset(Data, select=c("Country.code", "Latitude", "Longitude")) Long <- Coord$Longitude[!is.na(Coord$Longitude)] Lat <- Coord$Latitude[!is.na(Coord$Latitude)] coordnew<-data.frame('Long2' = c(Long), 'Lat2' = c(Lat)) coordinates(coordnew) <- c("Long2", "Lat2") # set spatial coordinates proj4string(coordnew) <- crs.geo # define projection system of our data ####################################################### ######### Preamble ########### ####################################################### library(sp) # classes for spatial data library(raster) # grids, rasters library(maptools) #includes readShapePoly function library(rgdal) setwd("E:/Users/pdfrenne/WorldClim") ############# Reading of the Worldclim data #tavg, 30 sec resolution, WorldClim v2 # Fick, S.E. and R.J. Hijmans, 2017. Worldclim 2: New 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology. tmean1 <- raster("wc2.0_30s_tavg_01.tif") # Tmin for January t1vals<-extract(tmean1,coordnew) tmean2 <- raster("wc2.0_30s_tavg_02.tif") # Tmin for Feb t2vals<-extract(tmean2,coordnew) tmean3 <- raster("wc2.0_30s_tavg_03.tif") # Tmin for March t3vals<-extract(tmean3,coordnew) tmean4 <- raster("wc2.0_30s_tavg_04.tif") # Tmin for Apr t4vals<-extract(tmean4,coordnew) tmean5 <- raster("wc2.0_30s_tavg_05.tif") # Tmin for may t5vals<-extract(tmean5,coordnew) tmean6 <- raster("wc2.0_30s_tavg_06.tif") # Tmin for jun t6vals<-extract(tmean6,coordnew) tmean7 <- raster("wc2.0_30s_tavg_07.tif") # Tmin for jul t7vals<-extract(tmean7,coordnew) tmean8 <- raster("wc2.0_30s_tavg_08.tif") # Tmin for aug t8vals<-extract(tmean8,coordnew) tmean9 <- raster("wc2.0_30s_tavg_09.tif") # Tmin for sep t9vals<-extract(tmean9,coordnew) tmean10 <- raster("wc2.0_30s_tavg_10.tif") # Tmin for oct t10vals<-extract(tmean10,coordnew) tmean11 <- raster("wc2.0_30s_tavg_11.tif") # Tmin for noc t11vals<-extract(tmean11,coordnew) tmean12 <- raster("wc2.0_30s_tavg_12.tif") # Tmin for Dec t12vals<-extract(tmean12,coordnew) #And then add temperature values (tXXvals) to the dataframe "Data" Data$t1vals<-t1vals Data$t2vals<-t2vals Data$t3vals<-t3vals Data$t4vals<-t4vals Data$t5vals<-t5vals Data$t6vals<-t6vals Data$t7vals<-t7vals Data$t8vals<-t8vals Data$t9vals<-t9vals Data$t10vals<-t10vals Data$t11vals<-t11vals Data$t12vals<-t12vals tmean<-(t1vals+t2vals+t3vals+t4vals+t5vals+t6vals+t7vals+t8vals+t9vals+t10vals+t11vals+t12vals)/12 Data$tmean<-tmean #calculate anomaly diff<-Data$Temp.open-Data$tmean Data$diff<-diff str(Data) #Split in different datasets for Tmean, Tmin and Tmax dta.Mean<-subset(Data,Variable.statistic=="Tmean") dta.Max<-subset(Data,Variable.statistic=="Tmax") dta.Min<-subset(Data,Variable.statistic=="Tmin") ####################################################################################### ########## and now the actual analyses and plots ########### ###################################################################################### par(mfrow=c(2,2)) #Supplementary Fig. 5 plot(Data$Buffering~Data$diff,xlab="Macroclimate temperature anomaly (°C)",ylab="Temperature offset (°C)", main="All data") abline(h=0,col='red', lwd=2, lty="dashed") abline(a=-1.11530,b=-0.25308,col='blue', lwd=2) mtext("a", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Max$Buffering~dta.Max$diff,xlab="Macroclimate temperature anomaly (°C)",ylab="Temperature offset (°C)", main="Max. temperatures") abline(h=0,col='red', lwd=2, lty="dashed") abline(a=-0.88338,b=-0.32906,col='blue', lwd=2) mtext("b", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Mean$Buffering~dta.Mean$diff,xlab="Macroclimate temperature anomaly (°C)",ylab="Temperature offset (°C)", main="Mean temperatures") abline(h=0,col='red', lwd=2, lty="dashed") abline(a=-1.17082,b=-0.18357,col='blue', lwd=2) mtext("c", side=3, line=1, adj=0, font=2, cex=2) plot(dta.Min$Buffering~dta.Min$diff,xlab="Macroclimate temperature anomaly (°C)",ylab="Temperature offset (°C)", main="Min temperatures") abline(h=0,col='red', lwd=2, lty="dashed") abline(a=0.62435,b=-0.07519,col='blue', lwd=2) mtext("d", side=3, line=1, adj=0, font=2, cex=2) #mixed effect model #Supplementary Table 4 results<-matrix(nrow = 4, ncol = 9) colnames(results)<-c('var','intercept','estimate','Chisq','p','ntotal','nstudies','R2m','R2c') results[,1]<-c('Max','Mean','Min','All') #Then repeat this 4 times, once for each dataset #change on this line here to dta.Min, dta.Mean, dta.Max or Data i=2 #1 = dta.Max, 2=dta.Mean, 3=dta.Min, 4=Data intermediateData<-dta.Mean intermediateData<-intermediateData[!is.na(intermediateData[,41]),] #omit all NAs in diff H2<-lmer(Buffering~diff + (1|First.author),data=intermediateData,na.action=na.omit,REML=T) H1<-lmer(Buffering~diff + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) H0<-lmer(Buffering~1 + (1|First.author),data=intermediateData,na.action=na.omit,REML=F) # summary(H2) # anova(H1,H0) results[i,2]<-fixef(H2)[1] results[i,3]<-fixef(H2)[2] results[i,4]<-anova(H0,H1)$Chisq[2] results[i,5]<-anova(H0,H1)$Pr[2] results[i,6]<-sum(!is.na(intermediateData$Temp.open)) results[i,7]<-nrow(coef(H1)$First.author) results[i,8]<-r.squaredGLMM(H2)[1] results[i,9]<-r.squaredGLMM(H2)[2] #And then write away results<-data.frame(results) write.table(results,"clipboard",sep="\t", col.names=NA) ###################################################################################### # THIS SCRIPTS EXTRACT VALUES FROM VARIABLES CAPTURING TOPOGRAPHIC HETEROGENEITY # Jonathan Lenoir (29.11.2018) ###################################################################################### library(raster) library(rgdal) # # Path to files # root <- "C:/Users/Admin/Documents/xxx/" # # Import the database and make a spatial dataframe for spatial extraction # id_cords <- read.csv(paste(root, "Microclimate-database.csv", sep=""), header=TRUE, sep=",")[, c(1,6,7,9)] coordinates(id_cords) <- c('Longitude','Latitude') crs(id_cords) <- "+proj=longlat +ellps=WGS84 +datum=WGS84" # # Load raster layers derived from Global Multi-resolution Terrain Elevation Data 2010 (GMTED) at 250 m resolution # See this article from Amatulli et al. (2018): https://www.nature.com/articles/sdata201840 # Standard deviation of elevation values aggregated at 1 km resolution ele_sd_1km <- raster(paste(root, "elevation_1KMsd_GMTEDsd.tif", sep="")) # Median of topographic position index (TPI) values at 1 km resolution # TPI is the difference between the elevation of a focal cell and the mean of its 8 surrounding cells # Positive and negative values correspond to ridges and valleys, respectively, while zero values correspond to flat areas tpi_md_1km <- raster(paste(root, "tpi_1KMmd_GMTEDmd.tif", sep="")) # Prepare the database for spatial extraction # id_cords@data[, "RASTER_CELL_ID"] <- cellFromXY(ele_sd_1km, coordinates(id_cords)) id_cords@data[, "X_CELL"] <- xFromCell(ele_sd_1km, id_cords$RASTER_CELL_ID) id_cords@data[, "Y_CELL"] <- yFromCell(ele_sd_1km, id_cords$RASTER_CELL_ID) # # Extract values from raster layers # id_cords@data[, "ele_sd_1km"] <- extract(ele_sd_1km, id_cords$RASTER_CELL_ID) id_cords@data[, "tpi_md_1km"] <- extract(tpi_md_1km, id_cords$RASTER_CELL_ID) # str(id_cords@data) # # str(id_cords@data) # # Save data into a .txt file # getwd() setwd(root) write.csv(id_cords@data, file="topographic_heterogeneity.csv") # q()