############################################################################################################################### ## Underprediction of extirpation and colonisation following climate and land-use change using species distribution models ##
###############################################################################################################################

# Alistair Auffret, March 2024

#### Bring in data
## Species data
spe.df<-read.csv("SpeciesData.csv") # species-level prediction accuracy and traits

## Grid-square data
sq.df<-read.csv("GridSquareData.csv", stringsAsFactors = FALSE) # grid-square level prediction accuracy and environmental change


mod.list<-list() # make list for all the tables
#############################################
############ Model evaluation ###############
#############################################

### TSS: overall performance using different evaluation methods

# Traditional cross validation
mean(spe.df$TSS_cross) # mean 0.65
range(spe.df$TSS_cross) # range 0.42, 0.87

# Independent temporal validation
# Data-derived thresholds for probability of occurrence being counted as presence in SDM predictions.
mean(spe.df$TSS_ind) # mean 0.28
range(spe.df$TSS_ind) # range 0, 0.72

# Relationship between validation methods. Output corresponds to Table S1.
m.tss.corr<-glm(TSS_ind~TSS_cross, data= spe.df)
mod.list[["tss.corr"]]<-data.frame(mod="tss.corr",pred=rownames(summary(m.tss.corr)$coefficients), summary(m.tss.corr)$coefficients, r2=with(summary(m.tss.corr), 1 - deviance/null.deviance))

# Sensitivity analyses with lower (data-derived minus 0.1 on 0-1 scale) thresholds for probability of occurrence being counted as presence in SDM predictions
mean(spe.df$TSS_ind.upr) # mean 0.29
range(spe.df$TSS_ind.upr) # range 0, 0.75

# Sensitivity analyses with upper (data-derived plus 0.1) thresholds for probability of occurrence being counted as presence in SDM predictions
mean(spe.df$TSS_ind.lwr) # mean 0.27
range(spe.df$TSS_ind.lwr) # range 0, 0.68

# AUC validation (temporal)
mean(spe.df$AUC)
range(spe.df$AUC)


#############################################
############ Species models #################
#############################################

### Species-level prediction of colonisations, extirpations and persistences using independent validation
# Models have predictive accuracy as a response, and species habitat and climate associations as predictors.
# Correlation between species mean temperature index and species temperature range index requires separate models
# Model weights are the number of observed colonisations, persistences and extirpations for those species.

## Summary statistics
# Colonisations
mean(spe.df$pred_col) # mean 0.67
fivenum(spe.df$pred_col) # median 0.78, IQR 0.35-0.99
# Extirpations
mean(spe.df$pred_ext) # mean 0.27
fivenum(spe.df$pred_ext) # median 0.13, IQR 0-0.48
# Persistences
mean(spe.df$pred_per) # mean 0.82
fivenum(spe.df$pred_per) # median 0.93, IQR 0.71-1


### Modelling - Original analyses using data-derived thresholds for probability of occurrence being counted as presence in SDM predictions. Outputs correspond to Table S2.

## Colonisations
# Species mean temperature index
m.col.range<-glm(pred_col~grass.spec+forest.spec+SCI_TempRange, family="quasibinomial", weights=n.obs.col,data=spe.df)
mod.list[["col.stri"]]<-data.frame(mod="col.stri",pred=rownames(summary(m.col.range)$coefficients), summary(m.col.range)$coefficients, r2=with(summary(m.col.range), 1 - deviance/null.deviance))

# Species mean temperature index
m.col.mean<-glm(pred_col~grass.spec+forest.spec+SCI_TempMean, family="quasibinomial", weights=n.obs.col,data=spe.df)
mod.list[["col.smti"]]<-data.frame(mod="col.smti",pred=rownames(summary(m.col.mean)$coefficients),summary(m.col.mean)$coefficients, r2=with(summary(m.col.mean), 1 - deviance/null.deviance))

## Persistences
# Species temperature range index
m.per.range<-glm(pred_per~grass.spec+forest.spec+SCI_TempRange, family="quasibinomial", weights=n.obs.per,data=spe.df)
mod.list[["per.stri"]]<-data.frame(mod="per.stri",pred=rownames(summary(m.per.range)$coefficients),summary(m.per.range)$coefficients, r2=with(summary(m.per.range), 1 - deviance/null.deviance))

# Species mean temperature index
m.per.mean<-glm(pred_per~grass.spec+forest.spec+SCI_TempMean, family="quasibinomial", weights=n.obs.per,data=spe.df)
mod.list[["per.smti"]]<-data.frame(mod="per.smti",pred=rownames(summary(m.per.mean)$coefficients),summary(m.per.mean)$coefficients, r2=with(summary(m.per.mean), 1 - deviance/null.deviance))

## Extirpations
# Species temperature range index
m.ext.range<-glm(pred_ext~grass.spec+forest.spec+SCI_TempRange, family="quasibinomial", weights=n.obs.ext,data=spe.df)
mod.list[["ext.stri"]]<-data.frame(mod="ext.stri",pred=rownames(summary(m.ext.range)$coefficients),summary(m.ext.range)$coefficients, r2=with(summary(m.ext.range), 1 - deviance/null.deviance))

# Species mean temperature index
m.ext.mean<-glm(pred_ext~grass.spec+forest.spec+SCI_TempMean, family="quasibinomial", weights=n.obs.ext,data=spe.df)
mod.list[["ext.smti"]]<-data.frame(mod="ext.smti",pred=rownames(summary(m.ext.mean)$coefficients),summary(m.ext.mean)$coefficients, r2=with(summary(m.ext.mean), 1 - deviance/null.deviance))



## Modelling - Analyses using data-derived thresholds for probability of occurrence being counted as presence in SDM predictions, including Shoener's D estimation of niche overlap. Outputs correspond to Table S3.

## Colonisations
# Species mean temperature index
m.col.range<-glm(pred_col~grass.spec+forest.spec+SCI_TempRange+Schoeners.D, family="quasibinomial", weights=n.obs.col,data=spe.df)
mod.list[["col.stri.niche"]]<-data.frame(mod="col.stri",pred=rownames(summary(m.col.range)$coefficients), summary(m.col.range)$coefficients, r2=with(summary(m.col.range), 1 - deviance/null.deviance))

# Species mean temperature index
m.col.mean<-glm(pred_col~grass.spec+forest.spec+SCI_TempMean+Schoeners.D, family="quasibinomial", weights=n.obs.col,data=spe.df)
mod.list[["col.smti.niche"]]<-data.frame(mod="col.smti",pred=rownames(summary(m.col.mean)$coefficients),summary(m.col.mean)$coefficients, r2=with(summary(m.col.mean), 1 - deviance/null.deviance))

## Persistences
# Species temperature range index
m.per.range<-glm(pred_per~grass.spec+forest.spec+SCI_TempRange+Schoeners.D, family="quasibinomial", weights=n.obs.per,data=spe.df)
mod.list[["per.stri.niche"]]<-data.frame(mod="per.stri",pred=rownames(summary(m.per.range)$coefficients),summary(m.per.range)$coefficients, r2=with(summary(m.per.range), 1 - deviance/null.deviance))

# Species mean temperature index
m.per.mean<-glm(pred_per~grass.spec+forest.spec+SCI_TempMean+Schoeners.D, family="quasibinomial", weights=n.obs.per,data=spe.df)
mod.list[["per.smti.niche"]]<-data.frame(mod="per.smti",pred=rownames(summary(m.per.mean)$coefficients),summary(m.per.mean)$coefficients, r2=with(summary(m.per.mean), 1 - deviance/null.deviance))

## Extirpations
# Species temperature range index
m.ext.range<-glm(pred_ext~grass.spec+forest.spec+SCI_TempRange+Schoeners.D, family="quasibinomial", weights=n.obs.ext,data=spe.df)
mod.list[["ext.stri.niche"]]<-data.frame(mod="ext.stri",pred=rownames(summary(m.ext.range)$coefficients),summary(m.ext.range)$coefficients, r2=with(summary(m.ext.range), 1 - deviance/null.deviance))

# Species mean temperature index
m.ext.mean<-glm(pred_ext~grass.spec+forest.spec+SCI_TempMean+Schoeners.D, family="quasibinomial", weights=n.obs.ext,data=spe.df)
mod.list[["ext.smti.niche"]]<-data.frame(mod="ext.smti",pred=rownames(summary(m.ext.mean)$coefficients),summary(m.ext.mean)$coefficients, r2=with(summary(m.ext.mean), 1 - deviance/null.deviance))



## Modelling - Sensitivity analyses looping through lower (data-derived minus 0.1 on 0-1 scale) and upper (data-derived plus 0.1) thresholds for probability of occurrence being counted as presence in SDM predictions. Data-derived thresholds not included in loop because those model objects are needed for plotting figures (below). Outputs correspond to Tables S5 and S6.

for(thresh in c(".lwr",".upr")){

## Colonisations
# Species mean temperature index
	
m.col.range.sens<-glm(get(paste0("pred_col",thresh))~grass.spec+forest.spec+SCI_TempRange, family="quasibinomial", weights=n.obs.col,data=spe.df)
mod.list[[paste0("col.stri",thresh)]]<-data.frame(mod=paste0("col.stri",thresh),pred=rownames(summary(m.col.range.sens)$coefficients), summary(m.col.range.sens)$coefficients, r2=with(summary(m.col.range.sens), 1 - deviance/null.deviance))

# Species mean temperature index
m.col.mean.sens<-glm(get(paste0("pred_col",thresh))~grass.spec+forest.spec+SCI_TempMean, family="quasibinomial", weights=n.obs.col,data=spe.df)
mod.list[[paste0("col.smti",thresh)]]<-data.frame(mod=paste0("col.smti",thresh),pred=rownames(summary(m.col.mean.sens)$coefficients),summary(m.col.mean.sens)$coefficients, r2=with(summary(m.col.mean.sens), 1 - deviance/null.deviance))

## Persistences
# Species temperature range index
m.per.range.sens<-glm(get(paste0("pred_per",thresh))~grass.spec+forest.spec+SCI_TempRange, family="quasibinomial", weights=n.obs.per,data=spe.df)
mod.list[[paste0("per.stri",thresh)]]<-data.frame(mod=paste0("per.stri",thresh),pred=rownames(summary(m.per.range.sens)$coefficients),summary(m.per.range.sens)$coefficients, r2=with(summary(m.per.range.sens), 1 - deviance/null.deviance))

# Species mean temperature index
m.per.mean.sens<-glm(get(paste0("pred_per",thresh))~grass.spec+forest.spec+SCI_TempMean, family="quasibinomial", weights=n.obs.per,data=spe.df)
mod.list[[paste0("per.smti",thresh)]]<-data.frame(mod=paste0("per.smti",thresh),pred=rownames(summary(m.per.mean.sens)$coefficients),summary(m.per.mean.sens)$coefficients, r2=with(summary(m.per.mean.sens), 1 - deviance/null.deviance))

## Extirpations
# Species temperature range index
m.ext.range.sens<-glm(get(paste0("pred_ext",thresh))~grass.spec+forest.spec+SCI_TempRange, family="quasibinomial", weights=n.obs.ext,data=spe.df)
mod.list[[paste0("ext.stri",thresh)]]<-data.frame(mod=paste0("ext.stri",thresh),pred=rownames(summary(m.ext.range.sens)$coefficients),summary(m.ext.range.sens)$coefficients, r2=with(summary(m.ext.range.sens), 1 - deviance/null.deviance))

# Species mean temperature index
m.ext.mean.sens<-glm(get(paste0("pred_ext",thresh))~grass.spec+forest.spec+SCI_TempMean, family="quasibinomial", weights=n.obs.ext,data=spe.df)
mod.list[[paste0("ext.smti",thresh)]]<-data.frame(mod=paste0("ext.smti",thresh),pred=rownames(summary(m.ext.mean.sens)$coefficients),summary(m.ext.mean.sens)$coefficients, r2=with(summary(m.ext.mean.sens), 1 - deviance/null.deviance))

}


#############################################
########### Grid-square models ##############
#############################################



### Grid-square level prediction of colonisations, extirpations and persistences using independent validation
# Models have predictive accuracy as a response, and species habitat and climate associations as predictors.
# Model weights are the number of observed colonisations, persistences and extirpations for those grid squares.

# Sequential regression - adding residual latitude effect removing the effect of climate
sq.df$lat.res<-resid(lm(lat~temp_change, data=sq.df))


## Models for data-derived predicted presence. Outputs correspond to Table S4.
## Predicting persistence

m.per.grid<-glm(frac_per~temp_change+micro_sd+grass_aba+lat.res+pcnm1+pcnm2, family="quasibinomial", weights=per, data=sq.df)
mod.list[["grid.per"]]<-data.frame(mod="grid.per",pred=rownames(summary(m.per.grid)$coefficients),summary(m.per.grid)$coefficients, r2=with(summary(m.per.grid), 1 - deviance/null.deviance))

## Predicting colonisation
m.col.grid<-glm(frac_col~temp_change+micro_sd+grass_aba+lat.res+pcnm1+pcnm2, family="quasibinomial", weights=col, data=sq.df)
mod.list[["grid.col"]]<-data.frame(mod="grid.col",pred=rownames(summary(m.col.grid)$coefficients),summary(m.col.grid)$coefficients, r2=with(summary(m.col.grid), 1 - deviance/null.deviance))

## Predicting extirpations
m.ext.grid<-glm(frac_ext~temp_change+micro_sd+grass_aba+lat.res+pcnm1+pcnm2, family="quasibinomial", weights=ext, data=sq.df)
mod.list[["grid.ext"]]<-data.frame(mod="grid.ext",pred=rownames(summary(m.ext.grid)$coefficients),summary(m.ext.grid)$coefficients, r2=with(summary(m.ext.grid), 1 - deviance/null.deviance))


## Model loop for lower and upper thresholds for predicted presence (sensitivity analysis). Outputs correspond to Tables S6 and S7.

for(thresh in c("_lwr","_upr")){

## Predicting persistence
m.per.grid.sens<-glm(get(paste0("frac_per",thresh))~temp_change+micro_sd+grass_aba+lat.res+pcnm1+pcnm2, family="quasibinomial", weights=per, data=sq.df)
mod.list[[paste0("grid.per",thresh)]]<-data.frame(mod=paste0("grid.per",thresh),pred=rownames(summary(m.per.grid.sens)$coefficients),summary(m.per.grid.sens)$coefficients, r2=with(summary(m.per.grid.sens), 1 - deviance/null.deviance))

## Predicting colonisation
m.col.grid.sens<-glm(get(paste0("frac_col",thresh))~temp_change+micro_sd+grass_aba+lat.res+pcnm1+pcnm2, family="quasibinomial", weights=col, data=sq.df)
mod.list[[paste0("grid.col",thresh)]]<-data.frame(mod=paste0("grid.col",thresh),pred=rownames(summary(m.col.grid.sens)$coefficients),summary(m.col.grid.sens)$coefficients, r2=with(summary(m.col.grid.sens), 1 - deviance/null.deviance))

## Predicting extirpations
m.ext.grid.sens<-glm(get(paste0("frac_ext",thresh))~temp_change+micro_sd+grass_aba+lat.res+pcnm1+pcnm2, family="quasibinomial", weights=ext, data=sq.df)
mod.list[[paste0("grid.ext",thresh)]]<-data.frame(mod=paste0("grid.ext",thresh),pred=rownames(summary(m.ext.grid.sens)$coefficients),summary(m.ext.grid.sens)$coefficients, r2=with(summary(m.ext.grid.sens), 1 - deviance/null.deviance))
}

mod.list # have a look at the model outputs


#############################################
################# Figures ###################
#############################################

library(visreg) # bring in library for confidence bands, predicted values etc.

# Write function for colour transparency for use in the figures
transp<-function(col,alpha){
  trans.col<-adjustcolor(col,alpha.f = alpha)
  return(trans.col)
}

### Figure 2. TSS cross-validation vs. TSS independent validation

vr.tss<-visreg(m.tss.corr, plot=FALSE) # create visreg object
pdf("Figure2.pdf",width=8, height=8,useDingbats=F) # for pdf creation, if desired.
par(mar=c(6,6,2,4))
plot(visregRes ~ TSS_cross, data=vr.tss$res, bty="n", pch=16, col=transp("#D35FB7",0.4), ylab="TSS temporal", xlab="TSS cross-validation",cex.axis=1.4,  cex.lab=1.5)
lines(visregFit ~ TSS_cross, data=vr.tss$fit, col="#D35FB7", lwd=3)
polygon(c(vr.tss$fit$TSS_cross,rev(vr.tss$fit$TSS_cross)), c(vr.tss$fit$visregLwr,rev(vr.tss$fit$visregUpr)),col=transp("#D35FB7",0.5), lty=0)
dev.off()


# Create visreg objects for remaining (logistic) models, converting predicted values back to 0-1.
plot.mods<-c("m.per.mean", "m.per.range", "m.col.mean", "m.col.range","m.ext.mean", "m.ext.range", "m.per.grid")
for(mod in plot.mods){  
  for(predvar in names(get(mod)$coefficients)[2:length(names(get(mod)$coefficients))]){
temp.out <- visreg(get(mod), plot=FALSE)[[which(names(get(mod)$coefficients)==predvar)-1]]
temp.out[["fit"]][,c("visregFit","visregUpr","visregLwr")] <- sapply(temp.out[["fit"]][,c("visregFit","visregUpr","visregLwr")],plogis)
temp.out[["res"]][,"visregRes"] <- plogis(temp.out[["res"]][,"visregRes"])
assign(paste(substr(mod,3,nchar(mod)),predvar, sep="."),temp.out)
}}


### Figure 3. Effects of species associations on predictive accuracy

pdf("Figure3.pdf",width=15, height=15,useDingbats=F) # create pdf if wanted

# Set up parameters
par(mfrow=c(3,3), mar=c(6,6,2,2))

# Colours
grass.col<-"#E1BE6A"
forest.col<-"#40B0A6"
range.col<-"#E66100"
mean.col<-"#5D3A9B"

# Sizes of things
cax<-2.2
clab<-2.4
cleg<-1.8
xspec<-0.5
xsmti<-4
xstri<-20
ypos<-1.05
abcex<-2.6

## Colonisation
# Habitat specialisation
plot(visregRes ~ grass.spec, data=col.range.grass.spec$res, bty="n", pch=19, col=transp(grass.col,0.4), ylab="Colonisation prediction accuracy", xlab="Habitat specialisation",cex.axis=cax, cex.lab=clab)
lines(visregFit ~ grass.spec, data=col.range.grass.spec$fit, col=grass.col, lwd=3)
polygon(c(col.range.grass.spec$fit$grass.spec,rev(col.range.grass.spec$fit$grass.spec)), c(col.range.grass.spec$fit$visregLwr,rev(col.range.grass.spec$fit$visregUpr)),col=transp(grass.col,0.5), lty=0)

points(visregRes ~ forest.spec, data=col.range.forest.spec$res, bty="n", pch=19, col=transp(forest.col,0.4))
lines(visregFit ~ forest.spec, data=col.range.forest.spec$fit, col=forest.col, lwd=3)
polygon(c(col.range.forest.spec$fit$forest.spec,rev(col.range.forest.spec$fit$forest.spec)), c(col.range.forest.spec$fit$visregLwr,rev(col.mean.forest.spec$fit$visregUpr)),col=transp(forest.col,0.5), lty=0)

legend(4,0.35, c("Grassland", "Forest"), pch=19, col=c(transp(grass.col,0.7),transp(forest.col,0.5)), cex=cleg)
text(xspec,ypos,"(a)", font=2, cex=abcex, xpd=TRUE)

# Species temperature range index
plot(visregRes ~ SCI_TempRange, data=col.range.SCI_TempRange$res, bty="n", pch=19, col=transp(range.col,0.30),ylab="Colonisation prediction accuracy", xlab="Species temperature range index (°C)",cex.axis=cax,  cex.lab=clab)
lines(visregFit ~ SCI_TempRange, data=col.range.SCI_TempRange$fit, col=range.col, lwd=3)
polygon(c(col.range.SCI_TempRange$fit$SCI_TempRange,rev(col.range.SCI_TempRange$fit$SCI_TempRange)), c(col.range.SCI_TempRange$fit$visregLwr,rev(col.range.SCI_TempRange$fit$visregUpr)),col=transp(range.col,0.5), lty=0)
text(xstri,ypos,"(b)", font=2, cex=abcex, xpd=TRUE)

# Species mean temperature index
plot(visregRes ~ SCI_TempMean, data=col.mean.SCI_TempMean$res, bty="n", pch=19, col=transp(mean.col,0.30),ylab="Colonisation prediction accuracy", xlab="Species mean temperature index (°C)",cex.axis=cax,  cex.lab=clab)
lines(visregFit ~ SCI_TempMean, data=col.mean.SCI_TempMean$fit, col=mean.col, lwd=3)
polygon(c(col.mean.SCI_TempMean$fit$SCI_TempMean,rev(col.mean.SCI_TempMean$fit$SCI_TempMean)), c(col.mean.SCI_TempMean$fit$visregLwr,rev(col.mean.SCI_TempMean$fit$visregUpr)),col=transp(mean.col,0.5), lty=0)
text(xsmti,ypos,"(c)", font=2, cex=abcex, xpd=TRUE)


## Persistence
# Habitat specialisation
plot(visregRes ~ grass.spec, data=per.range.grass.spec$res, bty="n", pch=19, col=transp(grass.col,0.4), ylab="Persistence prediction accuracy", xlab="Habitat specialisation",cex.axis=cax,  cex.lab=clab)
lines(visregFit ~ grass.spec, data=per.range.grass.spec$fit, col=grass.col, lwd=3)
polygon(c(per.range.grass.spec$fit$grass.spec,rev(per.range.grass.spec$fit$grass.spec)), c(per.range.grass.spec$fit$visregLwr,rev(per.range.grass.spec$fit$visregUpr)),col=transp(grass.col,0.5), lty=0)

points(visregRes ~ forest.spec, data=per.range.forest.spec$res, bty="n", pch=19, col=transp(forest.col,0.4))
lines(visregFit ~ forest.spec, data=per.range.forest.spec$fit, col=forest.col, lwd=3)
polygon(c(per.range.forest.spec$fit$forest.spec,rev(per.range.forest.spec$fit$forest.spec)), c(per.range.forest.spec$fit$visregLwr,rev(per.mean.forest.spec$fit$visregUpr)),col=transp(forest.col,0.5), lty=0)

legend(4,0.35, c("Grassland", "Forest"), pch=19, col=c(transp(grass.col,0.7),transp(forest.col,0.5)), cex=cleg)
text(xspec,ypos,"(d)", font=2, cex=abcex, xpd=TRUE)

# Species temperature range index
plot(visregRes ~ SCI_TempRange, data=per.range.SCI_TempRange$res, bty="n", pch=19, col=transp(range.col,0.30),ylab="Persistence prediction accuracy", xlab="Species temperature range index (°C)",cex.axis=cax,  cex.lab=clab)
lines(visregFit ~ SCI_TempRange, data=per.range.SCI_TempRange$fit, col=range.col, lwd=3)
polygon(c(per.range.SCI_TempRange$fit$SCI_TempRange,rev(per.range.SCI_TempRange$fit$SCI_TempRange)), c(per.range.SCI_TempRange$fit$visregLwr,rev(per.range.SCI_TempRange$fit$visregUpr)),col=transp(range.col,0.5), lty=0)
text(xstri,ypos,"(e)", font=2, cex=abcex, xpd=TRUE)

# Species mean temperature index
plot(visregRes ~ SCI_TempMean, data=per.mean.SCI_TempMean$res, bty="n", pch=19, col=transp(mean.col,0.30),ylab="Persistence prediction accuracy", xlab="Species mean temperature index (°C)",cex.axis=cax,  cex.lab=clab)
lines(visregFit ~ SCI_TempMean, data=per.mean.SCI_TempMean$fit, col=mean.col, lwd=3)
polygon(c(per.mean.SCI_TempMean$fit$SCI_TempMean,rev(per.mean.SCI_TempMean$fit$SCI_TempMean)), c(per.mean.SCI_TempMean$fit$visregLwr,rev(per.mean.SCI_TempMean$fit$visregUpr)),col=transp(mean.col,0.5), lty=0)
text(xsmti,ypos,"(f)", font=2, cex=abcex, xpd=TRUE)


## Extirpation
# Habitat specialisation
plot(visregRes ~ grass.spec, data=ext.range.grass.spec$res, bty="n", pch=19, col=transp(grass.col,0.4), ylab="Extirpation prediction accuracy", xlab="Habitat specialisation",cex.axis=cax,  cex.lab=clab)
lines(visregFit ~ grass.spec, data=ext.range.grass.spec$fit, col=grass.col, lwd=3)
polygon(c(ext.range.grass.spec$fit$grass.spec,rev(ext.range.grass.spec$fit$grass.spec)), c(ext.range.grass.spec$fit$visregLwr,rev(ext.range.grass.spec$fit$visregUpr)),col=transp(grass.col,0.5), lty=0)

points(visregRes ~ forest.spec, data=ext.range.forest.spec$res, bty="n", pch=19, col=transp(forest.col,0.4))
lines(visregFit ~ forest.spec, data=ext.range.forest.spec$fit, col=forest.col, lwd=3)
polygon(c(ext.range.forest.spec$fit$forest.spec,rev(ext.range.forest.spec$fit$forest.spec)), c(ext.range.forest.spec$fit$visregLwr,rev(ext.mean.forest.spec$fit$visregUpr)),col=transp(forest.col,0.5), lty=0)

legend(4,0.9, c("Grassland", "Forest"), pch=19, col=c(transp(grass.col,0.7),transp(forest.col,0.5)), cex=cleg)
text(xspec,ypos,"(g)", font=2, cex=abcex, xpd=TRUE)

# Species temperature range index
plot(visregRes ~ SCI_TempRange, data=ext.range.SCI_TempRange$res, bty="n", pch=19, col=transp(range.col,0.30),ylab="Extirpation prediction accuracy", xlab="Species temperature range index (°C)",cex.axis=cax,  cex.lab=clab)
lines(visregFit ~ SCI_TempRange, data=ext.range.SCI_TempRange$fit, col=range.col, lwd=3)
polygon(c(ext.range.SCI_TempRange$fit$SCI_TempRange,rev(ext.range.SCI_TempRange$fit$SCI_TempRange)), c(ext.range.SCI_TempRange$fit$visregLwr,rev(ext.range.SCI_TempRange$fit$visregUpr)),col=transp(range.col,0.5), lty=0)
text(xstri,ypos,"(h)", font=2, cex=abcex, xpd=TRUE)

# Species mean temperature index
plot(visregRes ~ SCI_TempMean, data=ext.mean.SCI_TempMean$res, bty="n", pch=19, col=transp(mean.col,0.30),ylab="Extirpation prediction accuracy", xlab="Species mean temperature index (°C)",cex.axis=cax,  cex.lab=clab)
lines(visregFit ~ SCI_TempMean, data=ext.mean.SCI_TempMean$fit, col=mean.col, lwd=3)
polygon(c(ext.mean.SCI_TempMean$fit$SCI_TempMean,rev(ext.mean.SCI_TempMean$fit$SCI_TempMean)), c(ext.mean.SCI_TempMean$fit$visregLwr,rev(ext.mean.SCI_TempMean$fit$visregUpr)),col=transp(mean.col,0.5), lty=0)
text(xsmti,ypos,"(i)", font=2, cex=abcex, xpd=TRUE)
dev.off()



### Figure 4. Overall grid-square predictive accuracy

# Set up parameters
transp.val <- 0.25
ext.col <- "#FFB000"
per.col <- "#DC267F"
col.col <-  "#648FFF"

pdf("Figure4.pdf",width=7, height=5,useDingbats=F)
boxplot(sq.df$frac_col, sq.df$frac_ext, sq.df$frac_per, outline=FALSE, col="seashell2", border=FALSE,boxwex=0.6,staplewex=0, lty=1, lwd=0.1, boxlwd=0.1,frame.plot=FALSE,frame=FALSE, cex.axis=1.4, xaxt="n",yaxt="n", cex.lab=1.5, horizontal=TRUE,notch=TRUE)

axis(1, tick=TRUE, lwd=2, cex.axis=1.4)
title(xlab="Fraction correctly predicted", cex.lab=1.5)

points(cbind(sq.df$frac_per,jitter(rep(3,length(sq.df$frac_per)),1,0.25)),col=transp(per.col,transp.val), pch=16,lwd=2)
points(cbind(sq.df$frac_col,jitter(rep(1,length(sq.df$frac_col)),1,0.25)),col=transp(col.col,transp.val), pch=16,lwd=2)
points(cbind(sq.df$frac_ext,jitter(rep(2,length(sq.df$frac_ext)),1,0.25)),col=transp(ext.col,transp.val), pch=16,lwd=2)

boxplot(sq.df$frac_col, sq.df$frac_ext, sq.df$frac_per, outline=FALSE, col=transp("seashell2",0.1), border="black",ylab="",boxwex=0.6,staplewex=0, notch=TRUE, lty=1, lwd=3., boxlwd=3,frame.plot=FALSE, frame=FALSE, xaxt="n", yaxt="n", horizontal=TRUE, add=TRUE)

text(mean(fivenum(sq.df$frac_col)[c(2,4)]),1.4,"Colonisations", font=2, cex=1)
text(mean(fivenum(sq.df$frac_ext)[c(2,4)]),2.4,"Extirpations", font=2, cex=1)
text(mean(fivenum(sq.df$frac_per)[c(2,4)]),3.4,"Persistences", font=2, cex=1)
dev.off()


### Figure 5. Prediction of persistence according to experienced climate and land use change

# Set up parameters
ypos<-1.05
xaba<-0.01
xmat<-1.15
abcex<-1.75

pdf("Figure5.pdf",width=11, height=5,useDingbats=F)
par(mfrow=c(1,2), mar=c(6,6,2,2))
## Grassland abandonment
plot(visregRes ~ grass_aba, data=per.grid.grass_aba$res, bty="n", pch=16, col=transp("#1AFF1A",0.4),ylab="Predicted persistences", xlab="Grassland abandonment",cex.axis=1.4,  cex.lab=1.5)
lines(visregFit ~ grass_aba, data=per.grid.grass_aba$fit, col="#1AFF1A", lwd=3)
polygon(c(per.grid.grass_aba$fit$grass_aba,rev(per.grid.grass_aba$fit$grass_aba)), c(per.grid.grass_aba$fit$visregLwr,rev(per.grid.grass_aba$fit$visregUpr)),col=transp("#1AFF1A",0.7), lty=0)
text(xaba,ypos,"(a)", font=2, cex=abcex, xpd=TRUE)

## Temperature change
plot(visregRes ~ temp_change, data=per.grid.temp_change$res, bty="n", pch=16, col=transp("#4B0092",0.30),ylab="Predicted persistences", xlab="Change in mean annual temperature (°C)",cex.axis=1.4,  cex.lab=1.5)
lines(visregFit ~ temp_change, data=per.grid.temp_change$fit, col="#4B0092", lwd=3)
polygon(c(per.grid.temp_change$fit$temp_change,rev(per.grid.temp_change$fit$temp_change)), c(per.grid.temp_change$fit$visregLwr,rev(per.grid.temp_change$fit$visregUpr)),col=transp("#4B0092",0.5), lty=0)
text(xmat,ypos,"(b)", font=2, cex=abcex, xpd=TRUE)

dev.off()


### FIGURE S1
pdf("FigureS1.pdf",width=11, height=5,useDingbats=F)
par(mar=c(7,4,0,2))

boxplot(spe.df$Schoeners.D, outline=FALSE, ylim=c(0,1),col=NULL, border=FALSE,boxwex=0.6,staplewex=0, lty=1, lwd=0.1, boxlwd=0.1,frame.plot=FALSE, frame=FALSE, cex.axis=1.4, xaxt="n",yaxt="n", cex.lab=1.5,  horizontal=TRUE)

axis(1, tick=TRUE, lwd=2, cex.axis=1.4,  font=2)
title(xlab="Niche overlap (Schoener's D)", cex.lab=1.5, font.lab=2)

points(cbind(spe.df$Schoeners.D[spe.df$Schoeners.D.p>0.05],jitter(rep(0.75,length(spe.df$Schoeners.D[spe.df$Schoeners.D.p>0.05])),1,0.15)),col="forestgreen", pch=21,lwd=2)
points(cbind(spe.df$Schoeners.D[spe.df$Schoeners.D.p<=0.05],jitter(rep(0.75,length(spe.df$Schoeners.D[spe.df$Schoeners.D.p<0.05])),1,0.15)),col="forestgreen", bg="forestgreen", pch=21,lwd=2)

dev.off()


### FIGURE S2. Summary of MESS analysis.
pdf("FigureS2.pdf",width=12, height=5,useDingbats=F)
par(mfrow=c(1,2), mar=c(5,5,2,4))

hist(sq.df$MESS, main="", xlab="MESS value",cex.lab=1.4)
text(-135,210,"(a)", font=2, cex=1.5, xpd=TRUE)

hist(sq.df$t.max.mod, xlim=c(14.5,19.5), ylim=c(0,230),col="#FFC20A", main="", xlab="Mean of monthly maximum temperatures (°C)",cex.lab=1.4)
hist(sq.df$t.max.hist, add=TRUE, col="#0C7BDC")
legend(18.15, 230, c("1961-1970", "2001-2010"), col=c("#0C7BDC","#FFC20A" ), bty="n", pch=15, pt.cex=1.25, cex=1)
text(14.75,210,"(b)", font=2, cex=1.5, xpd=TRUE)

dev.off()

table(sq.df$MESS.n.var) # 1015 squares with one out of range value, 23 with two (t max and open - some historical containing no open).



cat("\n\n")
sessionInfo()