library(BIEN)
library(dismo)
library(ENMeval)
library(maxnet)
library(raster)
library(rworldxtra)
library(sf)
library(tidyverse)

• Entender que modelan las SDM
• Diferenciando Maxent de otros algoritmos
  • Parametros de regularización
  • Pseudoaucencias vs Background
• Respuestas (Lineal, Cuadrática, Bisagra, Producto)
• Recomendaciones según número de presencias

• Distribución potencial, no real
• No considera disperción
• No es tan fácil dividir en presencia y ausencia

• Paquete BIEN para plantas
• Para animales y otros rgbif, pero hay que hacerse cuenta

N_pumilio <- BIEN_occurrence_species(species = "Nothofagus pumilio")

N_pumilio_SF <- N_pumilio %>% st_as_sf(coords = c(3, 2), crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
Hull <- N_pumilio_SF %>% st_union() %>% st_convex_hull()
Buffer <- Hull %>% st_buffer(dist = 1) %>% st_as_sf()
data("countriesHigh")
SA <- countriesHigh %>% st_as_sf() %>% st_make_valid() %>% st_crop(Buffer)
Bioclimatic <- getData(name = "worldclim", var = "bio", res = 0.5, 
    lon = -70, lat = -50)
## capas climaticas
Bioclimatic <- Bioclimatic %>% crop(Buffer) %>% trim()
## Cortamos el Tile
names(Bioclimatic) <- str_remove_all(names(Bioclimatic), "_43")

set.seed(2020)
Pres <- N_pumilio %>% dplyr::select(longitude, latitude) %>% 
    mutate(Pres = 1)
Aus <- dismo::randomPoints(mask = Bioclimatic[[1]], n = 5000) %>% 
    as.data.frame() %>% rename(longitude = x, latitude = y) %>% 
    mutate(Pres = 0)
Pres_Aus <- bind_rows(Pres, Aus)
Pres_Aus_Sf <- Pres_Aus %>% st_as_sf(coords = c(1, 2), crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")

Condiciones <- raster::extract(Bioclimatic, Pres_Aus_Sf) %>% 
    as.data.frame() %>% bind_cols(Pres_Aus)

• La forma mas sencilla maxnet

set.seed(2020)
Mod1 <- maxnet(p = Condiciones$Pres, data = Condiciones[, 1:19], 
    regmult = 1, maxnet.formula(p = Condiciones$Pres, data = Condiciones[, 
        1:19], classes = "lq"))

• Multiplicador de regularización (Número)
• Classes entre “lqpt”

plot(Mod1, c("bio1", "bio2", "bio3"))

plot(Mod1, type = "cloglog", c("bio1", "bio2", "bio3", "bio12"))

Mod2 <- maxnet(p = Condiciones$Pres, data = Condiciones[, 1:19], 
    regmult = 1, maxnet.formula(p = Condiciones$Pres, data = Condiciones[, 
        1:19], classes = "l"))
Mod3 <- maxnet(p = Condiciones$Pres, data = Condiciones[, 1:19], 
    regmult = 1, maxnet.formula(p = Condiciones$Pres, data = Condiciones[, 
        1:19], classes = "lh"))

plot(Mod2, c("bio1", "bio2", "bio3", "bio12"))

plot(Mod3, c("bio1", "bio2", "bio3", "bio12"))

• Si las presencias son < 10 “l”
• Si las presencias son < 15 “lq”
• Si las presencias son < 80 “lqh”
• Si las presencias son >= 80 “lqph”

Mod4 <- maxnet(p = Condiciones$Pres, data = Condiciones[, 1:19], 
    regmult = 1)

plot(Mod4, c("bio1", "bio2", "bio3", "bio12"))

• Maxent es como un modelo binomial
• Parte con una estimación para todas las variable
• Va penalizando cada variable según su importancia
• Quedamos con variables de valor 0
• Se corta con crossvalidation

Prediction <- predict(Bioclimatic, Mod4, type = "cloglog")

Eval <- dismo::evaluate(p = Pres[, 1:2], a = Aus[, 1:2], model = Mod4, 
    x = Bioclimatic, type = "cloglog")
EvalDF <- Eval@confusion %>% as.data.frame %>% mutate(Threshold = Eval@t) %>% 
    mutate(TP_TN = (tp/nrow(N_pumilio)) + (tn/5000))
head(EvalDF)

##    tp   fp fn tn  Threshold TP_TN
## 1 105 5000  0  0 -1.000e-04     1
## 2 105 5000  0  0 -9.999e-05     1
## 3 105 5000  0  0 -9.998e-05     1
## 4 105 5000  0  0 -9.997e-05     1
## 5 105 5000  0  0 -9.996e-05     1
## 6 105 5000  0  0 -9.995e-05     1

EvalThres <- EvalDF %>% dplyr::filter(TP_TN == max(TP_TN))

Prediction <- Prediction %>% mutate(Binary = ifelse(layer >= 
    EvalThres$Threshold, "Presencia", "Ausencia"))

• Tenemos que hacer crossvalidation
• Probar distintos valores de regularización
• Algunos candidatos de Features
• ENMeval

Results <- ENMevaluate(occ = N_pumilio[, c(3, 2)], env = Bioclimatic, 
    RMvalues = c(0.75, 1, 1.25), n.bg = 5000, method = "randomkfold", 
    overlap = F, kfolds = 5, bin.output = T, fc = c("L", "LQ", 
        "LQP"), rasterPreds = T)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |======================================================================| 100%

Models <- Results@results
Models$ID <- 1:nrow(Models)
Models <- Models %>% arrange(AICc)

# Indice del mejor modelo
Models$ID[1]

## [1] 6

# Mejor modelo
BestModels <- Results@models[[Models$ID[1]]]
## Prediccion del mejor modelo
Prediction <- predict(Bioclimatic, BestModels, type = "cloglog") %>% 
    as("SpatialPixelsDataFrame") %>% as.data.frame()