【R语言实例】物种分布模型介绍

Goodsta

关注

已于 2024-09-08 17:54:09 修改 · 1.3w 阅读

原创

本文为博客 VIP 文章，开通 VIP 后可阅读全文

开通 VIP

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

1. 背景知识

物种分布模型(species distribution models, SDM)是数量生态学里的一个流行的分析工具。SDM普遍使用数理模型评估全球生态变化对于物种迁移的潜在影响。它的优势体现在：有很多容易使用的软件工具与参考指南，对数据的要求较低。

1.1 SDM 输入

SDM要求输入有地理坐标参考的生物多样性观测。比如，个体坐标、物种的presence, 物种数量、物种的richness, 响应变量等。此外，地理图层、环境信息也可以作为SDM的输入。这些信息通常以数码格式保存。

1.2 在线数据库

GBIF the Global Biodiversity Information Facility
OBIS Ocean Biodiversity Information System
Movebank animal tracking data
WorldClim global climate data

利用这些数据，我们可以建立统计机器学习模型，描述特定地点的生物多样性观测与气候条件的关系，然后把模型投射到可利用的环境图层上，在时空上做预测。

在这里插入图片描述

2. SDM 实例

2.1 一个例子：Ring Ouzel

环颈鸫(Ring Ouzel), 鸫属的一种鸟类，生活在欧洲，与乌鸫有亲缘关系。其成年雄性通体黑色，胸部有白色月牙标记，有黄色的喙。雌性与其相似，但颜色较淡，胸部无白色月牙。

在这里插入图片描述
在这个实例里，我们想评估气候变化对栖息在瑞士的环颈鸫种群的影响。环颈鸫主要栖息地在瑞士北方的阿尔卑斯山脉，对气候变暖比较敏感，另外还有其它的影响指标。上世纪90年代以来，环颈鸫的种群密度已经下降，而在海拔2000米以上的群体仍然维持稳定。在邻近国家，种群假定是稳定的。
瑞士鸟类协会收集了两类环颈鸫分布数据。这里，我们研究气候对环颈鸫分布的影响。为此，选取了5个最重要的预测变量。

bio5 = maximum temperature of the warmest month,
bio2 = mean diurnal range,
bio14 = precipitation of driest month,
std = standard deviation of vegetation height,
rad = annual total radiation.

为了理解不同的统计模型/算法的区别，我们拟合简单的广义线性模型与随机森林。模型输出的是当前与未来分布概率图，潜在的分布二值图谱。

2.2 数据准备

在这一步，收集处理真实的环颈鸫生物多样性与环境数据。对于生物多样性数据，物种的presence信息可以直接观测到，而absence信息作为presence数据的对比数据很难获得。这种情况下，我们需要充分的background数据，或者pseudo-absence数据。准备好后，首先导入数据。

avi_dat <- read.table('data/Data_SwissBreedingBirds.csv', header=T, sep=',')

summary(avi_dat)

该数据框由56个鸟类物种的2,535条presence-absence信息记录、52个环境预测变量组成。其中，数据的70%用作单一物种分布建模，其余30%用作检验预测。下面，我们缩减数据框到相关的变量列。

avi_cols <- c('Turdus_torquatus', 'bio_5', 'bio_2', 'bio_14', 'std', 'rad', 'blockCV_tile')

avi_df <- data.frame(avi_dat)[avi_cols]

summary(avi_df)

在这里插入图片描述
为了定位对应当前气候的物种分布，以及未来气候下的潜在分布图，我们使用WorldClim数获得生物气候变量值。该数据可以直接从R环境下载。

library(raster)

# Please note that you have to set download=T if you haven't downloaded the data before:
bio_curr <- getData('worldclim', var='bio', res=0.5, lon=5.5, lat=45.5, path='data')[[c(2,5,14)]]

# Please note that you have to set download=T if you haven't downloaded the data before:
bio_fut <- getData('CMIP5', var='bio', res=0.5, lon=5.5, lat=45.5, rcp=45, model='NO', year=50, path='data', download=F)[[c(2,5,14)]]

我们使用背景mask定位瑞士坐标，这需要重新投射worldclim层。

# A spatial mask of Switzerland in Swiss coordinates
bg <- raster('/vsicurl/https://damariszurell.github.io/SDM-Intro/CH_mask.tif')

# The spatial extent of Switzerland in Lon/Lat coordinates is roughly:
ch_ext <- c(5, 11, 45, 48)

# Crop the climate layers to the extent of Switzerland
bio_curr <- crop(bio_curr, ch_ext)

# Re-project to Swiss coordinate system and clip to Swiss political bounday
bio_curr <- projectRaster(bio_curr, bg)
bio_curr <- resample(bio_curr, bg)
bio_curr <- mask(bio_curr, bg)
names(bio_curr) <- c('bio_2', 'bio_5', 'bio_14')

# For storage reasons the temperature values in worldclim are multiplied by 10. For easier interpretability, we change it back to °C.
bio_curr[[1]] <- bio_curr[