# Supporting R code

# Thin by 2km
A_ferreus_unfiltered<-read.csv("")
A_ferreus_filtered <- thin(loc.data=A_ferreus_unfiltered, lat.col = "latitude", long.col = "longitude", spec.col = "species", thin.par=2, reps=100, locs.thinned.list.return = FALSE, write.files = TRUE, max.files = 5, out.dir="F:/Tom_EnvLayers/RapoportsRuleSalamander/A.ferreus/thin", out.base = "A_ferreus_thinned_data", write.log.file = TRUE, log.file = "A_ferreus_thin_log.txt")

# Optimize ENM settings
bio_1 <- raster("")
bio_2 <- raster("")
bio_4 <- raster("")
bio_7 <- raster("")
bio_8 <- raster("")
bio_9 <- raster("")
bio_12 <- raster("")
bio_15 <- raster("")
bio_16 <- raster("")
env_7<-stack(bio_1, bio_2, bio_4, bio_7, bio_12, bio_15, bio_16)
env_6<-stack(bio_1, bio_4, bio_8, bio_9, bio_12, bio_15)

bg_ferreus<-randomPoints(env_7[[1]], n=10000)
bg_ferreus<-as.data.frame(bg_ferreus)
ferreus_blocks<-get.block(A_ferreus, bg_ferreus)
eval_ferreus_7<-ENMevaluate(occ=A_ferreus, env=env_7, bg.coords=bg_ferreus, method='block', RMvalues=c(0.5, 1, 1.5, 2), fc=c('LQ', 'LQH', 'LQHP', 'LQHPT'), rasterPreds=TRUE)
eval_ferreus_7@results
# By deltaAICc, LQHP_2. By AUC, LQ_1
eval_ferreus_6<-ENMevaluate(occ=A_ferreus, env=env_6, bg.coords=bg_ferreus, method='block', RMvalues=c(0.5, 1, 1.5, 2), fc=c('LQ', 'LQH', 'LQHP', 'LQHPT'), rasterPreds=TRUE)

# Null model via randomly matching mid-latitude and latitudinal extents.
all_latspan<-read.csv("mid_lat_temperate.csv", header=TRUE, row.names="species")

myfunction1<-function(x){
  null_range<-sample(x$lat_span, 169, replace=FALSE) # 169 species in total
  null_midpoint<-sample(x$lat_mid, 169, replace=FALSE)
  null_min<-null_midpoint-(0.5*null_range)
  null_max<-null_midpoint+(0.5*null_range)
  null_bind<-as.data.frame(cbind(null_range, null_midpoint, null_min, null_max))
  null_bind2<-subset(null_bind, null_min > 25.5) # lower latitudinal limit
  null_bind2<-as.data.frame(null_bind2)
  null_bind3<-subset(x=null_bind2, null_max < 57.5) # upper latitudinal limit
  null_bind3<-as.data.frame(null_bind3)
  Null_model<-cor.test(null_bind3$null_range, null_bind3$null_midpoint, method="spearman")
}

results1<-c()
for (i in 1:1000){output <- myfunction1(all_latspan)
results1<-c(results1, output$estimate)}
hist(results1, breaks=100, main=NULL, xlab="Correlation Coefficient", xlim=c(-0.4, 0.4))
abline(v=0.34, lty=2)
length(results1[abs(results1) < 0.34]) # compare null values of 'r' to empirical value

# Null model by removing species from analysis compared to when only northernmost species are removed
## asmple 25 at a time
myfunction2<-function(x){
  a <- sample_n(tbl=x, size=144, replace=TRUE)
  b <- a$lat_span
  c <- a$lat_mid
  d <- cor.test(x=b, y=c, method="spearman")
  d
}

for (i in 1:1000){myfunction2(all_latspan)}

f<-myfunction2(all_latspan)
f$estimate

results2<-c()
for (i in 1:1000){output2 <- myfunction2(x=all_latspan)
results2<-c(results2, output2$estimate)}
hist(results2, breaks=100, main=NULL, xlab="Correlation Coefficient", xlim=c(0, 0.7))
abline(v=0.15, lty=2)
length(results2[(results2) < 0.15])