SPATIAL DISTRIBUTION MODELS

Introduction

Install Packages

if (!require("rspat")) remotes::install_github("rspatial/rspat") ###Loading required package: rspat ###Loading required package: terra ###terra 1.7.62 
if (!require("predicts")) install.packages("predicts") ###Loading required package: predicts 
if (!require("geodata")) install.packages("geodata") ###Loading required package: geodata 

call terra and rspat package

library(terra) 
library(rspat) 
bf <- spat_data("bigfoot") 

Observation

dim(bf) 

[1] 3092    3

head(bf) 

• Purpose: Checks data dimensions and previews data.
• Output Use: Understand the structure of your dataset. ### Plot

plot(bf[,1:2], cex=0.5, col="red") 
library(geodata) 
wrld <- geodata::world(path=".") 
bnds <- wrld[wrld$NAME_0 %in% 
               c("Canada", "Mexico", "United States"), ] 
lines(bnds) 

Purpose: Plots species occurrence points. Output Use: Visualizes spatial distribution.

Visualization

wc <- geodata::worldclim_global("bio", res=10, ".") 
plot(wc[[c(1, 12)]], nr=2) 

climate data

bfc <- extract(wc, bf[,1:2]) 
head(bfc, 3) 

remove the first column with the ID

bfc <- bfc[,-1] 

plot the species’ distribution

plot(bfc[ ,"wc2.1_10m_bio_1"], bfc[, "wc2.1_10m_bio_12"], col="red", xlab="Annual mean temperature (.C)", ylab="Annual precipitation (mm)") 

Background data

ext_bf <- ext(vect(bf[, 1:2])) + 1 
ext_bf 

SpatExtent : -157.75, -63.4627, 24.141, 70.5 (xmin, xmax, ymin, ymax)

take 5000 random samples

set.seed(0) 
window(wc) <- ext_bf 
bg <- spatSample(wc, 5000, "random", na.rm=TRUE, xy=TRUE) 
head(bg) 

Plot spatSample

plot(bg[, c("x", "y")]) 

Data

bg <- bg[, -c(1:2)] 

Plot Data

plot(bg[,1], bg[,12], xlab="Annual mean temperature (∞C)", ylab="Annual precipitation (mm)", cex=.8) 
points(bfc[,1], bfc[,12], col="red", cex=.6, pch="+") 

legend("topleft", c("observed", "background"), col=c("red", "black"), pch=c("+", "o"),pt.cex=c(.6, .8)) 

East vs West

#eastern points 
bfe <- bfc[bf[,1] > -102, ] 
#western points 
bfw <- bfc[bf[,1] <= -102, ] 
dw <- rbind(cbind(pa=1, bfw), cbind(pa=0, bg)) 
de <- rbind(cbind(pa=1, bfe), cbind(pa=0, bg)) 
dw <- data.frame(dw) 
de <- data.frame(na.omit(de)) 
dim(dw) ###[1] 6224 20 

[1] 6224   20

dim(de) ###[1] 6866 20 

[1] 6866   20

Fit a model

CART

library(rpart) 
cart <- rpart(pa~., data=dw) 

complexity paramete

printcp(cart) 

Regression tree:
rpart(formula = pa ~ ., data = dw)

Variables actually used in tree construction:
[1] wc2.1_10m_bio_10 wc2.1_10m_bio_12 wc2.1_10m_bio_14 wc2.1_10m_bio_18 wc2.1_10m_bio_19 wc2.1_10m_bio_3 
[7] wc2.1_10m_bio_4 

Root node error: 983.29/6224 = 0.15798

n= 6224 

        CP nsplit rel error  xerror     xstd
1 0.322797      0   1.00000 1.00019 0.019351
2 0.080521      1   0.67720 0.68771 0.019687
3 0.073325      2   0.59668 0.59759 0.015480
4 0.068645      3   0.52336 0.52163 0.015231
5 0.027920      4   0.45471 0.47026 0.014758
6 0.014907      5   0.42679 0.44778 0.015136
7 0.010869      6   0.41188 0.44082 0.015406
8 0.010197      7   0.40102 0.43500 0.015330
9 0.010000      8   0.39082 0.43197 0.015319

plot

plotcp(cart) 

Fit the model

cart <- rpart(pa~., data=dw, cp=0.02) 

plot

library(rpart.plot) 
rpart.plot(cart, uniform=TRUE, main="Regression Tree") 

use the model to show how attractive the climate is for this species

x <- predict(wc, cart) 
x <- mask(x, wc[[1]]) 
x <- round(x, 2) 
plot(x, type="class", plg=list(x="bottomleft")) 

Random Forest

Test and Train Data

set.seed(123) 
i <- sample(nrow(dw), 0.2 * nrow(dw)) 
test <- dw[i,] 
train <- dw[-i,] 
fpa <- as.factor(train[, 'pa']) 

RandomForest model

library(randomForest) ###randomForest 4.7-1.1 ###Type rfNews() to see new features/changes/bug fixes. 
crf <- randomForest(train[, 2:ncol(train)], fpa) 
crf

Call:
 randomForest(x = train[, 2:ncol(train)], y = fpa) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 4

        OOB estimate of  error rate: 7.19%
Confusion matrix:
     0   1 class.error
0 3832 165  0.04128096
1  193 790  0.19633774

variable importance plot

varImpPlot(crf) 

use regression ( tune a parameter )

trf <- tuneRF(train[, 2:ncol(train)], train[, "pa"]) 

Warning: The response has five or fewer unique values.  Are you sure you want to do regression?

mtry = 6  OOB error = 0.05594974 
Searching left ...

Warning: The response has five or fewer unique values.  Are you sure you want to do regression?

mtry = 3    OOB error = 0.05612125 
-0.003065481 0.05 
Searching right ...

Warning: The response has five or fewer unique values.  Are you sure you want to do regression?

mtry = 12   OOB error = 0.05485775 
0.01951734 0.05 

trf 

   mtry   OOBError
3     3 0.05612125
6     6 0.05594974
12   12 0.05485775

mt <- trf[which.min(trf[,2]), 1] 
mt 

[1] 12

tuneRF and values of mt represen

rrf <- randomForest(train[, 2:ncol(train)], train[, "pa"], mtry=mt, ntree=250) 

Warning: The response has five or fewer unique values.  Are you sure you want to do regression?

plot(rrf) 

Predict

Regression

rp <- predict(wc, rrf, na.rm=TRUE)
plot(rp) 

optimal threshold

library(predicts) 
eva <- pa_evaluate(predict(rrf, test[test$pa==1, ]), predict(rrf, test[test$pa==0, ])) 
eva 

@stats

@thresholds

@tr_stats

ROC

plot(eva, "ROC") 

Box plot and predictive value

par(mfrow=c(1,2)) 
plot(eva, "boxplot") 
plot(eva, "density") 

use the “max specificity + sensitivity” threshold

tr <- eva@thresholds 
tr ## max_kappa max_spec_sens no_omission equal_prevalence equal_sens_spec 

## 1 0.4469667 0.3219856 0.00425973 0.1952333 0.2167239
plot(rp > tr$max_spec_sens) 

Classification

classification Random Forest model

rc <- predict(wc, crf, na.rm=TRUE) 
plot(rc) 

Add probability

rc2 <- predict(wc, crf, type="prob", na.rm=TRUE) 
plot(rc2, 2) 

Extrapolation

Distribution of the Eastern species

Plot

our model is really good in predicting the range of the West, but it cannot extrapolate at all to the East.

World Map

window(wc) <- NULL 
pm <- predict(wc, rrf, na.rm=TRUE) 
plot(pm) 
lines(wrld) 

estimate range shifts due to climate change

fut <- cmip6_world("CNRM-CM6-1", "585", "2061-2080", var="bio", res=10, path=".") 
names(fut)

 [1] "bio01" "bio02" "bio03" "bio04" "bio05" "bio06" "bio07" "bio08" "bio09" "bio10" "bio11" "bio12" "bio13"
[14] "bio14" "bio15" "bio16" "bio17" "bio18" "bio19"

names(fut) <- names(wc) 
window(fut) <- ext_bf 
pfut <- predict(fut, rrf, na.rm=TRUE) 
plot(pfut)