Here are presented, in a full and complete example, all main functions (starting with BIOMOD_[...]) of biomod2.
library(biomod2)
library(terra)
# Load species occurrences (6 species available)
data("DataSpecies")
head(DataSpecies)
# Select the name of the studied species
myRespName <- 'GuloGulo'
# Get corresponding presence/absence data
myResp <- as.numeric(DataSpecies[, myRespName])
# Get corresponding XY coordinates
myRespXY <- DataSpecies[, c('X_WGS84', 'Y_WGS84')]
# Load environmental variables extracted from BIOCLIM (bio_3, bio_4, bio_7, bio_11 & bio_12)
data("bioclim_current")
myExpl <- rast(bioclim_current)
# Format Data with true absences
myBiomodData <- BIOMOD_FormatingData(resp.var = myResp,
expl.var = myExpl,
resp.xy = myRespXY,
resp.name = myRespName)
myBiomodData
plot(myBiomodData)
Single or multiple set of pseudo-absences can be selected with the BIOMOD_FormatingData function, which calls the bm_PseudoAbsences function to do so. More examples are presented on the Secundary functions webpage.
# # Transform true absences into potential pseudo-absences
# myResp.PA <- ifelse(myResp == 1, 1, NA)
#
# # Format Data with pseudo-absences : random method
# myBiomodData.r <- BIOMOD_FormatingData(resp.var = myResp.PA,
# expl.var = myExpl,
# resp.xy = myRespXY,
# resp.name = myRespName,
# PA.nb.rep = 4,
# PA.nb.absences = 1000,
# PA.strategy = 'random')
#
# myBiomodData.r
# plot(myBiomodData.r)
# # Select multiple sets of pseudo-absences
#
# # Transform true absences into potential pseudo-absences
# myResp.PA <- ifelse(myResp == 1, 1, NA)
#
# # Format Data with pseudo-absences : random method
# myBiomodData.multi <- BIOMOD_FormatingData(resp.var = myResp.PA,
# expl.var = myExpl,
# resp.xy = myRespXY,
# resp.name = myRespName,
# PA.nb.rep = 4,
# PA.nb.absences = c(1000, 500, 500, 200),
# PA.strategy = 'random')
# myBiomodData.multi
# summary(myBiomodData.multi)
# plot(myBiomodData.multi)
# Print default modeling options
bm_DefaultModelingOptions()
# Create default modeling options
myBiomodOptions <- BIOMOD_ModelingOptions()
myBiomodOptions
# # Part (or totality) of the print can be copied and customized
# # Below is an example to compute quadratic GLM and select best model with 'BIC' criterium
# myBiomodOptions <- BIOMOD_ModelingOptions(
# GLM = list(type = 'quadratic',
# interaction.level = 0,
# myFormula = NULL,
# test = 'BIC',
# family = 'binomial',
# control = glm.control(epsilon = 1e-08,
# maxit = 1000,
# trace = FALSE)))
# myBiomodOptions
#
# # It is also possible to give a specific GLM formula
# myForm <- 'Sp277 ~ bio3 + log(bio10) + poly(bio16, 2) + bio19 + bio3:bio19'
# myBiomodOptions <- BIOMOD_ModelingOptions(GLM = list(myFormula = formula(myForm)))
# myBiomodOptions
### Model parameters can also be automatically tuned for your specific
### dataset with an optimization algorithm. The tuning can however
### be quite long. Duration for tuning all models sequentially
### with default optimization settings :
### on 1 x 2.5 GHz processor: approx. 45 min tuning all models
### on 8 x 2.5 GHz processor: approx. 15 min tuning all models
#
# # library(doParallel)
# # cl <- makeCluster(8)
# # doParallel::registerDoParallel(cl)
#
# time.seq <- system.time(
# bm.tuning <- BIOMOD_Tuning(bm.format = myBiomodData, ME.env = myExpl, ME.n.bg = ncell(myExpl))
# )
#
# # stopCluster(cl)
#
# plot(bm.tuning$tune.CTA.rpart)
# plot(bm.tuning$tune.CTA.rpart2)
# plot(bm.tuning$tune.RF)
# plot(bm.tuning$tune.ANN)
# plot(bm.tuning$tune.MARS)
# plot(bm.tuning$tune.FDA)
# plot(bm.tuning$tune.GBM)
# plot(bm.tuning$tune.GAM)
#
# # Get tuned modeling options
# myBiomodOptions <- bm.tuning$models.options
Several cross-validation methods are available and can be selected with the BIOMOD_Modeling function, which calls the bm_CrossValidation function to do so. More examples are presented on the Secundary functions webpage.
# # k-fold selection
# cv.k <- bm_CrossValidation(bm.format = myBiomodData,
# strategy = "kfold",
# nb.rep = 2,
# k = 3)
#
# # stratified selection (geographic)
# cv.s <- bm_CrossValidation(bm.format = myBiomodData,
# strategy = "strat",
# k = 2,
# balance = "presences",
# strat = "x")
# head(cv.k)
# head(cv.s)
# Model single models
myBiomodModelOut <- BIOMOD_Modeling(bm.format = myBiomodData,
bm.options = myBiomodOptions,
modeling.id = 'AllModels',
CV.strategy = 'random',
CV.nb.rep = 2,
CV.perc = 0.8,
var.import = 3,
metric.eval = c('TSS','ROC'))
# seed.val = 123)
# nb.cpu = 8)
myBiomodModelOut
# Get evaluation scores & variables importance
get_evaluations(myBiomodModelOut)
get_variables_importance(myBiomodModelOut)
# Represent evaluation scores & variables importance
bm_PlotEvalMean(bm.out = myBiomodModelOut)
bm_PlotEvalBoxplot(bm.out = myBiomodModelOut, group.by = c('algo', 'algo'))
bm_PlotEvalBoxplot(bm.out = myBiomodModelOut, group.by = c('algo', 'run'))
bm_PlotVarImpBoxplot(bm.out = myBiomodModelOut, group.by = c('expl.var', 'algo', 'algo'))
bm_PlotVarImpBoxplot(bm.out = myBiomodModelOut, group.by = c('expl.var', 'algo', 'run'))
bm_PlotVarImpBoxplot(bm.out = myBiomodModelOut, group.by = c('algo', 'expl.var', 'run'))
# Represent response curves
bm_PlotResponseCurves(bm.out = myBiomodModelOut,
models.chosen = get_built_models(myBiomodModelOut)[c(1:3, 12:14)],
fixed.var = 'median')
bm_PlotResponseCurves(bm.out = myBiomodModelOut,
models.chosen = get_built_models(myBiomodModelOut)[c(1:3, 12:14)],
fixed.var = 'min')
bm_PlotResponseCurves(bm.out = myBiomodModelOut,
models.chosen = get_built_models(myBiomodModelOut)[3],
fixed.var = 'median',
do.bivariate = TRUE)
# Model ensemble models
myBiomodEM <- BIOMOD_EnsembleModeling(bm.mod = myBiomodModelOut,
models.chosen = 'all',
em.by = 'all',
em.algo = c('EMmean', 'EMcv', 'EMci', 'EMmedian', 'EMca', 'EMwmean'),
metric.select = c('TSS'),
metric.select.thresh = c(0.7),
metric.eval = c('TSS', 'ROC'),
var.import = 3,
EMci.alpha = 0.05,
EMwmean.decay = 'proportional')
myBiomodEM
# Get evaluation scores & variables importance
get_evaluations(myBiomodEM)
get_variables_importance(myBiomodEM)
# Represent evaluation scores & variables importance
bm_PlotEvalMean(bm.out = myBiomodEM, group.by = 'full.name')
bm_PlotEvalBoxplot(bm.out = myBiomodEM, group.by = c('full.name', 'full.name'))
bm_PlotVarImpBoxplot(bm.out = myBiomodEM, group.by = c('expl.var', 'full.name', 'full.name'))
bm_PlotVarImpBoxplot(bm.out = myBiomodEM, group.by = c('expl.var', 'algo', 'merged.by.run'))
bm_PlotVarImpBoxplot(bm.out = myBiomodEM, group.by = c('algo', 'expl.var', 'merged.by.run'))
# Represent response curves
bm_PlotResponseCurves(bm.out = myBiomodEM,
models.chosen = get_built_models(myBiomodEM)[c(1, 6, 7)],
fixed.var = 'median')
bm_PlotResponseCurves(bm.out = myBiomodEM,
models.chosen = get_built_models(myBiomodEM)[c(1, 6, 7)],
fixed.var = 'min')
bm_PlotResponseCurves(bm.out = myBiomodEM,
models.chosen = get_built_models(myBiomodEM)[7],
fixed.var = 'median',
do.bivariate = TRUE)
# Evaluate models with Boyce index and MPA
myBiomodPO <- BIOMOD_PresenceOnly(bm.mod = myBiomodModelOut,
bm.em = myBiomodEM)
myBiomodPO
# Evaluate models with Boyce index and MPA (using background data)
myBiomodPO <- BIOMOD_PresenceOnly(bm.mod = myBiomodModelOut,
bm.em = myBiomodEM,
bg.env = values(myExpl))
myBiomodPO
# Project single models
myBiomodProj <- BIOMOD_Projection(bm.mod = myBiomodModelOut,
proj.name = 'Current',
new.env = myExpl,
models.chosen = 'all',
metric.binary = 'all',
metric.filter = 'all',
build.clamping.mask = TRUE)
myBiomodProj
plot(myBiomodProj)
# Project ensemble models (from single projections)
myBiomodEMProj <- BIOMOD_EnsembleForecasting(bm.em = myBiomodEM,
bm.proj = myBiomodProj,
models.chosen = 'all',
metric.binary = 'all',
metric.filter = 'all')
# Project ensemble models (building single projections)
myBiomodEMProj <- BIOMOD_EnsembleForecasting(bm.em = myBiomodEM,
proj.name = 'CurrentEM',
new.env = myExpl,
models.chosen = 'all',
metric.binary = 'all',
metric.filter = 'all')
myBiomodEMProj
plot(myBiomodEMProj)
# Load environmental variables extracted from BIOCLIM (bio_3, bio_4, bio_7, bio_11 & bio_12)
data("bioclim_future")
myExplFuture = rast(bioclim_future)
# Project onto future conditions
myBiomodProjectionFuture <- BIOMOD_Projection(bm.mod = myBiomodModelOut,
proj.name = 'Future',
new.env = myExplFuture,
models.chosen = 'all',
metric.binary = 'TSS',
build.clamping.mask = TRUE)
# Load current and future binary projections
CurrentProj <- get_predictions(myBiomodProj, metric.binary = "TSS")
FutureProj <- get_predictions(myBiomodProjectionFuture, metric.binary = "TSS")
# Compute differences
myBiomodRangeSize <- BIOMOD_RangeSize(proj.current = CurrentProj,
proj.future = FutureProj)
myBiomodRangeSize$Compt.By.Models
plot(myBiomodRangeSize$Diff.By.Pixel)
# Represent main results
gg = bm_PlotRangeSize(bm.range = myBiomodRangeSize,
do.count = TRUE,
do.perc = TRUE,
do.maps = TRUE,
do.mean = TRUE,
do.plot = TRUE,
row.names = c("Species", "Dataset", "Run", "Algo"))