Species distribution model - African elephant
The goal is to build a SDM for the African elephant. A pre-processed dataset from Angelov, 2020 can be found in the EcoData package which is only available on github:
set.seed(1337)
if(!require(EcoData)) devtools::install_github(repo = "TheoreticalEcology/EcoData",
dependencies = FALSE, build_vignettes = FALSE)
#> Loading required package: EcoData
library(EcoData)
df = EcoData::elephant$occurenceData
head(df)
#> Presence bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12
#> 3364 0 -0.4981747 -0.2738045 0.5368968 -0.5409999 -0.36843571 0.2947850 -0.5260099 -1.2253960 0.2494100 -0.64527314 -0.06267842 0.6285371
#> 6268 0 0.6085908 -0.5568352 1.0340686 -1.2492050 -0.11835651 0.8221087 -0.8938475 0.4233787 0.7746249 0.09168503 0.94419518 1.1121516
#> 10285 0 -0.7973005 1.4648130 -1.0540532 2.0759423 0.07614953 -1.5860029 1.6284678 0.2768209 -1.5153122 -0.03648161 -1.44165748 -1.2351482
#> 2247 0 0.6385034 1.3435141 -0.1591439 -0.5107148 1.10425291 -0.1622288 0.8577603 0.4600181 0.5855475 0.54026827 0.68153250 0.5951165
#> 9821 0 0.6684160 -0.6781341 0.6363311 -0.9906170 0.15950927 0.9099960 -0.8062671 0.3867393 0.8586593 0.31597665 0.94419518 1.1003561
#> 1351 0 0.9675418 -0.6781341 -0.3580126 -0.3748202 0.77081398 0.8748411 -0.3858812 0.3134604 1.0477367 0.98885151 0.94419518 0.7287986
#> bio13 bio14 bio15 bio16 bio17 bio18 bio19
#> 3364 0.6807958 -0.29703736 -0.008455252 0.7124535 -0.2949994 -1.06812752 1.96201807
#> 6268 0.5918442 0.01619202 -0.884507980 0.5607328 0.3506918 1.22589281 -0.36205814
#> 10285 -1.3396742 -0.50585695 0.201797403 -1.3499999 -0.5616980 -0.42763181 -0.62895735
#> 2247 0.8714061 -0.55806185 0.236839512 1.1012378 -0.5616980 -0.20541902 -0.58378979
#> 9821 0.5537222 0.59044589 -1.024676416 0.6413344 0.7437213 0.06254347 -0.05409751
#> 1351 1.1255533 -0.50585695 0.236839512 1.2956300 -0.4494038 -0.90473576 2.47939193Presence is our response variable and we have the 19 bioclim variables as features/predictors.
Let’s split part of the data away so that we can use it at the end to evaluate our model:
Adjusting optimization parameters - Convergence
We will first try a simple DNN with default values and the binomial likelihood. We use 30% of the data as validation holdout to check for overfitting:
library(cito)
model = dnn(Presence~., data = df,
batchsize = 100L,
validation = 0.3, loss = "binomial",
verbose = FALSE)
We see that the training and test losses were still decreasing which means we didn’t train the model long enough. We could now either increase the number of epochs or increase the learning rate so that the model trains faster:
model = dnn(Presence~., data = df,
batchsize = 100L,
lr = 0.05,
validation = 0.3, loss = "binomial",
verbose = FALSE)
Much better! But still now enough epochs. Also, let’s see if we can further decrease the loss by using a wider and deeper neural network:
model = dnn(Presence~., data = df,
batchsize = 100L,
hidden = c(100L, 100L, 100L),
lr = 0.05,
validation = 0.3, loss = "binomial",
verbose = FALSE)
At the end of the training, the losses start to get jumpy, which can be a sign of potential overfitting. We can control that by adding a weak regularization (as we want a L2 regularization, we set alpha to 1.0):
model = dnn(Presence~., data = df,
batchsize = 100L,
epochs = 150L,
hidden = c(100L, 100L, 100L),
lr = 0.05,
lambda = 0.001,
alpha = 1.0,
validation = 0.3, loss = "binomial",
verbose = FALSE)
We will turn on now advanced features that help with the convergence and to reduce overfitting:
learning rate scheduler - reduces learning rate during training
early stopping - stop training when validation loss starts to increase
model = dnn(Presence~., data = df,
batchsize = 100L,
epochs = 150L,
hidden = c(100L, 100L, 100L),
lr = 0.05,
lambda = 0.001,
alpha = 1.0,
validation = 0.3, loss = "binomial",
verbose = FALSE,
lr_scheduler = config_lr_scheduler("reduce_on_plateau", patience = 7), # reduce learning rate each 7 epochs if the validation loss didn't decrease,
early_stopping = 14 # stop training when validation loss didn't decrease for 10 epochs
)
Great! We found now a model architecture and training procedure that fits and trains well. Let’s proceed to our final model
Train final model with bootstrapping to obtain uncertainties
We haven’t directly started with bootstrapping because it complicates the adjustment of the training procedure.
Uncertainties can be obtained by using bootstrapping. Be aware that this can be computational expensive:
model_boot = dnn(Presence~., data = df,
batchsize = 100L,
epochs = 150L,
hidden = c(100L, 100L, 100L),
lr = 0.05,
lambda = 0.001,
alpha = 1.0,
validation = 0.3, loss = "binomial",
verbose = FALSE,
lr_scheduler = config_lr_scheduler("reduce_on_plateau", patience = 7), # reduce learning rate each 7 epochs if the validation loss didn't decrease,
early_stopping = 14, # stop training when validation loss didn't decrease for 10 epochs
bootstrap = 20L,
bootstrap_parallel = 5L
)Predictions
We can use the model now for predictions:
The predictions are 2/3 dimensional because of the bootstrapping. Calculate the AUC interval:
hist(sapply(1:20, function(i) Metrics::auc(test$Presence, predictions[i,,])),
xlim = c(0, 1), main = "AUC of ensemble model", xlab = "AUC")
We can now predict the habitat suitability of the elephant (Note that spatial dependencies are required):
library(raster)
#> Loading required package: sp
library(sp)
library(rsample)
library(latticeExtra)
#> Loading required package: lattice
#>
#> Attaching package: 'lattice'
#> The following object is masked from 'package:EcoData':
#>
#> melanoma
library(sp)
library(ggplot2)
#> Need help getting started? Try the R Graphics Cookbook: https://r-graphics.org
#>
#> Attaching package: 'ggplot2'
#> The following object is masked from 'package:latticeExtra':
#>
#> layer
library(maptools)
#> Please note that 'maptools' will be retired during October 2023,
#> plan transition at your earliest convenience (see
#> https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
#> for guidance);some functionality will be moved to 'sp'.
#> Checking rgeos availability: FALSE
customPredictFun = function(model, data) {
return(apply(predict(model, data, reduce = "none"), 2:3, mean)[,1])
}
normalized_raster = EcoData::elephant$predictionData
predictions =
raster::predict(normalized_raster,
model_boot,
fun = customPredictFun)
habitat_plot =
spplot(predictions, colorkey = list(space = "left") )
habitat_plot
Moreover, we can visualize the uncertainty of our model, instead of calculating the average occurrence probability, we calculate for each prediction the standard deviation and visualize it:
customPredictFun_sd = function(model, data) {
return(apply(predict(model, data, reduce="none"), 2:3, sd)[,1])
}
predictions =
raster::predict(normalized_raster,
model_boot,
fun = customPredictFun_sd)
uncertain_plot =
spplot(predictions, colorkey = list(space = "left") )
uncertain_plot
Inference
Neural networks are often called black-box models but the tools of explainable AI (xAI) allows us to understand them - and actually infer properties similar to what a linear regression model can provide (the calculation can take some time…):
results = summary(model_boot)
results
#> Summary of Deep Neural Network Model
#>
#> ── Feature Importance
#> Importance Std.Err Z value Pr(>|z|)
#> bio1 → 0.5671 0.2955 1.92 0.05493 .
#> bio2 → 0.3368 0.2757 1.22 0.22188
#> bio3 → 0.1971 0.1093 1.80 0.07128 .
#> bio4 → 0.2896 0.1992 1.45 0.14594
#> bio5 → 0.6996 0.3356 2.09 0.03707 *
#> bio6 → 0.5888 0.4589 1.28 0.19940
#> bio7 → 0.1762 0.1117 1.58 0.11486
#> bio8 → 0.3758 0.1940 1.94 0.05270 .
#> bio9 → 2.5090 1.4337 1.75 0.08012 .
#> bio10 → 0.3349 0.2130 1.57 0.11585
#> bio11 → 0.2795 0.2177 1.28 0.19921
#> bio12 → 1.9087 0.6461 2.95 0.00314 **
#> bio13 → 0.3197 0.1760 1.82 0.06924 .
#> bio14 → 0.4169 0.2696 1.55 0.12200
#> bio15 → 0.8942 0.2640 3.39 0.00071 ***
#> bio16 → 0.3672 0.3013 1.22 0.22299
#> bio17 → 0.1733 0.1445 1.20 0.23038
#> bio18 → 0.2410 0.1145 2.10 0.03536 *
#> bio19 → 0.0645 0.0355 1.82 0.06914 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ── Average Conditional Effects
#> ACE Std.Err Z value Pr(>|z|)
#> bio1 → 0.09733 0.03226 3.02 0.00256 **
#> bio2 → -0.05851 0.03811 -1.54 0.12465
#> bio3 → -0.00413 0.03759 -0.11 0.91261
#> bio4 → 0.03089 0.03336 0.93 0.35445
#> bio5 → 0.11839 0.03709 3.19 0.00141 **
#> bio6 → 0.09144 0.04549 2.01 0.04443 *
#> bio7 → 0.00310 0.03884 0.08 0.93644
#> bio8 → -0.04861 0.03167 -1.54 0.12475
#> bio9 → -0.18590 0.05309 -3.50 0.00046 ***
#> bio10 → -0.03716 0.04174 -0.89 0.37333
#> bio11 → 0.02729 0.04083 0.67 0.50397
#> bio12 → -0.19637 0.03947 -4.97 6.5e-07 ***
#> bio13 → 0.07646 0.03118 2.45 0.01420 *
#> bio14 → 0.09271 0.04582 2.02 0.04302 *
#> bio15 → -0.09483 0.02640 -3.59 0.00033 ***
#> bio16 → -0.05575 0.04044 -1.38 0.16803
#> bio17 → 0.03714 0.04335 0.86 0.39166
#> bio18 → 0.02949 0.01743 1.69 0.09069 .
#> bio19 → 0.00413 0.02446 0.17 0.86604
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ── Standard Deviation of Conditional Effects
#> ACE Std.Err Z value Pr(>|z|)
#> bio1 → 0.1238 0.0321 3.86 0.00012 ***
#> bio2 → 0.0827 0.0331 2.50 0.01243 *
#> bio3 → 0.0548 0.0230 2.38 0.01734 *
#> bio4 → 0.0717 0.0258 2.77 0.00553 **
#> bio5 → 0.1366 0.0359 3.81 0.00014 ***
#> bio6 → 0.1116 0.0438 2.55 0.01081 *
#> bio7 → 0.0565 0.0206 2.74 0.00618 **
#> bio8 → 0.0946 0.0224 4.22 2.5e-05 ***
#> bio9 → 0.2123 0.0601 3.53 0.00041 ***
#> bio10 → 0.0792 0.0235 3.37 0.00076 ***
#> bio11 → 0.0859 0.0205 4.18 2.9e-05 ***
#> bio12 → 0.2167 0.0356 6.08 1.2e-09 ***
#> bio13 → 0.0877 0.0267 3.28 0.00102 **
#> bio14 → 0.1134 0.0380 2.98 0.00285 **
#> bio15 → 0.1093 0.0189 5.79 7.0e-09 ***
#> bio16 → 0.0808 0.0343 2.35 0.01856 *
#> bio17 → 0.0671 0.0291 2.30 0.02123 *
#> bio18 → 0.0949 0.0265 3.59 0.00033 ***
#> bio19 → 0.0380 0.0152 2.50 0.01257 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Bioclim9, 12, 14, and 16 have large significant average conditional effects ($\approx$ linear effects). We can visualize them using accumulated local effect plots:
