This code is based on the vignette of the dismo package (Hijmans & Elith 2011)

Load the required packages:

library(rgeos)
library(rgdal)
library(sp)
library(raster)
library(maptools)
library(dismo)

If we did not load our last workspace (from the ranges’ exercise), we need to do it

load("~/Dropbox/Madagascar_BITC/exercise_ranges.RData")

Ecological Niche Model (Bioclim)

Get environmental variables fro Africa(load them, we’ll only use one)

afr.bio1 <- raster ("Madagascar_BITC_data/bio_1.asc")
afr.bio12 <- raster ("Madagascar_BITC_data/bio_12.asc")
afr.env.vars <- stack(afr.bio1,afr.bio12)

Set the extent of our region (Madagascar) and “cut” the environmental variables to that extent

mdg.extent <- extent (43, 51, -26, -12)
mdg.env.vars <- crop(afr.env.vars, mdg.extent)

Create the Bioclim model (Ecological Niche Model: ENM)

# Convert our sp points into a `matrix` object
fossa_points_mat <- as.matrix(fossa_points_nonas)

MODEL FITTING

#Run the ENM model (Bioclim) (i.e. from the occurrrence data - geographic space -  to the environmental space and applying a function to relate these two spaces: DUALITY!!)
fossa_bioclim <- bioclim(mdg.env.vars, fossa_points_mat)

plot(fossa_bioclim, a=1, b=2, p=0.85)

We can used this model to explore how the species responses to these particular variables (based on the particular algorithm that we used and the relationship between the geographic and environmental space), using a response function that creates response plots for each variable, with the other variables at their median value.

response(fossa_bioclim)

MODEL PREDICTION

#Create a prediction of the suitable area based on the ENM model (i.e. from the environmental space to the geographical space: DUALITY!!!)
fossa.map <- predict (mdg.env.vars, fossa_bioclim)

#Plot the resulting prediction based on the ENM model
plot (fossa.map)
#Add the actual points occurrence of the species (records)
points(fossa_points_mat[,1],fossa_points_mat[,2],pch='x',)

How does it look?

Are there any regions where the species is predicted to be but it’s absent?

MODEL EVALUATION

Is it a good model? We need to evaluate it!

Data partition

#Let’s make a training and a testing set.
group <- kfold(fossa_points_mat, 3)
pres_train <- fossa_points_mat[group != 1, ]
pres_test <- fossa_points_mat[group == 1, ]

#Background data for training and a testing set. 
backg <- randomPoints(mdg.env.vars, n=100) 
colnames(backg) = c('lon', 'lat')
group <- kfold(backg, 3)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]

#The first layer in the Raster- Stack is used as a ’mask’. That ensures that random points only occur within the spatial extent of the rasters, and within cells that are not NA, and that there is only a single absence point per cell. Here we further restrict the background points to be within 12.5% of our specified extent ’ext’.
r = raster(mdg.env.vars, 1)
plot(!is.na(r), col=c('white', 'light grey'), legend=FALSE) 
points(backg_train, pch='-', cex=0.5, col='yellow')
points(backg_test, pch='-', cex=0.5, col='black')
points(pres_train, pch= '+', col='green')
points(pres_test, pch='+', col='blue')

Model fitting and evaluation using training data as well as all data

#Fit the Bioclim model only for the training points
fossa_bc2 <- bioclim(mdg.env.vars, pres_train)

#Evaluate the model, by providing presence and back- ground (absence) points, the model, and a RasterStack:
e <- evaluate(pres_test, backg_test, fossa_bc2, mdg.env.vars)
e

#What about the model with all data?
e2 <- evaluate(pres_test, backg_test, fossa_bioclim, mdg.env.vars)
e2

# max TPR+TNR: maximum of the sum of the sensitivity (true positive rate) and specificity (true negative rate). This is sometimes used as a threshold for setting cells to presence or absence.

Predicted values and Binary map

Use the predictor variables to make a prediction to a RasterLayer and then use a threshold to create a Binary Map

pb <- predict(mdg.env.vars, fossa_bc2, progress='') 
pb

#Plot the results from the prediction and the thresholded map
par(mfrow=c(1,2))
plot(pb, main='Bioclim raw values')
plot(wrld_simpl, add=TRUE, border='dark grey')
threshold <- e@t[which.max(e@TPR + e@TNR)]
plot(pb > threshold, main='presence/absence')
plot(wrld_simpl, add=TRUE, border='dark grey')
points(pres_train, pch='+')

#Before, we use the "recommended" threshold. We can use another one!
par(mfrow=c(1,2))
plot(pb, main='Bioclim raw values')
plot(wrld_simpl, add=TRUE, border='dark grey')
#Our own threshold
threshold <- 0.5
plot(pb > 0.5, main='presence/absence')
plot(wrld_simpl, add=TRUE, border='dark grey') 
points(pres_train, pch='+')

#Predict to the raster using the Bioclim model with all the data
pbOrig <- predict(mdg.env.vars, fossa_bioclim, progress='') 
pbOrig

#Plot the results from the prediction and the thresholded map
par(mfrow=c(1,2))
plot(pbOrig, main='Bioclim raw values')
plot(wrld_simpl, add=TRUE, border='dark grey')
threshold <- e2@t[which.max(e2@TPR + e2@TNR)]
plot(pbOrig > threshold, main='presence/absence')
plot(wrld_simpl, add=TRUE, border='dark grey') 
points(pres_train, pch='+')

MODEL TRANSFERABILITY

The exercise above (prediction) was only for Madagascar. But, can we predict to another region?

#Use the fitted model to predict onto another region (continental Africa)
fossa.map.afr <- predict (afr.env.vars, fossa_bioclim)

#Plot the model's prediction
plot (fossa.map.afr)

How does it look? Where are the favorable conditions for our species in continental Africa?

What about in the future?!

#Let's get the same environmental variables for the future (50 years)
afr.bio1.fut <- raster ("Madagascar_BITC_data/bio_1_fut.asc")
afr.bio12.fut <- raster ("Madagascar_BITC_data/bio_12_fut.asc")
afr.env.vars.fut <- stack(afr.bio1.fut,afr.bio12.fut)
#Now "cut"" them for Madagascar
mdg.env.vars.fut <- crop(afr.env.vars.fut, mdg.extent)

#Now, let's use the fitted model to predict onto another time (Madagascar in 50 years!)
### NOTE: the variable names need to be IDENTICAL!
#check the variable names
names(mdg.env.vars)
names(mdg.env.vars.fut)
#change if necessary
names(mdg.env.vars.fut) <- c("bio_1","bio_12")


#Predict to the future
fossa.map.fut <- predict (mdg.env.vars.fut, fossa_bioclim)

#Plot the model's prediction
plot (fossa.map.fut)

How does it look?

#Plot both maps to see the potential differences
par(mfrow=c(1,2))
plot(fossa.map)
plot(fossa.map.fut)

What do you think? Is the (environmetal suitability of) fossa going to change in 50 years? How much?