作者简介： 本文作者系大学统计学专业教师，多年从事统计学的教学科研工作，在随机过程、统计推 断、机器学习领域有深厚的理论积累与应用实践。个人主页

背景介绍

本文使用物种分布模型，对分布在美国俄亥俄州的两种响尾蛇类型：木纹响尾蛇(Timber Rattlesnake) , 森林响尾蛇(Crotalus horridus)建立地理分布模型。在分析模型之后，我们进一步讨论哪些气候变量对于这些物种的分布范围影响最大。我们也将预测2070年该物种的分布变化情况。

环境与物种数据

我们使用R语言实现代码。首先，加载必需的库。

library(raster)
library(rgdal)
library(maps)
library(mapdata)
library(dismo) 
library(rJava) 
library(maptools)
library(jsonlite)
1
2
3
4
5
6
7
8

建模使用的环境数据来自Worldclim. 我们使用2.5 arc seconds 经纬度当前与未来数据。

currentEnv=getData("worldclim", var="bio", res=2.5)
1

为了估计气候改变，我们使用RCP8.5, 假设温室气体排放持续增加到2100年。我们下载2070年的预测。

futureEnv=getData('CMIP5', var='bio', res=2.5, rcp=85, model='HE', year=70)
names(futureEnv)=names(currentEnv)
1
2

由于气候变量较多，有些彼此相关性强，我们在这里使用简单的变量。

currentEnv=dropLayer(currentEnv, c("bio2", "bio3", "bio4", "bio10", "bio11", "bio13", "bio14", "bio15"))
futureEnv=dropLayer(futureEnv, c("bio2", "bio3", "bio4", "bio10", "bio11", "bio13", "bio14", "bio15"))
1
2

我们保留的预测变量是：

接着，我们需要物种观测记录。这些记录来自Global Biodiversity Information Facility(GBIF). 使用dismo包的函数gbif()可以直接下载数据。

rattler=gbif('crotalus','horridus')
1

下载数据后，需要做初步的数据处理。例如，去掉缺失项、重复项。

rattler=subset(rattler, !is.na(lon) & !is.na(lat)) # the ! symbol can be used to mean 'not'
rattlerdups=duplicated(rattler[, c("lon", "lat")])
rattler <-rattler[!rattlerdups, ]
1
2
3

现在，让我们看一下初步的物种分布图层。我们在观测范围扩展1度。

data(wrld_simpl)
plot(wrld_simpl, xlim=c(min(rattler$lon)-1,max(rattler$lon)+1), ylim=c(min(rattler$lat)-1,max(rattler$lat)+1), axes=TRUE, col="light yellow")
points(rattler$lon, rattler$lat, col="orange", pch=20, cex=0.75)
1
2
3

rattler <- rattler[rattler$lon < -40 & rattler$lat > 25 , ]
map('worldHires', xlim=c(min(rattler$lon)-1,max(rattler$lon)+1), ylim=c(min(rattler$lat)-1,max(rattler$lat)+1), fill=TRUE, col="light yellow")
points(rattler$lon, rattler$lat, col="orange", pch=20, cex=0.75)
1
2
3

model.extent<-extent(min(rattler$lon)-10,max(rattler$lon)+10,min(rattler$lat)-10,max(rattler$lat)+10)
modelEnv=crop(currentEnv,model.extent)
modelFutureEnv=crop(futureEnv, model.extent)
plot(modelEnv[["bio1"]]/10, main="Annual Mean Temperature")
map('worldHires',xlim=c(min(rattler$lon)-10,max(rattler$lon)+10), ylim=c(min(rattler$lat)-10,max(rattler$lat)+10), fill=FALSE, add=TRUE)
points(rattler$lon, rattler$lat, pch="+", cex=0.2)
1
2
3
4
5
6

建模与评价

我们使用最大熵(Maximum Entropy, Maxent))算法拟合SDM. Maxent定义了最优的环境预测变量。

rattlerocc=cbind.data.frame(rattler$lon,rattler$lat) #first, just make a data frame of latitudes and longitudes for the model
fold <- kfold(rattlerocc, k=5) # add an index that makes five random groups of observations
rattlertest <- rattlerocc[fold == 1, ] # hold out one fifth as test data
rattlertrain <- rattlerocc[fold != 1, ] # the other four fifths are training data
rattler.me <- maxent(modelEnv, rattlertrain) #note we just using the training data
plot(rattler.me)
1
2
3
4
5
6

rattler.pred <- predict(rattler.me, modelEnv)
plot(rattler.pred, main="Predicted Suitability")
map('worldHires', fill=FALSE, add=TRUE)
points(rattler$lon, rattler$lat, pch="+", cex=0.2)
1
2
3
4

bg <- randomPoints(modelEnv, 1000)
e1 <- evaluate(rattler.me, p=rattlertest, a=bg, x=modelEnv)

plot(e1, 'ROC')

1
2
3
4
5

rattler.2070 = predict(rattler.me, modelFutureEnv)
plot(rattler.2070, main="Predicted Future Suitability")
map('worldHires', fill=FALSE, add=TRUE)
points(rattler$lon, rattler$lat, pch="+", cex=0.2)
1
2
3
4

rattler.change=rattler.2070-rattler.pred
plot(rattler.change, main="Predicted Change in Suitability")
map('worldHires', fill=FALSE, add=TRUE)
points(rattler$lon, rattler$lat, pch="+", cex=0.2)
1
2
3
4

rattlerChangePoints = extract(rattler.change, rattlerocc)
hist(rattlerChangePoints, main="")
abline(v=0, col="red")
1
2
3

完