Species Distribution Models
- 1.写在前面
- 2.物种分布模型介绍
- 3.输入数据准备及预处理
- 3.1.如何从GBIF网站上获取分布点数据（基于rgbif包）
- 3.2.分布点稀疏处理（基于spThin函数）
- 3.3.如何获取环境变量数据（基于getData函数）
- 3.4.环境变量多重共线性分析
- 4.Biomod2模型构建（v.3.5.1）
1.写在前面
我为什么想要写这个内容呢？ 其实我已经想了很久了，但是一直没有把所有的代码给理清楚，此外这将是我之后研究的主要方向。我对这个方向已经有了初步的了解，看到的满是绝望和不甘。就我个人看法，这个方向已经过了它的热潮，目前正处于低谷期，无人问暇。因此，今年上半年我曾经有一段时间非常痛苦，辗转反侧，我不知道方向还有研究的价值或者必要吗？已经有很多人做过这个方向了，我再来做还有用吗？能提升吗？实际上，可能大多数人只是把这个作为一个应用的工具（水论文），而不是深入的探究，因此我认为只是做应用层面的，这个方向已经饱和了，而我们想要继续深挖下去，就得更加了解其原理和意义，与实际情况相结合才能用更好的产出（如生物多样性保护等）。目前来说物种分布模型（Species Distribution Models，SDMs）已经有很多人在研究了，包括国内的中科院动物所乔慧捷老师、中科院青藏高原所的郭彦龙老师、以及法国科学院的Wilfried Thuiller教授（Biomod2包的发明者）。
现在，我把我的所有SDMs的构建流程进行记录，希望对相关研究方向的学者有所助益。
2.物种分布模型介绍
目前常用的SDMs有单一模型和集成模型，考虑到单个模型的优点和缺点，使用整体模型可能更安全。其中，Biomod2平台是目前常用的集成平台，它包括了十个不同类型的模型：
- 三个回归模型：广义线性模型（generalized linear model, GLM）;多元自适应回归样条（multivariate adaptive regression splines, MARS），广义加法模型（generalized additive models, GAM）；
- 五个机器学习模型：人工神经网络（artificial neural network, ANN），最大熵（maximum entropy, MaxEnt），随机森林（random forest, RF）和广义助推模型（generalized boosting model, GBM），分类树分析（classification tree analysis, CTA）、
- 一个分类模型：灵活判别分析（flexible discriminant analysis, FDA）
- 一个范围包络：表面分布区分室模型（surface range envelope, SRE）
如果你想继续了解SDMs的研究方向和基本知识，我推荐你阅读以下几篇 综诉类文献：
[1] 郭彦龙,赵泽芳,乔慧捷,等.物种分布模型面临的挑战与发展趋势[J].地球科学进展,2020,35(12):1292-1305.
[2] 刘晓彤,袁泉,倪健.中国植物分布模拟研究现状[J].植物生态学报,2019,43(04):273-283.
[3] 许仲林,彭焕华,彭守璋.物种分布模型的发展及评价方法[J].生态学报,2015,35(02):557-567.
3.输入数据准备及预处理
首先我们需要知道构建SDMs需要哪些要素：物种分布点数据（Species occurrence data）和环境变量数据（Environmental data），分布点数据可从几个地方进行获取，首先是全球生物多样性网络（GBIF）、中国数字植物标本馆（CVH）、NSII-中国国家标本资源平台以及自己的野外调查数据。
3.1.如何从GBIF网站上获取分布点数据（基于rgbif包）
##需要改四个地方的物种名 ----
setwd("D:/SDMs/Thesis/Data/")
library(mapdata)
library(sp)
library(sf)
library(raster)
library(ggplot2)
library(CoordinateCleaner)
library(rgbif)
library(dplyr)
library(beepr)# 这个包是用来进行声音提示的
#第一步：使用''rgbif'' 包从 GBIF 中搜索与目标物种小叶栎(Quercus chenii)相关的数据------
key <- name_suggest(q="Quercus chenii", rank='species')$data$key[1]
key #物种识别号：2879172
occ_count(taxonKey=key, georeferenced=TRUE,
basisOfRecord = "PRESERVED_SPECIMEN"
)
# 266，目前，GBIF 数据库中有 266 条基于博物馆的记录是具有经纬度信息的。
# 需要注意的是，因为数据库一直在更新，这个数字可能会变动，
# 因此使用本流程过程中可能会有部分数值有一些小范围的变动。
data <- occ_search(scientificName = "Quercus chenii", limit= 10000,
hasCoordinate = TRUE, basisOfRecord= "PRESERVED_SPECIMEN"
# ,country = "CN" # 将物种限定在中国境内
) #这一步需要 1-2 分钟完成，具体看网速和物种分布数量。
beepr::beep(8) # 上面代码运行结束会有声音提示
str(data$data)
data$meta # 查看有多少个物种分布点
names(data$data) #检查数据的列名，一共有162列数据
head(data$data[,c("countryCode","country")]) # 例如：检查 54 列和 57 列的列名
datasel<-data$data %>%
dplyr::select(species, decimalLongitude,
decimalLatitude,country,coordinateUncertaintyInMeters) %>%
rename(Species = species, Long = decimalLongitude, Lat =
decimalLatitude, Uncertain=coordinateUncertaintyInMeters) %>%
mutate (Species = as.factor(Species), country=
as.factor(country))
datasel # 现在有了一个包含266 行、5 列数据：Species、Long、Lat、country、Uncertain
datasel %>% group_by(country) %>% summarize(total= n()) %>%
arrange(desc(total)) # GBIF 检索的数据集中每个国家的记录数量
# write.csv(datasel, "Quercus acutissima Carruth.csv")
#第二步：使用 CoordinateCleaner 包删除有问题的分布记录--------------------------
map('worldHires', col=1:10)
# plot(world_map, col="tan2",bg="lightblue", border= NA) #图示地图。
points(datasel$Long, datasel$Lat, bg="red", col='white', pch=21, cex=1.0)
#将从 GBIF 下载的数据映射到地图上 (图 9)。
datasel %>% filter(Long==0, Lat==0) # 统计经纬度等于零的记录
problem_records<- clean_coordinates(datasel, lon = "Long", lat = "Lat",
species = "Species",
tests = c("capitals", "centroids",
"equal","institutions",
"zeros", "seas"))
# 检测记录点坐标是否围绕首都、国家的中心
# 是否落入海洋，为零，或在饲养动物的博物馆(机构)周围。
summary(problem_records) # 2个有问题的分布点被发现
datasel_problematic <- datasel[which(problem_records$.summary==
"FALSE"),] #只选择有问题的点。
# 查看有问题分布点的位置
map('worldHires', col=1:10)
points(datasel_problematic$Long, datasel_problematic$Lat, bg="red",
col='white', pch=22, cex=0.8) # 将有问题的点映射到地图上
# 从数据集中排除所有有问题的记录
datasel_clean <- datasel[which(problem_records$.summary== "TRUE"),]
str(datasel_clean) # 现在还剩下264个正常分布点数据
# 删除坐标中具有较高不确定性的记录。
# 看看有多少个体有这样的信息，并把值的范围图像化
summary(is.na(datasel_clean$Uncertain))
hist(datasel_clean$Uncertain/1000, breaks = 10,
main="Coordinate uncertainy in km", xlab= "Uncertainty in km")
# 注意，GBIF 提供的信息是以米为单位
# 但是为了评估不确定性，使用千米为单位更加方便。 (图 13)。
# 这个命令中用于选择“信息不确定”的点，或者是数值<10 km 以下的点。
datasel_cleanF<-datasel_clean %>% filter(is.na(Uncertain)| Uncertain/1000 < 10)
#最后一步是删除重复的记录。
datasel_cleanF <- datasel_cleanF %>% distinct(Long, Lat,.keep_all=T)
str(datasel_cleanF) # 图 14。
#以下代码可以图示这个最终数据集。
map('worldHires', col=1:10)
points(datasel_cleanF$Long, datasel_cleanF$Lat, bg="blue", col='white',
pch=21, cex=0.5) # 图 15
# 最后记得将数据保存下来
write.csv(datasel_cleanF, "Quercus chenii_gbif.csv")
如果你想要从CVH网站上获得分布点数据，就需要自己去查看里面的标本记录，然后手动记录下来，这个过程可能比较繁琐，目前没有相关的包能够获取上面的记录。
3.2.分布点稀疏处理（基于spThin函数）
为了防止数据空间自相关性对模拟结果的影响，通常需要对分布点数据进行稀疏处理，如10km范围内只保留一个分布点数据。目前常用的方法是在ArcGIS中安装SDMtoolbox （v.2.5） 小插件，并使用其中所提供的 “Spatially Rarefy Occurence Data” 功能进行处理。但是采用这样的方式有时候会报让人摸不着头脑的错误，所以我提供了一个使用R语言进行稀疏处理的方法，如下：
## thin----------------------------------------------
setwd("C:\\Users\\YP\\Desktop\\SDMData")
# library(wallace)
library(spThin) # v.0.2.0
library(readxl)
library(mapdata)
# run_wallace() ##运行run后可直接进入到浏览器界面
QuercusL <- read.csv("Quercus chenii_gbif.csv", sep = ",")
head(QuercusL)
summary(QuercusL)
# remove rows with duplicate coordinates
occs.dups <- duplicated(QuercusL[c('Long', 'Lat')])
occs <- QuercusL[!occs.dups,]
# remove NAs
occs <- occs[complete.cases(occs$Long, occs$Lat), ]
summary(occs)
output <- spThin::thin(occs, 'Lat', 'Long', 'Species',
thin.par = 10, reps = 100, # 稀疏距离参数thin.par设置为10km，并重复运行100次
locs.thinned.list.return = TRUE,
write.files = FALSE,
verbose = TRUE)
plotThin(output)
# 基于最优迭代构建最优输出
maxThin <- which(sapply(output, nrow) == max(sapply(output, nrow)))
# if there's more than one max, pick the first one
maxThin <- output[[ifelse(length(maxThin) > 1, maxThin[1], maxThin)]]
# subset occs to match only thinned occs
occs <- occs[as.numeric(rownames(maxThin)),]
write.csv(occs,"./Quercus chenii_gbif_occs.csv")
# plot
map('worldHires', 'China')
points(occs$Long, occs$Lat, bg="red", col='white', pch=21, cex=1.0)
分布点空间位置：
3.3.如何获取环境变量数据（基于getData函数）
目前，常用的大尺度空间数据主要由WorldClim提供，我之前写了一篇关于如何下载该数据的方法，感兴趣的可以看看：Worldclim（v1.4、v2.1）数据集使用介绍
在这里我可以提供另一种方法：使用R语言下载数据。代码如下：
##使用"Raster"包从 WorldClim 下载环境数据--------------------------------
# 目前该函数只提供了下载2.5, 5, 和 10 min的空间分辨率下的数据，且来源于worldclim1.4（CMIP5）
bioclim<-getData('worldclim', var="bio", res=5) #把 19 个5min分辨率的气候数据图层下载到本地
Future_bioclim2<-getData('CMIP5', var='tmin', res=10, rcp=85, model='AC', year=70) # 获取未来气候数据
names(bioclim) #查看存储在"bioclim"对象中的不同气候变量层的名称
我可能不太建议使用R代码来下载数据，因为这可能受到网速的限制，导致下载突然中断，所以直接到官网里面下载，这样也很方便。
当然还可以使用CHELSA数据集，这是由瑞士森林雪和景观研究所（WSL）提供的。现在也有很多研究使用该数据：
3.4.环境变量多重共线性分析
在构建SDMs时，会存在不同环境变量对分布分布点的贡献程度相似，即可看做这几个变量对物种分布的作用相同，因此为了降低模型拟合的不确定性，就需要从这几个变量中选择一个贡献度最高，最关键的环境变量，因此就需要使用多重共线性分析。常用的方法为Pearson相关性系数以及方差膨胀因子（VIF）。其中Pearson相关性系数和VIF计算方法如下：
library(raster)
library(pheatmap)
library(usdm)
setwd("C:/Users/YP/Desktop/SDMData")
quercus <- read.csv("Quercus chenii_gbif_occs.csv", sep = ",")
summary(quercus)
current_folder1 <- "C:/Users/YP/Desktop/Data/"
# 气候因子1
current_bio <- stack(paste0(current_folder1,
c("wc2.1_30s_bio_1","wc2.1_30s_bio_2",
"wc2.1_30s_bio_3","wc2.1_30s_bio_4",
"wc2.1_30s_bio_5","wc2.1_30s_bio_6",
"wc2.1_30s_bio_7","wc2.1_30s_bio_8",
"wc2.1_30s_bio_9","wc2.1_30s_bio_10",
"wc2.1_30s_bio_11","wc2.1_30s_bio_12",
"wc2.1_30s_bio_13","wc2.1_30s_bio_14",
"wc2.1_30s_bio_15","wc2.1_30s_bio_16",
"wc2.1_30s_bio_17","wc2.1_30s_bio_18",
"wc2.1_30s_bio_19"
),
".tif")
)
# 修改变量名
new_names1 <- paste0("bio", 1:19)
names(current_bio) <- new_names1
# 地形因子2
current_topo <- stack(paste0(current_folder1,
c("DEM_2023", "Aspect", "Slope"),
".tif")
)
new_names2 <- c("Elevation", "Aspect", "Slope")
names(current_topo) <- new_names2
# 土壤因子3
current_soil <- stack(paste0(current_folder1,
c("BD1","PH1", "TK1",
"TN1","TP1","soil_1km1"),
".tif")
)
new_names3 <- c("BD","pH","TK","TN","TP","SoilType")
names(current_soil) <- new_names3
# 将以上不同变量数据合并为一个叠加图层
(current <- stack(current_bio, current_topo, current_soil))
# 进行坐标转换，将分布点数据转换为环境数据空间投影格式，即WGS1984
longlat <- as.data.frame(quercus[, c(4:5)]) # Long Lat
coordinates(longlat) <- c("Long", "Lat") # 这里只能是Long在前，Lat在后，不然下一行代码会报错！！！
proj4string(longlat) <- CRS("+proj=longlat +datum=WGS84")
UTM_Proj <- CRS(st_crs(current)$proj4string)
xy <- spTransform(longlat, UTM_Proj)
# 检查坐标系是否一致
if (st_crs(current)$proj4string == st_crs(xy)$proj4string) {
print("坐标系转换结果：坐标系相同")
} else {
print("坐标系转换结果：坐标系不同")
}
# 从所有的记录中提取环境因子对应值
data <- raster::extract(current, xy)
data_na <- na.omit(data)
# 构建pearson相关分析(pearson>0.80)
cor_matrix <- cor(data_na, method = "pearson") # R在带stats包
cor_matrix
# dev.off() # 聚类热图绘制报错时可以使用，这可以清除之前可能干扰当前绘图的图形设置
pheatmap(cor_matrix,cluster_cols = F,display_numbers = T,main = "环境因子相关性分析")
# 计算方差膨胀因子VIF
vifstep(data_na) # usdm包:vifstep()
多重共线性诊断VIF： 一般文章的材料方法中对VIF的限定不同，有的限定在5以内，有的限定在10以内，也有的限定在20以内，这个根据自己的判断来定，没有十分明确的数值，数值越小越好，一般建议保留在VIF<10的变量。
文章：
-
1. Analysis of the distribution pattern of Chinese Ziziphus jujuba under climate change based on optimized biomod2 and MaxEnt models. Ecological Indicators：
-
2. Predicting the invasive trend of exotic plants in China based on the ensemble model under climate change: A case for three invasive plants of Asteraceae. Science of The Total Environment:
-
3. Projected degradation of quercus habitats in southern China under future global warming scenarios. Forest Ecol. Manag.:
-
4. Incorporating eco-evolutionary information into species distribution models provides comprehensive predictions of species range shifts under climate change. Science of The Total Environment:
Pearson相关性热图：
在得到Pearson相关性系数之后，如果我们要去除Pearson>0.80的变量，那么我们应该选择这两个变量中好解释的变量（不要为难自己，给自己挖坑），或者是你认为更加重要的变量（你就是这个领域的专家，你的知识就是最正确的）。
4.Biomod2模型构建（v.3.5.1）
目前最新的版本应该是v.4.2-1，最新版可能使用起来更加流畅和友好，但是我比较习惯使用v.3.5.1版本构建模型，这个是已经成熟使用过的版本，可能报错少一些，运行更加稳定。
