Tsetse-fly-Ensemble-Modelling-Kenya.R

# ------------------------------------------------------------
# Script: Tsetse Habitat Suitability Modelling - Kenya
# Purpose: Multicollinearity check and modelling purposes
# Context: Supports tsetse habitat suitability modelling in Kenya
# Author: Raphael Mongare
# R Version: 4.2.3
# ------------------------------------------------------------

# Install the relevant packages

install.packages('corrplot')
install.packages('ggpubr')
install.packages('PresenceAbsence')
install.packages('RColorBrewer')
install.packages('quickPlot')
install.packages('rgdal')
install.packages('vip')
install.packages('pdp')
install.packages('boot')
install.packages('raster')
install.packages('foreign')
install.packages('sdm')
install.packages('gbm')
install.packages('goeveg')
install.packages('mapview')


# Load the required packages
library(corrplot)
library(ggpubr)
library("PresenceAbsence")
library(RColorBrewer)
library(quickPlot)
library(vip)
library(pdp)
library(boot)
library(raster)
library(foreign)
library(sdm)
library(gbm)
library(rgdal)
library(goeveg)
library(mapview)
library(vioplot)
sdm::installAll()

#-------------------------------------------------------------------------------
## set the working directory####
setwd('C:Raw_Data_Annual\\Processed_V4')

## load the database which has the extracted predictor raster values####
d <- read.csv('Tsetse-predictors-Cleaned')


#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
## Use the following function to perfom VIF to check for multi-collinearity
## Copy and paste it in a new script and add it to source

## VIF Calculation(Add to source)
## T.A. GROEN (GROEN@ITC.NL)
## JUNE 2013
VIFcalc<-function(d) # d is the dataframe of your predictor variables
{
  result<-data.frame(var=c(names(d)),
                     VIF=numeric(length(d[1,])))
  for(i in (1:length(d[1,]))) 
  {
    result$VIF[i] <-1/(1-summary(lm(d[,i] ~ .,data=d[,names(d)!=names(d)[i]]))$r.squared)
  }
  return(result)
}

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------


## Value of inflation (VIF) ####
# drop the variables based on importance to tsetse flies##
## Do each variable separetly per year###

## Compute VIF for predictors yearly (2017 - 2020) and combine the retained non-collinear
## predictors with the static variables in the end. Retain variables with VIF = <10.
####   2017    ####
# vif_2017 <- VIFcalc(data.frame(d$maxSkinTEMP17,d$minSkinTEMP17,d$medSkinTEMP17,
#                                d$maxSoilTEMP17,d$minSoilTEMP17,d$medSoilTEMP17,
#                                d$maxAIRTEMP17,d$minAIRTEMP17,d$medAIRTEMP17,
#                                d$maxLSTNight17,d$minLSTNight17,d$medLSTNight17,
#                                d$maxLSTDay17,d$minLSTDay17,d$medLSTDay17,
#                                d$sumTOTPREC17,d$sumSurfRunoff17,
#                                d$maxNDVI17,d$minNDVI17,d$medNDVI17,
#                                d$maxSbSMOIST17,d$minSbSMOIST17,d$medSbSMOIST17,
#                                d$maxSSMOIST17,d$minSSMOIST17,d$medSSMOIST17))

vif_2017 <- VIFcalc(data.frame(d$maxSkinTEMP17,
                               d$minSoilTEMP17,
                               d$minLSTNight17,
                               d$maxLSTDay17,d$minLSTDay17,
                               d$sumTOTPREC17,d$sumSurfRunoff17,
                               d$maxNDVI17,d$minNDVI17,d$medNDVI17,
                               d$maxSbSMOIST17,
                               d$maxSSMOIST17,d$minSSMOIST17))

vif_2017

####   2018    ####

# vif_2018 <- VIFcalc(data.frame(d$maxSkinTEMP18,d$minSkinTEMP18,d$medSkinTEMP18,
#                                d$maxSoilTEMP18,d$minSoilTEMP18,d$medSoilTEMP18,
#                                d$maxAIRTEMP18,d$minAIRTEMP18,d$medAIRTEMP18,
#                                d$maxLSTDay18,d$minLSTDay18,d$medLSTDay18,
#                                d$maxLSTNight18,d$minLSTNight18,d$medLSTNight18,
#                                d$sumTOTPREC18,d$sumSurfRunoff18,
#                                d$maxNDVI18,d$minNDVI18,d$medNDVI18,
#                                d$maxSbSMOIST18,d$minSbSMOIST18,d$medSbSMOIST18,
#                                d$maxSSMOIST18,d$minSSMOIST18,d$medSSMOIST18))

vif_2018 <- VIFcalc(data.frame(
  d$maxSoilTEMP18,
  d$minAIRTEMP18,
  d$maxLSTDay18,d$minLSTDay18,
  d$maxLSTNight18,d$minLSTNight18,
  d$sumTOTPREC18,d$sumSurfRunoff18,
  d$maxNDVI18,d$minNDVI18,d$medNDVI18,
  d$maxSbSMOIST18,
  d$maxSSMOIST18,d$minSSMOIST18,d$medSSMOIST18))
vif_2018

####   2019    ####

# vif_2019 <- VIFcalc(data.frame(d$maxSkinTEMP19,d$minSkinTEMP19,d$medSkinTEMP19,
#                                d$maxSoilTEMP19,d$minSoilTEMP19,d$medSoilTEMP19,
#                                d$maxAIRTEMP19,d$minAIRTEMP19,d$medAIRTEMP19,
#                                d$maxLSTDay19,d$minLSTDay19,d$medLSTDay19,
#                                d$maxLSTNight19,d$minLSTNight19,d$medLSTNight19,
#                                d$sumTOTPREC19,d$sumSurfRunoff19,
#                                d$maxNDVI19,d$minNDVI19,d$medNDVI19,
#                                d$maxSbSMOIST19,d$minSbSMOIST19,d$medSbSMOIST19,
#                                d$maxSSMOIST19,d$minSSMOIST19,d$medSSMOIST19))

vif_2019 <- VIFcalc(data.frame(d$maxSkinTEMP19,d$minSkinTEMP19,
                               d$maxLSTDay19,d$minLSTDay19,
                               d$maxLSTNight19,d$minLSTNight19,
                               d$sumSurfRunoff19,
                               d$maxNDVI19,d$minNDVI19,d$medNDVI19,
                               d$maxSbSMOIST19,
                               d$maxSSMOIST19,d$minSSMOIST19,d$medSSMOIST19))
vif_2019

####   2020    ####

# vif_2020 <- VIFcalc(data.frame(d$maxSkinTEMP20,d$minSkinTEMP20,d$medSkinTEMP20,
#                                d$maxSoilTEMP20,d$minSoilTEMP20,d$medSoilTEMP20,
#                                d$maxAIRTEMP20,d$minAIRTEMP20,d$medAIRTEMP20,
#                                d$maxLSTDay20,d$minLSTDay20,d$medLSTDay20,
#                                d$maxLSTNight20,d$minLSTNight20,d$medLSTNight20,
#                                d$sumTOTPREC20,d$sumSurfRunoff20,
#                                d$maxNDVI20,d$minNDVI20,d$medNDVI20,
#                                d$maxSbSMOIST20,d$minSbSMOIST20,d$medSbSMOIST20,
#                                d$maxSSMOIST20,d$minSSMOIST20,d$medSSMOIST20))

vif_2020 <- VIFcalc(data.frame(d$maxSkinTEMP20,
                               d$minAIRTEMP20,
                               d$maxLSTDay20,d$minLSTDay20,
                               d$maxLSTNight20,d$minLSTNight20,
                               d$sumTOTPREC20,d$sumSurfRunoff20,
                               d$maxNDVI20,d$minNDVI20,d$medNDVI20,
                               d$maxSbSMOIST20,d$minSbSMOIST20,
                               d$maxSSMOIST20))
vif_2020


### Address multicollinearity for the retained variables across the years + static variables
### Retain variables with VIF = <10
vif_17_20 <- VIFcalc(data.frame(d$maxSkinTEMP17,
                                d$minSoilTEMP17,
                                d$minLSTNight17,
                                d$minLSTDay17,
                                d$sumSurfRunoff17,
                                d$maxNDVI17,d$minNDVI17,
                                d$maxSbSMOIST17,
                                d$maxLSTDay18,d$minLSTDay18,
                                d$minLSTNight18,
                                d$sumTOTPREC18,d$sumSurfRunoff18,
                                d$maxNDVI18,d$minNDVI18,
                                d$maxSSMOIST18,d$minSSMOIST18,
                                d$maxLSTDay19,d$minLSTDay19,
                                d$minLSTNight19,
                                d$maxNDVI19,d$minNDVI19,
                                d$minLSTDay20,
                                d$minLSTNight20,
                                d$maxNDVI20,d$minNDVI20,
                                d$minSbSMOIST20,
                                d$sand_combat,
                                d$NewTWI,d$Kenya_Human_Pop_2020_1km,
                                d$Cattle_Density_2020,d$EucDist_ProtectedAreas,
                                d$Goat))

vif_17_20

#-------------------------------------------------------------------------------
## Backward stepwise regression to retain statistically significant predictors
glm <- glm(Occ ~ maxSkinTEMP17+minSkinTEMP17+minLSTNight17
           +minLSTDay17+sumSurfRunoff17+maxNDVI17+minNDVI17+maxSbSMOIST17
           +minSSMOIST17+medSSMOIST17+maxLSTDay18+minLSTDay18
           +minLSTNight18+sumSurfRunoff18+maxNDVI18+minNDVI18+medNDVI18
           +minSbSMOIST18+maxSSMOIST18+minSSMOIST18+maxLSTDay19+minLSTDay19
           +minLSTNight19+maxNDVI19+minNDVI19+medSbSMOIST19+maxSSMOIST19
           +minLSTDay20+minLSTNight20+maxNDVI20+minNDVI20+minSbSMOIST20  
           +Elevation+sand_combat+NewTWI+Kenya_Human_Pop_2020_1km
           +Cattle_Density_2020+Sheep+Goat,data=d, family= binomial(link = "logit"))

glm.step <- step(glm) # Retain statistically significant variables

glm.vip <- vi(glm.step) #Calculate the importance score of variables
glm.vip


# View(glm.step)
# summary(glm.step)
# summary(glm.vip)

#-------------------------------------------------------------------------------
## Use the variables retained after applying the step function
## Create a sdm model data
sdm.kenya <-sdmData(formula = Occ ~ maxNDVI17 + minNDVI17 + minSSMOIST17 + 
                      medSSMOIST17 + maxLSTDay18 + minLSTDay18 + minLSTNight18 + 
                      maxNDVI18 + minNDVI18 + medNDVI18 + minSSMOIST18 + 
                      minLSTDay19 + minLSTNight19 + maxNDVI19 + medSbSMOIST19 + 
                      minLSTDay20 + maxNDVI20 + sand_combat + NewTWI + 
                      Kenya_Human_Pop_2020_1km + Cattle_Density_2020 + Sheep + 
                      Goat+Elevation,train=d)


# Create a model objects for (glm, svm, maxent and rf) using the sdm function
sdmEnsmb.genus <-sdm(Occ~.,data=sdm.kenya,methods=c('glm', 'svm', 'maxent', 'rf'), replication='cv',test.percent=30, cv.folds=10)
accglm.kenya<- getEvaluation(sdmEnsmb.genus,wtest="test", opt="max(se+sp)",stat=c("sensitivity","specificity", "AUC","TSS","prevalence","Kappa"))
colMeans(accglm.kenya[sapply(accglm.kenya, is.numeric)])
accglm.kenya

# Save the model object as RDS file for easy retrieval: it takes a lot of time 
# to create the model object for the 4 algorithms

#### saveRDS(sdmEnsmb.genus, 'sdmEnsmb.genus')

# Get the mean variable importance for all the fitted models
vi.kenyaglm <- getVarImp(sdmEnsmb.genus, wtest="test")
vi.kenyaglm  
plot(vi.kenyaglm, geom = 'boxplot')

# Get the ROC curves

glm.roc <- roc(sdmEnsmb.genus, method = 'glm',
               wtest= c('train', 'test'), smooth=T,legend=T)

svm.roc <- roc(sdmEnsmb.genus, method = 'svm',
               wtest= c('train', 'test'), smooth=T,legend=T)

maxent.roc <- roc(sdmEnsmb.genus, method = 'maxent',
               wtest= c('train', 'test'), smooth=T,legend=T)

rf.roc <- roc(sdmEnsmb.genus, method = 'rf',
               wtest= c('train', 'test'), smooth=T,legend=T)

roc(sdmEnsmb.genus,method = c('glm', 'svm', 'maxent', 'rf'),smooth = T, legend=F)

#-------------------------------------------------------------------------------
# Generate the response curves for each variable
curvar <-c('maxNDVI17','minNDVI17','minSSMOIST17', 
               'medSSMOIST17','maxLSTDay18','minLSTDay18','minLSTNight18', 
               'maxNDVI18','minNDVI18','medNDVI18','minSSMOIST18', 
               'minLSTDay19','minLSTNight19','maxNDVI19','medSbSMOIST19', 
               'minLSTDay20','maxNDVI20','sand_combat','NewTWI', 
               'Kenya_Human_Pop_2020_1km','Cattle_Density_2020','Sheep',
               'Goat','Elevation')

curves1 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('Sheep'))
plot(curves1, cex.lab = 8, plot.background = element_erect(fill = 'transparent', color = NA))
ggsave(curves1, filename = 'Sheep.png', bg='transparent') 

curves2 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('Kenya_Human_Pop_2020_1km'))
ggsave(curves2, filename = 'Human population.png', bg='transparent') 

curves3 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('medSSMOIST17'))
ggsave(curves3, filename = 'medSSMOIST17.png', bg='transparent') 


curves4 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('minSSMOIST17'))
ggsave(curves4, filename = 'minSSMOIST17.png', bg='transparent') 

curves5 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('Cattle_Density_2020'))
ggsave(curves5, filename = 'Cattle.png', bg='transparent') 

curves6 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('maxNDVI17'))
ggsave(curves6, filename = 'maxNDVI17.png', bg='transparent') 

curves7 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('Elevation'))
ggsave(curves7, filename = 'Elevation.png', bg='transparent')

curves8 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('maxLSTDay18'))
ggsave(curves8, filename = 'maxLSTDay18.png', bg='transparent')

curves9 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('medSbSMOIST19'))
ggsave(curves9, filename = 'medSbSMOIST19.png', bg='transparent')

curves10 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('medNDVI18'))
ggsave(curves10, filename = 'medNDVI18.png', bg='transparent')

curves11 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('minLSTNight18'))
ggsave(curves11, filename = 'minLSTNight18.png', bg='transparent')

curves12 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('Goat'))
ggsave(curves12, filename = 'Goat.png', bg='transparent')

curves13 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('minNDVI18'))
ggsave(curves13, filename = 'minNDVI18.png', bg='transparent')

curves14 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('minLSTDay18'))
ggsave(curves14, filename = 'minLSTDay18.png', bg='transparent')

curves15 <- rcurve(sdmEnsmb.genus, smooth=T, n = c('minLSTDay20'))
ggsave(curves15, filename = 'minLSTDay20.png', bg='transparent')

#-------------------------------------------------------------------------------
## Stack the raster variables for prediction ###
Var.stack <-stack('maxNDVI17.tif','minNDVI17.tif','minSSMOIST17.tif', 
               'medSSMOIST17.tif','maxLSTDay18.tif','minLSTDay18.tif','minLSTNight18.tif', 
               'maxNDVI18.tif','minNDVI18.tif','medNDVI18.tif','minSSMOIST18.tif', 
               'minLSTDay19.tif','minLSTNight19.tif','maxNDVI19.tif','medSbSMOIST19.tif', 
               'minLSTDay20.tif','maxNDVI20.tif','sand_combat.tif','NewTWI.tif', 
               'Kenya_Human_Pop_2020_1km.tif','Cattle_Density_2020.tif','Sheep.tif',
               'Goat.tif','Elevation.tif')
plot(Var.stack)

#ensmbGenus.pred <- readRDS('ensmbGenus.pred')

start_time <- Sys.time()
ensmbGenus.pred <- predict(sdmEnsmb.genus,newdata=Var.stack,mean=TRUE)#,filename='glm_combat_annual_UV4test1.tif'
end_time <- Sys.time()
end_time - start_time
plot(ensmbGenus.pred)

# Save the trained model file for quick retrieval as it takes time to make predictions
#saveRDS(ensmbGenus.pred,'ensmbGenus.pred')

# Create an ensemble model for the 4 models: GLM, SVM, MaxEnt and RF

start_time <- Sys.time()
ensMod.genus <- ensemble(sdmEnsmb.genus, Varglm, mean=TRUE, setting = list(method = 'weighted', stat = 'AUC'))
end_time <- Sys.time()
end_time - start_time
plot(ensMod.genus)

ensmbEval <- getEvaluation(sdmEnsmb.genus, wtest="test", opt="max(se+sp)",stat=c("sensitivity","specificity", "AUC","TSS","prevalence","Kappa"))
colMeans(ensmbEval[sapply(accglm.kenya, is.numeric)])
ensmbEval

## Save the ensemble model
#saveRDS(ensMod.genus,'ensMod.genus')


#-------------------------------------------------------------------------------
## Generate an area of applicability to determine uncertinity the spatial predictions
## Area of Applicability Estimation using the ensemble.novel function of the
## BiodiveristyR package

install.packages('BiodiveristyR')
library(BiodiveristyR)

## Subset the areas of reference from the variables used in the prediction
sel.dat.pal <- subset(d, select = c('maxSkinTEMP17', 'minSkinTEMP17', 'minLSTDay17', 
                                    'sumSurfRunoff17', 'maxNDVI17', 'maxSbSMOIST17', 'minSSMOIST17', 
                                    'maxLSTDay18', 'minLSTNight18', 'sumSurfRunoff18', 'minNDVI18', 
                                    'minSbSMOIST18', 'maxSSMOIST18', 'minSSMOIST18', 'maxLSTDay19', 
                                    'maxNDVI19', 'minNDVI19', 'medSbSMOIST19', 'minLSTNight20', 'maxNDVI20', 
                                    'minSbSMOIST20', 'Elevation', 'sand_combat', 'NewTWI', 'Kenya_Human_Pop_2020_1km', 
                                    'Cattle_Density_2020', 'Sheep'))



View(sel.dat)

# Compute the novel envtal conditions (areas where some of the envtal conditions are 
# outside the range of envtal conditions of reference)
# Create the novel object which provides the minima-maxima reference values used 
# in training the model
novel.test.pal <- ensemble.novel.object(sel.dat.pal, name="noveltest")
novel.test.pal

# Create the map with novel conditions using the stacked predictor variables (raster files)
# and the novel object with the reference conditions used in model training
novel.raster.pal <- ensemble.novel(x = Var.stack, novel.object = novel.test.pal)
novel.raster.pal

# Project the data to WGS coordinate reference system
#new_crs <- 'EPSG: 4326'
#novel.rast.pal.Prj <- projectRaster(novel.raster.pal, crs = CRS(new_crs))

# Plot the novel conditions map (darkgreen = Non-novel, gray = novel)
plot(novel.raster.pal, col=c("darkgreen","gray"),plg=list(x="topleft",box.col="black",bty="o", bg = 'white', main="AOA"))


# Select the coordinates columns of the occurrence records and export as a csv
# This to use this information to overlay the novel map and check if all the 
# records are within the non-novel zone, as expected

ptdF <- subset(d, select = c('LONG', 'LAT'))
head(ptdF)

write.csv(ptdF, 'PallidipesTsetsePoints.csv')

# Save the novel map for further visualization
writeRaster(novel.raster.pal, 'Ensemble-AOA_Novel.tif', overwrite = T)

#-------------------------------------------------------------------------------

## Visualizing the dynamic spatial objects 

# Define the color palette
cl <- colorRampPalette(c('white','white','white','darkseagreen','yellow','orange','red'))

#projection(ensMod) <- projection(kenya)

plot(ensMod, col = cl(200))

mapview(ensMod, col = cl(200), maxpixels =  967430)
mapview(glmkenya.pred, col.regions = cl(200), maxpixels = 3869720)
ncol(glmkenya.pred)
nrow(glmkenya.pred)
install.packages('stars')
library(stars)

library(utils)

m = mapview(ensMod, col = cl(200))
html_fl = tempfile(fileext = ".html")
png_fl = tempfile(fileext = ".png")

## create standalone .html file
mapshot(m, url = html_fl)
browseURL(html_fl)

#-------------------------------------------------------------------------------
## Save the model outputs to a director for further visualization
  # Select the individual models from the predictions
glm.mod <- ensmbGenus.pred$sp_1.m_glm.re_cros
svm.mod <- ensmbGenus.pred$sp_1.m_svm.re_cros
max.mod <- ensmbGenus.pred$sp_1.m_maxent.re_cros
rf.mod <- ensmbGenus.pred$sp_1.m_rf.re_cros
Plot(rf.mod) # to confirm and assess the model

## Write each raster file to a directory
writeRaster(rf.mod, filename = '4-RF-Glossina-Model.tif', overwrite=TRUE)

### Generate optimal model thresholds based on validation dataset ####
set.seed(1234) # For reproducibility
d$randomocc<-runif(length(d$Occ)) # Generate uniform random points between 0 and 1 attachin a new column

####Divide training and testing
trainKenya<-d[d$randomocc<0.70,] ## the training database (70%)
evaluKenya<-d[d$randomocc>=0.70,] ## the evaluation database (30%)

evaluKenya$p.logisticKenya<-predict(sdmEnsmb.genus ,type="response",newdata=evaluKenya) # use 'response' to generate probabilities
                                                                                        # add a new column "p.logisticKenya" to evaluKenya
pa.validateKenya<-data.frame(length(evaluKenya$Occ),
                             PA=evaluKenya$Occ,
                             logistic=evaluKenya$p.logisticKenya) # Evaluates how well the model predicts

## Use the optimal.thresholds function to generate classification thresholds
optKenya<- as.data.frame(optimal.thresholds(pa.validateKenya))
data.frame(ID=optKenya[,1], Means=rowMeans(optKenya[,-1]))



#-----------Generate Variable Importance Plots----------------------------------------------------------------------------
## Variable importance
sdm.VIP <- vi(
  ensmbGenus.pred,
  method = "permute", #c("model", "firm", "permute", "shap"),
  feature_names = NULL,
  abbreviate_feature_names = NULL,
  sort = TRUE,
  decreasing = TRUE,
  scale = TRUE,
  rank = TRUE,
  train = 'Var.stack'
)

#-------------------------------------------------------------------------------
## Creating a correlation plot for non-collinear variables

install.packages('car')
library(car)
d.sel<- d[,c('maxNDVI17','minNDVI17','minSSMOIST17',
             'medSSMOIST17','maxLSTDay18','minLSTDay18','minLSTNight18',
             'maxNDVI18','minNDVI18','medNDVI18','minSSMOIST18',
             'minLSTDay19','minLSTNight19','maxNDVI19','medSbSMOIST19',
             'minLSTDay20','maxNDVI20','sand_combat','NewTWI',
             'Kenya_Human_Pop_2020_1km','Cattle_Density_2020','Sheep',
             'Goat','Elevation')]


corrplot(d.sel, is.corr = F, method= 'square', type = 'lower', order = 'original' )
corrplot(d.sel)
corrplot(cor(d.sel, method = 'number', type = 'lower'))
corrplot(cor(d.sel, method = 'kendall'), method = 'square', type = 'lower', 
         tl.col = 'black', tl.cex = 0.9, cl.cex = 0.8, cl.ratio = 0.25, cl.offset = 0.5, 
         number.cex = 0.5, order = 'hclust')


#--------------------------------END--------------------------------------------