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