The output interpretation and test of robustness

Overview

In this module we will build the final model by averaging the 3 ensemble models (one for each replicate) obtained in module 3. A weighted average will be obtained by using the TSS scores associated to each ensemble model that will be extracted and plotted with the ROC scores as well.

Obtaining the final model and Interpretation of the results

In this section, based on six steps, we will get an avereged model using TSS scores associated to each ensemble model obtained in module 3.

Obtaining the final model and Interpretation of the results - Step 1/6

First of all, we need to load the output of the ensemble models.

l<-list.files(getwd(), "ensemble.models.out", recursive=T, full.names=T)
> for(i in 1:length(l))
+ {
+ load(l[i])
+ }
>
> myEM<-lapply(grep("ensemble.models.out", ls(), value=T), get)

Loading required package: biomod2
Loading required package: sp
Loading required package: raster
Loading required package: parallel
Loading required package: reshape
biomod2 3.1-54 loaded.

Now, we will get the TSS and ROC scores with the get_evaluations() function of biomod2.

> library(biomod2)
> TSS<-unlist(lapply(myEM, get_evaluations), recursive=F)

By default, when asked to build weighted ensemble models, biomod2 calculates an ensemble for each evaluation metric used as weight. In this case, we will consider just the ensemble obtained averaging the single models’ predictions with the TSS as weight. The “TSS” object in the strings below will be used later on as the weight vector for the final ensemble.

> ROC<-sapply(TSS[grep("TSS", names(TSS))], "[[", c(1))
> TSS<-sapply(TSS[grep("TSS", names(TSS))], "[[", c(2))

Have a look at the TSS scores

> TSS

Obtaining the final model and Interpretation of the results - Step 2/6

Now we load the six ensembles, 3 for the current and 3 for future climate conditions, and perform the weighted mean with the “TSS” object as weight.

> m<-list.files(getwd(), "ensemble.gri", recursive=T, full.names=T)
> m_current<-lapply(grep("current", m, value=T), stack)
> m_future<-lapply(grep("future", m, value=T), stack)

Here we have 2 esemble models for each replicates. We want only those obtained using TSS as evaluation metrics

> m_current<-lapply(m_current, "[[", c(2))
> m_future<-lapply(m_future, "[[", c(2))

Finally we compute the final ensemble model by averaging the three ensemble models by weighting them according to the respective TSS values

> current_final<-weighted.mean(stack(m_current), w=TSS)
> future_final<-weighted.mean(stack(m_future), w=TSS)

To obtain maps of the final models, we provided in the tutorial a script with a palette function the reproduce the same colour scale used by the MAXENT software. The function requires the minimum and maximum value of the raster intended to be plotted and the interval of values used to shift from a colour to another.

> source("script/palette_MAXENT.R")
> plot(current_final, col=palette_MAXENT(min=0, max=1000, interval=100)[[1]], breaks=palette_MAXENT(0, 1000, 100)[[2]])

Obtaining the final model and Interpretation of the results - Step 3/6

It might be convenient to plot the output in a binary format (presence/absence). The threshold for presence was calculated using the TSS metrics as specified in the BIOMOD_EnsembleForecasting() function (see biomod2 help for details). Being these maps made of zeros and ones, we will not produce a weighted mean, but we will sum them and retain as cells of predicted presence just those where at least 2 of the 3 binary maps predicted the species as present. We first load the binary maps:

> b<-list.files(getwd(), "ensemble_TSSbin.gri", recursive=T, full.names=T)
> b_current<-lapply(grep("current", b, value=T), stack)
> b_future<-lapply(grep("future", b, value=T), stack)

Retain only maps obtained by using TSS as evaluation metrics:

> b_current<-lapply(b_current, "[[", c(2))
> b_future<-lapply(b_future, "[[", c(2))

Finally, we sum the cells that predicted the presence of the species in at least 2 replicates

> b_current<-sum(stack(b_current))
> values(b_current)[values(b_current)%in%c(0:1)]<-0
> values(b_current)[values(b_current)%in%c(2:3)]<-1
> b_future<-sum(stack(b_future))
> values(b_future)[values(b_future)%in%c(0:1)]<-0
> values(b_future)[values(b_future)%in%c(2:3)]<-1

Plot the resulting map of presence/absence under current climate

> plot(b_current, col=c("darkgrey", "darkgreen"))

Obtaining the final model and Interpretation of the results - Step 4/6

and under future climate conditions

> plot(b_current, col=c("darkgrey", "darkgreen"))

Obtaining the final model and Interpretation of the results - Step 5/6

According to these four maps the probability to find P. krameri in its native range will be lower in many areas in the future. By contrast, new potential suitable area will be present outside its native range.
 
Now we can graphically evaluate the effect of climate change on the potential distribution of the species making a histogram of the values of the probability maps for the current and future climates:

> library(scales)
> par(mfrow=c(1,3))
> hist(values(current_final), xlim=c(0,1000), xlab="Probability of presence*1000", col="darkblue", main="Predicted probabilies of presence")
> hist(values(future_final), xlim=c(0,1000), xlab="Probability of presence*1000", col=alpha("red",0.7), main="Future climate", add=T)

In addition, it is possible to compare the number of cells of predicted presence for the current and future climates, in order to evaluate if the species is predicted to colonize new areas.

> barplot(c(table(values(b_current))[2], table(values(b_future))[2]), col="darkgreen", main="Area of predicted presence",
+names=c("Current", "Future"), ylab="number of cells")

Obtaining the final model and Interpretation of the results - Step 6/6

Finally, a boxplot can be used to evaluate if the future predicted distribution would change its prevalent latitude. First of all, the rasterToPoints() function can be used to extract the coordinates of each raster cell for both the current and the future predicted distributions.
Then, it is necessary to extract just the latitudes of the cells with the species predicted as present and finally, to launch the boxplot.

> rtp_current<-rasterToPoints(b_current)
> rtp_current<-rtp_current[rtp_current[,3]==1,2]
> rtp_future<-rasterToPoints(b_future)
> rtp_future<-rtp_future[rtp_future[,3]==1,2]
> boxplot(list(rtp_current, rtp_future), col="grey", main="Latitude of predicted presence",
+names=c("Current", "Future"), ylab="decimal degrees")

The first panel on the left of the figure above clearly shows how the probability values for the future climate (the red bars) result shifted toward higher values with respect of the current climate (the blue bars). Furthermore, the species will reduce its potential suitable habitat in the future (green histograms in the middle). Despite reduced new areas will became potentially hospitable for this species, especially in northern latitudes (boxplot on the right).

Evaluate the robustness of the projections: the MESS index

In this section, based on two steps, we will use Multivariate Environmental Similarity Surfaces (MESS) for projections and evaluate their robustness.

Evaluate the robustness of the projections: the MESS index - Step 1/2

The projections made outside the model training range (the native range) should be treated carefully. An artefact in the extrapolation is possible and it is very important to evaluate where the extrapolation could be particularly critical. Elith et al. (2010) proposed a Multivariate Environmental Similarity Surfaces (MESS) index to take into account this issue. The MESS is an index of similarity reporting the closeness of a point described by a set of environmental attributes to the distribution of these attributes within a population of reference points.
In the strings below, we will load one of the biomod output, called “formated.input.data”, from where it is possible to get the coordinates of the points used for the training of the model (we will perform the MESS analysis just for the first replicate). Then we will re-load the rasters of the environmental variables selected in the previous models and finally apply the mess() function provided by the dismo R package.

> ## load single models files
> s<-list.files(getwd(), "formated.input.data", recursive=T, full.names=T, all.files=T)<%6%5%> > load(s[1])
> s<-grep("data", ls(), value=T)<%6%5%> > mydata<-get(s)<%6%5%>
> # extract training dataset
> dataset<-mydata@coord<%6%5%>
> # perform MESS analysis (it may takes several minutes)
> library(dismo)
> library(gtools)

The mess function requires the stack of the environmental variables (vars_proj) and the values of the same variables extracted at the occurrence points:

> r <- getData("worldclim", var = "bio", res = 2.5)<%6%5%> > vars_training<-stack(mixedsort(list.files("predictors\\", "gri", full.names=T)))<%6%5%> > r<-r[[which(names(r)%in%names(vars_training))]]<%6%5%> > vars_proj<-stack(aggregate(r, 5, progress="text"))<%6%5%>

Evaluate the robustness of the projections: the MESS index - Step 2/2

Now we can calculate the MESS index and plot the map

> MESS<-mess(X=vars_proj, V=extract(vars_training, dataset), full=F)
> color.high<-colorRampPalette(c("lightgrey","blue"))
> color.low<-colorRampPalette(c("red","lightgrey"))
> plot(MESS, col=c(color.low(abs(minValue(MESS))),color.high(maxValue(MESS))))

In the figure plot, the red values indicate areas where the extrapolation effect results more serious. According to this map most of the areas outside the native range where the probability to find a suitable area for the species is high (according to the ensemble models) is not affected by serious problems. However there are still many regions where the projection must be interpreted carefully.

Final remarks

The procedure showed in this tutorial is not the result of a planned research but is just an exemplary workflow for SDM. We showed how to download occurrences data, how to clean them, how to select environmental variables and how to fit models. Finally, we showed how to obtain good maps and test useful for a publication. You must be aware that the dataset that we used has many limitations and the results obtained must not be over interpreted. It is possible that you obtain a different result when you will run the scripts provided in the tutorial. The scope of the tutorial was not to show the results of a research but just to provide a tool to start to understand how SDMs works in practice. After the reading of this tutorial It is important that you do some experimentation: try to modify settings in the procedure and see what’s happen. Take this tutorial for what it is. Do NOT use it as a black box trying to force your data in. The results will be probably incorrect. It is very important that you read recent papers and other tutorials to understand how SDMs works not only in practice but also the theory that is behind species distribution modelling.