Iturbide et al 2016 (submitted)

data(Oak_phylo2) data(biostack)

# Extract raster values at grid coordinates
ac <- xyFromCell(biostack$bio13, 1:ncell(biostack$bio13))
ex <- extract(biostack$bio13, 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")

bc <- rep(list(c(-10, 35, 30, 70)), times = length(Oak_phylo2))
del <- delimit(bounding.coords = bc, grid = sp_grid_world , names = names(Oak_phylo2))
bg <- del$bbs.grid

#TS method

#step1
unsuitable.bg <-OCSVMprofiling(xy = Oak_phylo2,
                                     varstack = biostack,
                                     bbs.grid = bg)

#step2
ext <- bgRadio(xy = Oak_phylo2, bounding.coords = bc,
                     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 = Oak_phylo2, bg.grids = ext,
                              exclusion.buffer = 0.249,
                              prevalence = 0.5, kmeans = FALSE) 
}

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

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

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

for(i in 1:10){
  modirs <- allModeling(data = presaus.TS[[i]], varstack = biostack,
                        
                        k = 10, algorithm = "maxent", destdir = dirs[i],
                        
                        projection = CRS("+proj=longlat +init=epsg:4326"))
}

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

for(i in 1:10){
  auc_mars_TS[[i]] <-loadTestValues(data = presaus.TS[[i]], test = "auc",sourcedir = dirs[i],
                                    algorithm = "mars")
  auc_glm_TS[[i]] <-loadTestValues(data = presaus.TS[[i]], test = "auc",sourcedir = dirs[i],
                                   algorithm = "glm")
  auc_maxent_TS[[i]] <-loadTestValues(data = presaus.TS[[i]], test = "auc", sourcedir = dirs[i],
                                      algorithm = "maxent")
  
}

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

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