Iturbide et al 2016 (submitted)

In this tutorial we are using the same presence data used in other sections of this wiki, this is distributions of Quercus robur phylogenies H11 and H5 in Europe. Regarding climatic variables, a subset of the data used in the manuscript is used (E-OBS for the baseline data and the regional climate projections from SMHI_ECHAM5). This worked example is fully reproducible, as the data is available with the mopa package.

data(Oak_phylo2)
data(biostackENSEMBLES)

# Extract raster values at grid coordinates
ac <- xyFromCell(biostackENSEMBLES$baseline$bio2, 1:ncell(biostackENSEMBLES$baseline$bio2))
ex <- extract(biostackENSEMBLES$baseline$bio2, ac)

# Convert to a Spatial object and define projection
sp_grid_world <- SpatialPoints(ac[-which(is.na(ex)), ])
projection(sp_grid_world) <- CRS("+proj=longlat +init=epsg:4326")

del <- delimit(bounding.coords = c(-10, 35, 30, 70), grid = sp_grid_world)
bg <- del$bbs.grid

#TS method

#step1
unsuitable.bg <-OCSVMprofiling(xy = presences,
                               varstack = biostackENSEMBLES$baseline,
                               bbs.grid = bg)

#step2
ext <- bgRadio(xy = presences, bounding.coords = c(-10, 35, 30, 70),
               bg.absence = unsuitable.bg$absence,
               start = 0.249, by = 0.249,
               unit = "decimal degrees")

pa_TS <- list()
for (i in 1:10){
  pa_TS[[i]] <-PseudoAbsences(xy = presences, bg.grids = ext,
                              exclusion.buffer = 0.249,
                              prevalence = 0.5, kmeans = FALSE) 
}

presaus.TS <- list()
for(i in 1:10){
presaus.TS[[i]] <- bindPresAbs(presences = presences, absences = pa_TS[[i]])
}

dirs <- paste("mydirectory/modTS/round",1:10, sep="")
sapply(1:10, function(x){dir.create(dirs[x])})

for(i in 1:10){
  modmarsTS <- allModeling(data = presaus.TS[[i]], varstack = biostackENSEMBLES$baseline,
                          
                          k = 10, algorithm = "mars", destdir = dirs[i],
                          
                          projection = CRS("+proj=longlat +init=epsg:4326"))
  
  modglmTS <- allModeling(data = presaus.TS[[i]], varstack = biostackENSEMBLES$baseline,
                          
                          k = 10, algorithm = "glm", destdir = dirs[i],
                          
                          projection = CRS("+proj=longlat +init=epsg:4326"))
}

auc_mars_TS <- list()
auc_glm_TS <- list()

for(i in 1:10){
  auc_mars_TS[[i]] <-loadTestValues(models = modmarsTS, test = "auc",sourcedir = dirs[i])
  auc_glm_TS[[i]] <-loadTestValues(models = modglmTS, test = "auc", sourcedir = dirs[i])
}

ind.mars.TS <- list()
ind.glm.TS <- list()

for(i in 1:10){
  ind.mars.TS[[i]] <- indextent(testmat = auc_mars_TS[[i]], diagrams = TRUE)
  ind.glm.TS[[i]] <- indextent(testmat = auc_glm_TS[[i]], diagrams = TRUE)
}

mars.models <- list()
glm.models <- list()

for(i in 1:10){
mars.models[[i]] <- loadDefinitiveModel(ind.mars.TS[[i]], slot = "allmod")
glm.models[[i]] <- loadDefinitiveModel(ind.glm.TS[[i]], slot = "allmod")
}

#predictions
prob.mars <- list()
prob.glm <- list()
for(i in 1:10){
  prob.mars[[i]] <- predictAllmod(model = mars.models[[i]], environment = biostackENSEMBLES$future_71_100_SMHI_ECHAM5)
  prob.glm[[i]] <- predictAllmod(glm.models[[i]], biostackENSEMBLES$future_71_100_SMHI_ECHAM5)
}

prob.data <- list("mars" = prob.mars, "glm" = prob.glm)

datarange <- list()
for (i in 1:length(prob.data)){
  maps <- list()
  data <- prob.data[[i]]
  H11 <- list()
  H5 <- list()
    for(j in 1:10){
      H11[[j]] <- data[[j]]$H11
      H5[[j]] <- data[[j]]$H5
    }
    names(H11) <- paste("round", 1:10, sep="")
    names(H5) <- names(H11)
    H11 <- calc(stack(H11), max, na.rm=T) - calc(stack(H11), min, na.rm=T)
    H5 <- calc(stack(H5), max, na.rm=T) - calc(stack(H5), min, na.rm=T)
    H11@data@values[which(is.infinite(H11@data@values))] <-NA
    H5@data@values[which(is.infinite(H5@data@values))] <-NA
    st <- stack(H11, H5)
    names(st) <- c("H11", "H5")  
    datarange[[i]] <- st
}
names(datarange) <- names(prob.data)

spplot(datarange$mars, at = seq(0,1,.1), col.regions = rev(gray.colors(10, end = 1, start = 0, gamma = 1.2)), strip = T)

spplot(datarange$glm, at = seq(0,1,.1), col.regions = rev(gray.colors(10, end = 1, start = 0, gamma = 1.2)), strip = T)

par(mfrow = c(1,2))
gridboxes <- length(na.omit(datarange$mars$H11@data@values))

#phylogeny H11
mars <- cbind(1:gridboxes *100/gridboxes , sort(na.omit(datarange$mars$H11@data@values)))
glm <- cbind(1:gridboxes *100/gridboxes , sort(na.omit(datarange$glm$H11@data@values)))
plot(mars[,2], mars[,1],xlim= c(0,1), ylim = c(0,100), col = "red", typ = "l", ylab = "gridboxes (%)", xlab = "range", main = "H11")
lines(glm[,2], glm[,1], typ = "l")
segments(x0 = seq(0.2,0.8,0.2), y0 = 0, x1 = seq(0.2,0.8,0.2), y1 = 100, col = "grey", lty = 3)
segments(x0 = 0, y0 = seq(20,80,20), x1 = 1, y1 = seq(20,80,20), col = "grey", lty = 3)
mars.5 <- round(100 - mars[which(mars[,2]>0.5),][1,][1], digits = 1)
glm.5 <- round(100 - glm[which(glm[,2]>0.5),][1,][1], digits = 1)
mars.25 <- round(100 - mars[which(mars[,2]>0.25),][1,][1], digits = 1)
glm.25 <- round(100 - glm[which(glm[,2]>0.25),][1,][1], digits = 1)
segments(x0 = 0.5, y0 = 0, x1 = 0.5, y1 = 100, col = "blue")
segments(x0 = 0.25, y0 = 0, x1 = 0.25, y1 = 100, col = "green")
legend(x = 0.51, y = 45, legend = paste(c(mars.5, glm.5),"   ",c(mars.25, glm.25), sep=" "),text.col = "green", lty = 1, col = c("red", "black"), border = F, box.lwd = 0, cex = .8, bty = "n")
legend(x = 0.51, y = 45, legend = c(mars.5, glm.5), text.col = "blue", lty = 1, col = c("red","black"), border = F, box.lwd = 0, cex = .8, bty = "n")

#phylogeny H5  
mars <- cbind(1:gridboxes *100/gridboxes , sort(na.omit(datarange$mars$H5@data@values)))
glm <- cbind(1:gridboxes *100/gridboxes , sort(na.omit(datarange$glm$H5@data@values)))
plot(mars[,2], mars[,1], xlim= c(0,1), ylim = c(0,100), col = "red", typ = "l", ylab = "gridboxes (%)", xlab = "range", main = "H5")
lines(glm[,2], glm[,1], typ = "l")
segments(x0 = seq(0.2,0.8,0.2), y0 = 0, x1 = seq(0.2,0.8,0.2), y1 = 100, col = "grey", lty = 3)
segments(x0 = 0, y0 = seq(20,80,20), x1 = 1, y1 = seq(20,80,20), col = "grey", lty = 3)
mars.5 <- round(100 - mars[which(mars[,2]>0.5),][1,][1], digits = 1)
glm.5 <- round(100 - glm[which(glm[,2]>0.5),][1,][1], digits = 1)
mars.25 <- round(100 - mars[which(mars[,2]>0.25),][1,][1], digits = 1)
glm.25 <- round(100 - glm[which(glm[,2]>0.25),][1,][1], digits = 1)
segments(x0 = 0.5, y0 = 0, x1 = 0.5, y1 = 100, col = "blue")
segments(x0 = 0.25, y0 = 0, x1 = 0.25, y1 = 100, col = "green")
legend(x = 0.51, y = 45, legend = paste(c(mars.5, glm.5),"   ",c(mars.25, glm.25), sep=" "),text.col = "green", lty = 1, col = c("red", "black"), border = F, box.lwd = 0, cex = .8, bty = "n")
legend(x = 0.51, y = 45, legend = c(mars.5, glm.5), text.col = "blue", lty = 1, col = c("red","black"), border = F, box.lwd = 0, cex = .8, bty = "n")

legend(x = 0.5, y = 95, legend = c("mars", "glm"), lty = 1, col = c("red","black"), bty = "n",
       border = F, box.lwd = 0, cex = 1)

mars.H11 <- list()
glm.H11 <- list()
for(i in 1:10){
  mars.H11[[i]] <- prob.mars[[i]]$H11
  glm.H11[[i]] <- prob.glm[[i]]$H11
}
names(mars.H11) <- as.character(1:10)
names(glm.H11) <- names(mars.H11)
mars.H11 <- stack(mars.H11)
glm.H11 <- stack(glm.H11)

ov.mars.H11 <- nicheOver(stack = mars.H11, metric = "D")
library(lattice)
levelplot(ov.mars.H11)

ind <- which(ov.mars.H11 == min(ov.mars.H11, na.rm = T), arr.ind = T)
extremes <- stack(list(prob.mars[[ind[1]]]$H11, prob.mars[[ind[2]]]$H11))
spplot(extremes, strip = F)

ov.glm.H11 <- nicheOver(stack = glm.H11, metric = "D")
ind <- which(ov.glm.H11 == min(ov.glm.H11, na.rm = T), arr.ind = T)
extremes <- stack(list(prob.glm[[ind[1]]]$H11, prob.glm[[ind[2]]]$H11))
spplot(extremes, strip = F)