Abstract

This script aims to show you how you can use outputs from response.plot2() function to produce a fully customizable response plot within ggplot2 framework.

This example is designed with Vulpes vulpes distribution from biomod2 internal data.

The script

# load reauired libraries
suppressPackageStartupMessages(library(biomod2))
suppressPackageStartupMessages(library(raster))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))


# species occurrences
sp.dat <- 
  read.csv(
    system.file(
      "external/species/mammals_table.csv",
      package="biomod2"
    ), 
    row.names = 1,
    stringsAsFactors = FALSE)

# the name of studied species
sp.name <- 'VulpesVulpes'

# the presence/absences data for our species 
sp.resp <- sp.dat %>% pull(sp.name) %>% as.numeric()

# the XY coordinates of species data
sp.xy <- sp.dat %>% select(X_WGS84, Y_WGS84)


# Environmental variables extracted from BIOCLIM (bio_3, bio_4, bio_7, bio_11 & bio_12)
sp.expl <- 
  raster::stack( 
    system.file(
      "external/bioclim/current/bio3.grd", 
      package="biomod2"
    ),
    system.file( 
      "external/bioclim/current/bio4.grd", 
      package="biomod2"
    ),
    system.file(
      "external/bioclim/current/bio7.grd", 
      package="biomod2"
    ),
    system.file( 
      "external/bioclim/current/bio11.grd", 
      package="biomod2"
    ),
    system.file(
      "external/bioclim/current/bio12.grd", 
      package="biomod2"
    )
  )

# 1. Formatting Data
bm.dat <- 
  BIOMOD_FormatingData(
    resp.var = sp.resp,
    expl.var = sp.expl,
    resp.xy = sp.xy,
    resp.name = sp.name
  )
## 
## -=-=-=-=-=-=-=-=-=-=-=-= VulpesVulpes Data Formating -=-=-=-=-=-=-=-=-=-=-=-=
## 
## > No pseudo absences selection !
##       ! No data has been set aside for modeling evaluation
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
# 2. Defining Models Options using default options.
bm.opt <- BIOMOD_ModelingOptions()

# 3. Doing Modelisation
bm.mod <- 
  BIOMOD_Modeling(
    bm.dat,
    models = c('RF','GLM'), 
    models.options = bm.opt, 
    NbRunEval = 2, 
    DataSplit = 80, 
    Prevalence = 0.5,
    VarImport = 0, 
    models.eval.meth = c('TSS'),
    SaveObj = TRUE,
    do.full.models = FALSE
  )
## 
## 
## Loading required library...
## 
## Checking Models arguments...
## 
## Creating suitable Workdir...
## 
##  > Automatic weights creation to rise a 0.5 prevalence
## 
## 
## -=-=-=-=-=-=-=-=-=-=-= VulpesVulpes Modeling Summary -=-=-=-=-=-=-=-=-=-=-=
## 
##  5  environmental variables ( bio3 bio4 bio7 bio11 bio12 )
## Number of evaluation repetitions : 2
## Models selected : RF GLM 
## 
## Total number of model runs : 4 
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## 
## 
## -=-=-=- Run :  VulpesVulpes_AllData 
## 
## 
## -=-=-=--=-=-=- VulpesVulpes_AllData_RUN1 
## 
## Model=Breiman and Cutler's random forests for classification and regression
##  Evaluating Model stuff...
## Model=GLM ( quadratic with no interaction )
##  Stepwise procedure using AIC criteria
##  selected formula : VulpesVulpes ~ bio7 + I(bio11^2) + I(bio4^2) + bio11 + I(bio3^2) + 
##     bio3 + I(bio7^2) + bio4 + I(bio12^2)
## <environment: 0x5616ee728920>
## 
##  Evaluating Model stuff...
## 
## -=-=-=--=-=-=- VulpesVulpes_AllData_RUN2 
## 
## Model=Breiman and Cutler's random forests for classification and regression
##  Evaluating Model stuff...
## Model=GLM ( quadratic with no interaction )
##  Stepwise procedure using AIC criteria
##  selected formula : VulpesVulpes ~ bio7 + I(bio11^2) + I(bio4^2) + bio11 + I(bio3^2) + 
##     bio3 + I(bio7^2) + bio4
## <environment: 0x5616f4fbf020>
## 
##  Evaluating Model stuff...
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
# 4. Plot response curves

# 4.1 Load the models for which we want to extract the predicted response curves
bm.mod.names <- 
  BIOMOD_LoadModels(
    bm.mod, 
    models = c('GLM', 'RF')
  )
# 4.2 get 2D response plots data
rp.dat <- 
  response.plot2(
    models  = bm.mod.names,
    Data = get_formal_data(bm.mod, 'expl.var'),
    show.variables = get_formal_data(bm.mod, 'expl.var.names'),
    data_species = get_formal_data(bm.mod, 'resp.var'),
    do.bivariate = FALSE,
    fixed.var.metric = 'median',
    col = c("blue", "red"),
    legend = FALSE,
    plot = FALSE
  )
# 4.3 use ggplot to display teh results
rp.dat %>%
  ## transform the pred.name to extract model, cv run and data info
  mutate(
    species = pred.name %>% strsplit('_') %>% sapply(function(x) x[1]),
    pa.dat = pred.name %>% strsplit('_') %>% sapply(function(x) x[2]),
    cv.rep = pred.name %>% strsplit('_') %>% sapply(function(x) x[3]),
    model = pred.name %>% strsplit('_') %>% sapply(function(x) x[4])
  ) %>%
  ggplot(
    aes(
      x = expl.val,
      y = pred.val,
      colour = model,
      group = pred.name
    )
  ) +
  geom_line(size = 1) +
  facet_wrap(~ expl.name, scales = 'free_x') + 
  labs(
    x = '',
    y = 'probability of occurence',
    colour = 'model type'
  ) + 
  scale_color_brewer(type = 'qual', palette = 4) +
  theme_minimal() +
  theme(
    legend.position = 'bottom'
  )