In biomodhub/biomod2: Ensemble Platform for Species Distribution Modeling

knitr::opts_chunk$set(echo = TRUE, comment = NA)
options(width = 100)
library(biomod2)
library(terra)
library(data.table)
library(kableExtra)
library(DT)
library(patchwork)

bm.calib = NULL
hasEvalMod = hasEvalEns = FALSE
hasVarImpMod = hasVarImpEns = FALSE
hasBmFiles = !is.na(bm.files)
hasBmMod = !is.na(bm.mod)
hasBmEns = !is.na(bm.ens)
hasBmForm = !is.na(bm.form)
ind.mod = ind.ens = ind.proj = NA
color1 = params.color$color1
color2 = params.color$color2
color3 = params.color$color3

bm.calib = get_calib_lines(bm.mod)

hasPA = inherits(bm.form, "BIOMOD.formated.data.PA")
hasDataEval = ifelse("has.data.eval" %in% slotNames(bm.form), bm.form@has.data.eval, FALSE)

# form.summ = summary(bm.form, calib.lines = bm.calib) ## TODO CORRECT
form.summ = data.frame("dataset" = "initial",
                       "run" = NA,
                       "PA" = NA,
                       "Presences" = length(which(bm.form@data.species == 1)),
                       "True_Absences" = length(which(bm.form@data.species == 0)),
                       "Pseudo_Absences" = 0,
                       "Undefined" = length(which(is.na(bm.form@data.species))))

nb.1 = form.summ$Presences[which(form.summ$dataset == "initial")]
nb.0 = form.summ$True_Absences[which(form.summ$dataset == "initial")]
nb.NA = form.summ$Pseudo_Absences[which(form.summ$dataset == "calibration" &
                                          form.summ$run == "allRun" &
                                          form.summ$PA == "PA1")]

# ind.mod = which(bm.files$file == sub("[.]/", "", bm.mod@link))
ind.mod = which(bm.files$file == sub(paste0(".*", sp.name, "/"), paste0(sp.name, "/"), bm.mod@link))
ind.proj = which(bm.files$type == "projection" &
                   bm.files$refer_to == bm.files$refer_to[ind.mod])
name.proj = unique(bm.files$ID[ind.proj])

ind.ens = which(bm.files$type == "models" &
                        bm.files$level == "ensemble" &
                        bm.files$refer_to == bm.files$refer_to[ind.mod])

eval.mod = get_evaluations(bm.mod)
hasEvalMod = !is.null(eval.mod)

varImp.mod = get_variables_importance(bm.mod)
hasVarImpMod = !is.null(varImp.mod)
if (hasVarImpMod) {
  varImp.mod$var.imp = round(varImp.mod$var.imp, 3)
}

eval.ens = get_evaluations(bm.ens)
hasEvalEns = !is.null(eval.ens)

varImp.ens = get_variables_importance(bm.ens)
hasVarImpEns = !is.null(varImp.ens)
if (hasVarImpEns) {
  varImp.ens$var.imp = round(varImp.ens$var.imp, 3)
}

# name.var = get_formal_data(bm.mod, subinfo = "expl.var.names")
name.models.mod = get_built_models(bm.mod)
name.algo = sapply(name.models.mod, function(x) strsplit(x, "_")[[1]][4])
name.algo = sort(unique(name.algo))

name.models.ens = get_built_models(bm.ens)
name.ens = sapply(name.models.ens, function(x) strsplit(x, "_")[[1]][2])
name.ens = sort(unique(name.ens))

tab = knitr::kable(bm.files[, c("file", "path", "type", "level", "ID")], format = "html")
tab = kable_styling(tab, bootstrap_options = "condensed")
if (any(!is.na(c(ind.mod, ind.ens)))) {
  row_spec(tab, c(ind.mod, ind.ens), bold = TRUE, color = color1)
}

Formated data

The r sp.name species dataset contains r nb.1 presences and r nb.0 true absences.

cat('
Pseudo-absences have been drawn from the'
, form.summ$Undefined[which(form.summ$dataset == "initial")]
, ' available background points.
')

print(bm.form)

tab = knitr::kable(form.summ, format = "html")
kable_styling(tab, bootstrap_options = "hover")

form.plot = plot(bm.form, calib.lines = bm.calib, plot.type = "raster")

Single models

print(bm.mod)

Modeling options

print(get_options(bm.mod))

Evaluation

rmCol = c("full.name")
if (!hasDataEval) { rmCol = c("full.name", "evaluation") }
DT::datatable(eval.mod[, -which(colnames(eval.mod) %in% rmCol)])

pp1 = bm_PlotEvalMean(bm.out = bm.mod, dataset = "calibration", do.plot = FALSE)
pp2 = bm_PlotEvalMean(bm.out = bm.mod, dataset = "validation", do.plot = FALSE)
pp3 = bm_PlotEvalMean(bm.out = bm.mod, dataset = "evaluation", do.plot = FALSE)
pp1$plot + pp2$plot + pp3$plot

pp1 = bm_PlotEvalBoxplot(bm.out = bm.mod, dataset = "calibration"
                         , group.by = c("algo", "run"), do.plot = FALSE)
pp2 = bm_PlotEvalBoxplot(bm.out = bm.mod, dataset = "validation"
                         , group.by = c("algo", "run"), do.plot = FALSE)
pp3 = bm_PlotEvalBoxplot(bm.out = bm.mod, dataset = "evaluation"
                         , group.by = c("algo", "run"), do.plot = FALSE)
pp1$plot + pp2$plot + pp3$plot

Variables' importance

rmCol = c("full.name")
DT::datatable(varImp.mod[, -which(colnames(varImp.mod) %in% rmCol)])

pp = bm_PlotVarImpBoxplot(bm.out = bm.mod
                          , group.by = c("algo", "expl.var", "expl.var")
                          , do.plot = FALSE)
pp$plot

Response curves

if (length(name.models.mod) > 9) {
  # mod.ind = seq(1, length(name.models.mod), 9)
  # mod.ind = c(mod.ind, length(name.models.mod))
  # mod.ind = unique(mod.ind)
  # mod.sel = lapply(1:(length(mod.ind)-1), function(x) {
  #   name.models.mod[mod.ind[x]:mod.ind[x+1]]
  # })

  mod.sel = lapply(name.algo, function(x) {
    grep(x, name.models.mod, value = TRUE)
  })
} else {
  mod.sel = list("all")
}

for (mod.i in mod.sel) {
  pp = bm_PlotResponseCurves(bm.out = bm.mod
                             , fixed.var = "mean"
                             , models.chosen = mod.i
                             , do.bivariate = FALSE
                             , do.plot = FALSE
                             , do.progress = FALSE)
  pp = pp$plot + ylim(0, 1) + guides(color = guide_legend(ncol = 3))
  print(pp)
}

# require(foreach)
# require(confintr)
# require(ggplot2)
# 
# # tmp = resp_models$tab
# tmp = TAB
# tmp$pred.name = as.character(tmp$pred.name)
# tmp$PA = sapply(tmp$pred.name, function(x) strsplit(x, "_")[[1]][2])
# tmp$RUN = sapply(tmp$pred.name, function(x) strsplit(x, "_")[[1]][3])
# tmp$algo = sapply(tmp$pred.name, function(x) strsplit(x, "_")[[1]][4])
# tmp.split = split(tmp, list(tmp$expl.name, tmp$algo, tmp$id))
# tmp.mean = foreach(tabi = tmp.split, .combine = "rbind") %do%
#   {
#     res = unique(tabi[, c("id", "expl.name", "algo", "expl.val")])
#     if (length(unique(tabi$pred.val)) == 1) {
#       res$pred.val = res$conf.inf = res$conf.sup = mean(tabi$pred.val)
#     } else {
#       res$pred.val = mean(tabi$pred.val)
#       res$conf.inf = ci_mean(tabi$pred.val)$interval[1]
#       res$conf.sup = ci_mean(tabi$pred.val)$interval[2]
#     }
#     return(res)
#   }
# head(tmp.mean)
# 
# ggplot(tmp.mean, aes_string(x = "expl.val", y = "pred.val", group = "algo")) + 
#   geom_ribbon(aes(ymin = conf.inf, ymax = conf.sup), alpha = 0.5, fill = "grey60") +
#   geom_line(aes(color = algo)) + 
#   facet_wrap("expl.name", scales = "free_x") + 
#   xlab("") + 
#   ylab("") + 
#   ylim(0, 1) +
#   theme(legend.title = element_blank(), 
#         legend.key = element_rect(fill = "white"), legend.position = "bottom", 
#         axis.text.x = element_text(angle = 45, hjust = 1))

Ensemble models

print(bm.ens)

Evaluation

rmCol = c("full.name")
if (!hasDataEval) { rmCol = c("full.name", "evaluation") }
tab = knitr::kable(eval.ens[, -which(colnames(eval.ens) %in% rmCol)], format = "html")
kable_styling(tab, bootstrap_options = "hover")

pp1 = bm_PlotEvalMean(bm.out = bm.ens, dataset = "calibration", do.plot = FALSE)
pp2 = bm_PlotEvalMean(bm.out = bm.ens, dataset = "validation", do.plot = FALSE)
pp3 = bm_PlotEvalMean(bm.out = bm.ens, dataset = "evaluation", do.plot = FALSE)
pp1$plot + pp2$plot + pp3$plot

pp1 = bm_PlotEvalBoxplot(bm.out = bm.ens, dataset = "calibration"
                         , group.by = c("algo", "algo"), do.plot = FALSE)
pp2 = bm_PlotEvalBoxplot(bm.out = bm.ens, dataset = "validation"
                         , group.by = c("algo", "algo"), do.plot = FALSE)
pp3 = bm_PlotEvalBoxplot(bm.out = bm.ens, dataset = "evaluation"
                         , group.by = c("algo", "algo"), do.plot = FALSE)
pp1$plot + pp2$plot + pp3$plot

Variables' importance

rmCol = c("full.name")
DT::datatable(varImp.ens[, -which(colnames(varImp.ens) %in% rmCol)])

pp = bm_PlotVarImpBoxplot(bm.out = bm.ens
                          , group.by = c("algo", "expl.var", "expl.var")
                          , do.plot = FALSE)
pp$plot

Response curves

if (any(grepl("EMcv", name.models.ens))) {
  mod.sel = list(name.models.ens[!grepl("EMcv", name.models.ens)]
                 , name.models.ens[grepl("EMcv", name.models.ens)])
} else {
  mod.sel = list("all")
}

for (mod.i in mod.sel) {
  pp = bm_PlotResponseCurves(bm.out = bm.ens
                             , fixed.var = "mean"
                             , models.chosen = mod.i
                             , do.bivariate = FALSE
                             , do.plot = FALSE
                             , do.progress = FALSE)
  pp = pp$plot + guides(color = guide_legend(ncol = 1))
  print(pp)
}

Projections

tab = knitr::kable(bm.files[, c("file", "path", "type", "level", "ID")], format = "html")
tab = kable_styling(tab, bootstrap_options = "condensed")
if (any(!is.na(ind.proj))) {
  row_spec(tab, ind.proj, bold = TRUE, color = color3)
}

for (proj.i in name.proj) {
  ind.proj.i = intersect(ind.proj, which(bm.files$ID == proj.i))
  for (j in ind.proj.i) {
    bm.proj = get(load(file = bm.files$file[j]))
    if (inherits(try(bm.proj@data.type), "try-error")) {
      bm.proj@data.type <- "binary"
    }
    print(bm.proj)
    plot(bm.proj)
    ## + clamping mask ??
  }
}
rm(list = basename(bm.files$file))