Running the analysis

Overview

In this module we will use the dataset obtained in module 2 to train the SDMs (occurrences + bioclimatic variables). We produce an SDM for P. krameri under the current climatic conditions and then we will project the model’s results in space (at global level) and in time taking into account the current and the future (“worst” IPCC Fifth Report scenario for the 2070, see also IPCC, 2013) climate conditions. The final training of the model will be performed using an ensemble forecasting approach based on the use of multiple modelling algorithms and the average of the respective predictions (see Thuiller et al., 2009).

SMD algorithms and the ensemble modelling

Species distribution models can be fitted using different techniques. The most used ones can be dived in two categories: regression-based techniques (generalized linear model, generalized additive model and multivariate adaptive regression splines) and machine learning techniques (maximum entropy, genetic algorithm for rule set production, gradient boosting machines, random forest, support vector machines).
You can choice to use one of these techniques or more than one. Different techniques often lead to results that can be also very different. This can happen because each method relies to different assumptions. To know the advantage and limitation of each technique that you will use is a good starting point for a proper interpretation of the output. However, several studies demonstrated that results of different model techniques can be so variable limiting their usefulness in a forecasting framework forecasting (Araújo & New, 2006). A possible solution is to combine multiple models in an ensemble model that present advantages over single models, particularly when the aim of SDM is prediction and forecasting (Araújo & New, 2006).

Prepare raster layers for the projection in space and time

First load the necessary libraries and data

> library(biomod2)
> library(rgeos)
> library(raster)
> library(maptools)
> library(gtools)
>data(wrld_simpl)

and the training dataset and the bioclimatic layers (from training area) saved in module 2

> load("krameri.training")
> vars_training<-stack(mixedsort(list.files("predictors/", "gri", full.names=T)))

Then we create a second set of varibles for project at a global scale under current climatic conditions

> r <- getData("worldclim", var = "bio", res = 2.5)
> r<-r[[which(names(r)%in%names(vars_training))]]

For reasons of computational time we reduced the spatial resolution of the bioclimatic layers used for the projections from 5km to 25km.

> vars_proj<-stack(aggregate(r, 5, progress="text"))

Finally, we create a third set of bioclimatic variables to project our models at the global scale under future climate conditions (IPCC, 2013; downloaded from Worldclim).
We download and unzip the files creating a new folder named “future”:

> download.file("http://biogeo.ucdavis.edu/data/climate/cmip5/2_5m/cc85bi70.zip", "future.zip")
> unzip("future.zip", exdir="future")

Load layers:

> future<-stack(mixedsort(list.files("future","tif",full.names=T)))
> names(future)<-paste("bio", 1:19, sep="")
> future<-future[[which(names(future)%in%names(vars_training))]]

and downscale to 25 Km:

> future<-stack(aggregate(future, 5, progress="text"))

Now we are ready to fit the models

Fit the models

In this section we will see, in 14 steps, how to fit our models using an ensemble forecasting approach based on the use of multiple modelling algorithms and the average of the respective predictions

Fit the models - Step 1/14

First of all we create the directory were the results will be stored

> dir.create("biomod.output")

and set the new working directory

> setwd("biomod.output")

To fit our models and the ensemble model we will use the procedure implemented in the R package “biomod2”. The first step is to format the input data. In the module 2 we prepared a dataset including three replicates (to take into account the coordinates uncertainty in datasets). We need to train a complete model from each of them. A detailed description of the procedure will be shown for the first replicate, with the remaining two runs executed within a for loop. We will use the first replicate (i=1) of the krameri dataset. These are only presence data (resp.var=rep(1, nrow(krameri[[i]])), the explanatory variables are in the vars_training object. We will provide coordinates of the first replicate (resp.xy=krameri[[i]]) Together with occurrences data we need absence data. In our case we do not have true absence so we will use pseudo-absence. We set one replicate (PA.nb.rep=1) of 1000 pseudoabsences (PA.nb.absences=1000), randomly placed (PA.strategy='random') in the study area identified in module 2. Regarding the choice of pseudoabsences, see Barbet-Massin et al. (2012).

> i=1
> mydata<-BIOMOD_FormatingData(resp.var=rep(1, nrow(krameri[[i]])),
+               expl.var=vars_training,
+               PA.nb.rep=1,
+               PA.nb.absences=1000,
+               PA.strategy='random',
+               resp.xy=krameri[[i]],
+               resp.name=paste("krameri.", i, sep=""))

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= krameri.1 Data Formating -=-=-=-=-=-=-=-=-=-=-=-=-

! No data has been set aside for modeling evaluation
> Pseudo Absences Selection checkings...
> random pseudo absences selection
> Pseudo absences are selected in explanatory variables
-=-=-=-==-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Fit the models - Step 2/14

have a look to the input data for the SDM training

> mydata

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.formated.data.PA' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

sp.name = krameri.1

143 presences, 0 true absences and 1000 undifined points in dataset


7 explanatory variables

1 Pseudo Absences dataset available ( PA1 ) with 1000 absences in each (true abs + pseudo abs)

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-==-=-==-=-==-=-==-

Fit the models - Step 3/14

The second step is to set the parameters, by using the BIOMOD_ModelingOptions() function provided by the biomod2 packages, for the modelling algorithms that are going to be used. Here, we will use four different models: generalized linear model (GLM), gradient boosting machine (GBM), generalized additive model (GAM) and maximum entropy (MAXENT):

> myBiomodOption<-BIOMOD_ModelingOptions(GLM=list(type='polynomial'),
+                                    GBM=list(n.trees=1000),
+                                    GAM=list(k=4),
+                                    MAXENT=list(path_to_maxent.jar="../",
+                                                  maximumiterations=1000))

If you want to have a look at the parameters of the four models (here we reported only those of the models chosen) just type:

> myBiomodOption

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.Model.Options' -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
GLM = list( type = 'polynomial',
    interaction.level = 0,
    myFormula = NULL,
    test = 'AIC',
    family = binomial(link = 'logit'),
     mustart = 0.5,
     control = glm.control(epsilon = 1e-08, maxit = 50, trace = FALSE) ),

GBM = list( distribution = 'bernoulli',
    n.trees = 1000,
    interaction.depth = 7,
    n.minobsinnode = 5,
     shrinkage = 0.001,
     bag.fraction = 0.5,
    train.fraction = 1,
    cv.folds = 3,
    keep.data = FALSE,
    verbose = FALSE,
    perf.method = 'cv'),

GAM = list( algo = 'GAM_mgcv',
    type = 's_smoother',
    k = 4,
    interaction.level = 0,
    myFormula = NULL,
    family = binomial(link = 'logit'),
    method = 'GCV.Cp',
    optimizer = c('outer','newton'),
    select = FALSE,
    knots = NULL,
    paraPen = NULL,
    control = list(nthreads = 1, irls.reg = 0, epsilon = 1e-07, maxit = 100, trace = FALSE
, mgcv.tol = 1e-07, mgcv.half = 15, rank.tol = 1.49011611938477e-08
, nlm = list(ndigit=7, gradtol=1e-06, stepmax=2, steptol=1e-04, iterlim=200, check.analyticals=0)
, optim = list(factr=1e+07), newton = list(conv.tol=1e-06, maxNstep=5, maxSstep=2, maxHalf=30, use.svd=0)
, outerPIsteps = 0, idLinksBases = TRUE, scalePenalty = TRUE, keepData = FALSE) ),

MAXENT = list( path_to_maxent.jar = '/Users/paolo/Dropbox/SDM',
    memory_allocated = 512,
    maximumiterations = 1000,
    visible = FALSE,
    linear = TRUE,
    quadratic = TRUE,
    product = TRUE,
    threshold = TRUE,
    hinge = TRUE,
    lq2lqptthreshold = 80,
    l2lqthreshold = 10,
    hingethreshold = 15,
    beta_threshold = -1,
    beta_categorical = -1,
    beta_lqp = -1,
    beta_hinge = -1,
    defaultprevalence = 0.5)
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

Fit the models - Step 4/14

Now we are ready to compute the models. By using the BIOMOD_Modeling() function we tell to biomod2 which models we want to use, the percentage of the dataset used to train the model and the evaluation metrics used for the single models. Here we will use 70% of the occurrence data (DataSplit=70) to train the models while the remaining 30% will be use for validation. This process was repeated 3 times (NbRunEval=3). The ROC and TSS were the metrics used to evaluate model performances:

> mymodel<-BIOMOD_Modeling(data=mydata,
+     models=c('GLM','GAM','GBM','MAXENT'),
+     models.options=myBiomodOption,
+      NbRunEval=3,
+     DataSplit=70,
+      VarImport=2,
+     models.eval.meth=c('ROC', 'TSS'),
+      do.full.models=F)

Loading required library...

Checking Models arguments...

Creating suitable Workdir...

    ! Weights where automaticly defined for krameri.1_PA1 to rise a 0.5 prevalence !


-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= krameri.1 Modeling Summary -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

7 environmental variables ( bio1 bio2 bio7 bio12 bio15 bio17 bio19 )
Number of evaluation repetitions : 3
Models selected : GLM GAM GBM MAXENT

Total number of model runs : 12

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

-=-=-=- Run : krameri.1_PA1

-=-=-=--=-=-=- krameri.1_PA1_RUN1

Model=GLM polynomial with no interaction
    Stepwise procedure using AIC criteria
    selected formula : krameri.1 ~ bio1 + I(bio2^3) + bio7 + I(bio1^3) + bio15 + I(bio15^3) +
    I(bio15^2) + I(bio12^2)

    Evaluating Model stuff...
    Evaluating Predictor Contributions...

Model=GAM
    GAM_mgcv algorithm chosen

    Automatic formula generation...
    > GAM (mgcv) modelling...
    Evaluating Model stuff...
    Evaluating Predictor Contributions...

Model=Generalised Boosting Regression
    1000 maximum different trees and 3 Fold Cross-Validation
    Evaluating Model stuff...
    Evaluating Predictor Contributions...

Model=MAXENT
    Creating Maxent Temp Proj Data..
Running Maxent...
Getting predictions...
    Removing Maxent Temp Data..
    Evaluating Model stuff...
    Evaluating Predictor Contributions...
…

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Fit the models - Step 5/14

Check the summary

> mymodel

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= BIOMOD.models.out -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Modeling id : 1415915636

Species modeled : krameri.1

Considered variables : bio1 bio2 bio7 bio12 bio15 bio17 bio19


Computed Models : krameri.1_PA1_RUN1_GLM krameri.1_PA1_RUN1_GAM krameri.1_PA1_RUN1_GBM
krameri.1_PA1_RUN1_MAXENT krameri.1_PA1_RUN2_GLM krameri.1_PA1_RUN2_GAM krameri.1_PA1_RUN2_GBM
krameri.1_PA1_RUN2_MAXENT krameri.1_PA1_RUN3_GLM krameri.1_PA1_RUN3_GAM krameri.1_PA1_RUN3_GBM
krameri.1_PA1_RUN3_MAXENT


Failed Models : none

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

Fit the models - Step 6/14

As you can see in the output we run 12 different models (4 models run 3 times each). All work correctly done without errors. If you anre interested in you can print out the variable importance for each model in each replicate:

> get_variables_importance(mymodel)

, , RUN1, PA1

, , RUN2, PA1

, , RUN3, PA1

Fit the models - Step 7/14

Now we check the result in the training area by projecting the models.
To save time we fist downscale the bioclimatic variables in the training area to 5Km:

> vars_training_5km<-stack(aggregate(vars_training, 5, progress="text"))

then project:

> myproj<-BIOMOD_Projection(modeling.output=mymodel,
+      new.env= vars_training_5km,
+      proj.name="current",
+      selected.models= "all",
+      do.stack=T)

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Do Models Projections -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

    > Building clamping mask

    > Projecting krameri.1_PA1_RUN1_GLM ...
    > Projecting krameri.1_PA1_RUN1_GAM ...
    > Projecting krameri.1_PA1_RUN1_GBM ...
    > Projecting krameri.1_PA1_RUN1_MAXENT ...
    > Projecting krameri.1_PA1_RUN2_GLM ...
    > Projecting krameri.1_PA1_RUN2_GAM ...
    > Projecting krameri.1_PA1_RUN2_GBM ...
    > Projecting krameri.1_PA1_RUN2_MAXENT ...
    > Projecting krameri.1_PA1_RUN3_GLM ...
    > Projecting krameri.1_PA1_RUN3_GAM ...
    > Projecting krameri.1_PA1_RUN3_GBM ...
    > Projecting krameri.1_PA1_RUN3_MAXENT ...
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

Fit the models - Step 8/14

and finally plot 4 maps (only from the first run):

> source("../script/palette_MAXENT.R")
> map.current<-stack("krameri.1/proj_current/proj_current_krameri.1.gri")
> plot(map.current[[c(1,2,3,4)]], col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]])

The figure report the projections for the four models in the training area. As you can see the four models are similar in identify areas with higher and lower probability of presence but the probability associated to these areas can differ largely among models.

Fit the models - Step 9/14

> myproj<-BIOMOD_Projection(modeling.output=mymodel,
+             new.env=vars_proj,
+             proj.name="current",
+             selected.models=”all”,
+             do.stack=T)

we plot the maps:

> source("../script/palette_MAXENT.R")
> map.current<-stack("krameri.1/proj_current/proj_current_krameri.1.gri")
> plot(map.current, col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]])

As you can see in this figure there is a large variability of the projections among different models but also among the 3 runs of each model. This justify the use of an ensemble model for projection in space and time.
To do this we will use the function BIOMOD_EnsembleModeling().

Fit the models - Step 10/14

It is possible to specify if we want to use all the models (em.by="all") or just some of these, which metrics we want to use to evaluate performance (“eval.metric” argument) and to set a threshold (“eval.metric.quality.threshold” argument) to filter the single models according to their performance (see also Swets, 1988 and Allouche et al., 2006). Furthermore, we specified which metrics must be used to evaluate the ensemble model (“models.eval.meth” argument):

> myEM<-BIOMOD_EnsembleModeling(modeling.output=mymodel,
+             em.by="all",
+             eval.metric=c('ROC','TSS'),
+             eval.metric.quality.threshold=c(0.7,0.4),
+             models.eval.meth = c('ROC', 'TSS'),
+             prob.mean=F,
            + prob.mean.weight=T)

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Build Ensemble Models -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

> Evaluation & Weighting methods summary :
ROC over 0.7
TSS over 0.4


> TotalConsensus ensemble modeling
! Models projections for whole zonation required...
> Projecting krameri.1_PA1_RUN1_GLM ...
> Projecting krameri.1_PA1_RUN1_GAM ...
> Projecting krameri.1_PA1_RUN1_GBM ...
> Projecting krameri.1_PA1_RUN1_MAXENT ...
> Projecting krameri.1_PA1_RUN2_GLM ...
> Projecting krameri.1_PA1_RUN2_GAM ...
> Projecting krameri.1_PA1_RUN2_GBM ...
> Projecting krameri.1_PA1_RUN2_MAXENT ...
> Projecting krameri.1_PA1_RUN3_GLM ...
> Projecting krameri.1_PA1_RUN3_GAM ...
> Projecting krameri.1_PA1_RUN3_GBM ...
> Projecting krameri.1_PA1_RUN3_MAXENT ...

> Prababilities wegthing mean...
original models scores = 0.73 0.724 0.715 0.702 0.762 0.727 0.75 0.751 0.708 0.722 0.711 0.708
final models weights = 0.084 0.083 0.082 0.081 0.087 0.083 0.086 0.086 0.081 0.083 0.082 0.081
Evaluating Model stuff...
> Prababilities wegthing mean...
original models scores = 0.427 0.401 0.448 0.471 0.434
final models weights = 0.196 0.184 0.205 0.216 0.199
Evaluating Model stuff...
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

From the output we understand that all the 12 models of the first replicate passed the threshold and will be used to obtain the ensemble model. It is likely that in other replicates only a part of the models will be included in the ensemble model.

Fit the models - Step 11/14

We create a probability map for the ensemble model under current climatic conditions:

> myEF<-BIOMOD_EnsembleForecasting(EM.output=myEM,
+ projection.output=myproj,
+ total.consensus=T,
+ binary.meth=c("TSS"))

a projection under future climatic conditions:

> myproj<-BIOMOD_Projection(modeling.output=mymodel,
+ new.env=future,
+ proj.name="future",
+ selected.models=”all”,
+ do.stack=T)

and the ensemble model map :

> myEF<-BIOMOD_EnsembleForecasting(EM.output=myEM,
+ projection.output=myproj,
+ total.consensus=T,
+ binary.meth=c("TSS"))

Fit the models - Step 12/14

Let’s plot the resulting maps for the ensemble models obtained for projection at global scale under current and future climatic conditions:

> source("../script/palette_MAXENT.R")
> map.current<-stack("krameri.1/proj_current/proj_current_krameri.1_ensemble.gri")
> plot(map.current[[2]], col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]])

Fit the models - Step 13/14

> map.future<-stack("krameri.1/proj_future/proj_future_krameri.1_ensemble.gri")
> plot(map.future[[2]], col=palette_MAXENT(0, 1000, 100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]])

These two maps are the ensemble SDMs for current and future climatic conditions obtained from the first replicates. To obtain our final model we need to fit an ensemble SDM models using all the three replicates.

Fit the models - Step 14/14

So now we will use the other two replicates using the same parameters that we showed for the first replicate.

> for(i in 2:length(krameri))
> {
 mydata<-BIOMOD_FormatingData(resp.var=rep(1, nrow(krameri[[i]])), <%6%5%> expl.var=vars_training,
PA.nb.rep=1,
PA.nb.absences=1000,
PA.strategy='random',
resp.xy=krameri[[i]],
resp.name=paste("krameri.", i, sep=""))

 mymodel<-BIOMOD_Modeling(data=mydata, <%6%5%> models=c('GLM','GAM','GBM','MAXENT'),
models.options=myBiomodOption,
NbRunEval=3,
DataSplit=70,
VarImport=2,
models.eval.meth=c('ROC', 'TSS'),
do.full.models=F)

 myEM<-BIOMOD_EnsembleModeling(modeling.output=mymodel,<%6%5%> em.by="all",
eval.metric=c('ROC','TSS'),
eval.metric.quality.threshold=c(0.7,0.4),
models.eval.meth = c('ROC', 'TSS'),
prob.mean=F,
prob.mean.weight=T)

 ## to save time we select for projections only those models that passed the threshold

 sel.mod<-lapply(myEM@em.models, function(x)x@model)<%6%5%>  sel.mod<-unique(do.call(c, sel.mod[grep("wmean",names(sel.mod))]))<%6%5%>
 myproj<-BIOMOD_Projection(modeling.output=mymodel,<%6%5%> new.env=vars_proj,
proj.name="current",
selected.models=sel.mod,
do.stack=T)

 myEF<-BIOMOD_EnsembleForecasting(EM.output=myEM,<%6%5%> projection.output=myproj,
total.consensus=T,
binary.meth=c("TSS"))

 myproj<-BIOMOD_Projection(modeling.output=mymodel,<%6%5%> new.env=future,
proj.name="future",
selected.models=sel.mod,
do.stack=T)

 myEF<-BIOMOD_EnsembleForecasting(EM.output=myEM,<%6%5%> projection.output=myproj,
total.consensus=T,
binary.meth=c("TSS"))
}

At the end of this procedure we have 6 ensemble SDMs for P. krameri, 3 for current and 3 for future climatic conditions at a global scale. In the next module we will produce the 2 final ensemble SDMs by averaging these models, and we will do some test to evaluate the goodness of our models.

References

Araújo MB & New M. 2006. Ensemble forecasting of species distributions. TREE 22: 42-47


IPCC (2013) Climatic Change 2013: The Physical Science Basis. Cambridge University Press, Cambridge, United Kingdom.