Analyze the residuals of the single models

Source: R/bm_ModelAnalysis.R

This function return several graphs to help analyse the single models.

bm_ModelAnalysis(
  bm.mod,
  models.chosen = "all",
  do.plot = TRUE,
  color.by = "full.name"
)

Arguments

bm.mod

a BIOMOD.models.out

models.chosen

a vector containing model names to be kept, must be either all or a sub-selection of model names that can be obtained with the get_built_models function

do.plot

a logical, print the plots or not

color.by

a character between "full.name", "species", "PA", "RUN" or "algo" to select the color parameter of the plots.

Value

A list with:

  1. a dataframe with the observation data, the fitted values and the residuals for all models

  2. a dataframe with the compilation of the R scores

  3. A to E plots (see details)

Details

All the plots will be made for all the models, independently to the different models assumptions. It is up to the user to interpret the graphs in the light of the model assumptions.

  1. Plot A: residuals ~ observations number. This plot helps to detect one or several outliers. The x-axis only helps to find the outlier number.

  2. Plot B and C: These are two representation of the distribution of the residuals. If your residuals must follow a normal distribution, the points should follow the black line in the Q-Q plot and present a gaussian distribution on the histogram.

  3. Plot D residuals ~ fitted values. This plot helps to detect an heteroscedasticity of the residuals.

  4. Plot E : a plot of the different Rsquared values available (Rsquared, Rsquared_aj). An big gap between the calibration and validation values can be the sign of an overfitting.

Author

Hélène Blancheteau

Examples

library(terra)

# Load species occurrences (6 species available)
data(DataSpecies)
head(DataSpecies)

# Select the name of the studied species
myRespName <- 'VulpesVulpes'

# 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 <- terra::rast(bioclim_current)

# Transform binary data as count data
poissonDistri <- rpois(sum(myResp), 10)
myResp[myResp == 1] <- poissonDistri

# ---------------------------------------------------------------
file.out <- paste0(myRespName, "/", myRespName, ".AllModels.models.out")
if (file.exists(file.out)) {
  myBiomodModelOut <- get(load(file.out))
} else {

  # Format Data with true absences
  myBiomodData <- BIOMOD_FormatingData(resp.name = myRespName,
                                       resp.var = myResp,
                                       resp.xy = myRespXY,
                                       expl.var = myExpl)

  # Model single models
  myBiomodModelOut <- BIOMOD_Modeling(bm.format = myBiomodData,
                                      modeling.id = 'AllModels',
                                      models = c('GAM', 'GLM', 'MARS'),
                                      CV.strategy = 'random',
                                      CV.nb.rep = 2,
                                      CV.perc = 0.8,
                                      OPT.strategy = 'bigboss',
                                      metric.eval = c('Rsquared', 'Rsquared_aj'),
                                      seed.val = 42)
}


# ---------------------------------------------------------------
# bm_ModelAnalysis
analysis <- bm_ModelAnalysis(myBiomodModelOut, color.by = "RUN")

plot(analysis$plotA)

unlink(myRespName, recursive = TRUE)