【翻译】Species distribution modeling 6 Medeling methods
6.1 实例中使用的算法和数据的类型
物种分布建模已经用过了大量的算法。它们可以归类为‘框架（profile）’，‘回归（regression）’和‘机器学习（machine learning）’方法。框架方法只考虑‘存在（presence）’数据，而不是不存在或背景数据。回归和机器学习方法使用存在数据和不存在或背景数据。回归与机器学习方法之间的区别不太清晰（sharp），但它也许仍可以作为一种模型归类的途径。可以做出的另一个区别是仅存在和存在-不存在模型之间的。框架方法通常是仅存在的，其他方法也可以这样，屈居于它们是否使用调查-不存在（survey-abence）或者是假设不存在（pseudo-absence）/背景（background）数据。一个完全不同的模型类型由仅仅或主要使用已知发生的地理位置的模型组成，并且不依赖于这些位置处的预测变量的值。我们将这些模型称为‘地理模型（geographical models）’。我们在下方导论这些不同类型模型的例子。
首先我们重建到现在为止使用过的数据。
library(dismo) library(maptools) data(wrld_simpl) predictors <- stack(list.files(file.path(system.file(package="dismo"), 'ex'), pattern='grd$', full.names=TRUE )) file <- file.path(system.file(package="dismo"), "ex/bradypus.csv") bradypus <- read.table(file, header=TRUE, sep=',') bradypus <- bradypus[,-1] presvals <- extract(predictors, bradypus) set.seed(0) backgr <- randomPoints(predictors, 500) absvals <- extract(predictors, backgr) pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals))) sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals))) sdmdata[,'biome'] <- as.factor(sdmdata[,'biome'])
我们将会使用同样的数据去演示所有模型，除了一些不能使用分类变量的模型。因此对于这些模型，我们从预测器堆栈（predictors stack）中删除分类变量。
pred_nf <- dropLayer(predictors, 'biome')
我们使用Bradypus得到物种的存在数据。首先我们制作一个训练（training）和一个测试（testing）集。
set.seed(0) group <- kfold(bradypus, 5) pres_train <- bradypus[group != 1, ] pres_test <- bradypus[group == 1, ]
为了加快处理速度，让我们将预测限制到一个更为限制性的区域（定义成一个矩形区域）。
ext <- extent(-90, -32, -33, 23)
训练背景数据和测试集。RasterStack中的第一个图层被用作‘蒙板（mask）’。可以确保随机点仅在栅格（rasters）的空间范围内发生，单元格不是NA且每个单元格只有一个不存在点。这里我们进一步将背景点限制在我们制定范围‘ext’的12.5%以内。
set.seed(10)
backg <- randomPoints(pred_nf, n=1000, ext=ext, extf = 1.25)
colnames(backg) = c('lon', 'lat')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]
r <- raster(pred_nf, 1)
plot(!is.na(r), col=c('white', 'light grey'), legend=FALSE)
plot(ext, add=TRUE, col='red', lwd=2)
points(backg_train, pch='-', cex=0.5, col='yellow')
points(backg_test, pch='-', cex=0.5, col='black')
points(pres_train, pch= '+', col='green')
points(pres_test, pch='+', col='blue')
6.2 框架方法
这里描述了三种方法，有Bioclim、Domain和Mahal。这些方法在dismo包中被实现，而使用这3种模型的步骤是相同的。
6.2.1 Bioclim
BIOCLIM算法被广泛应用于物种分布建模。BIOCLIM是一个传统的‘气候包络模型（climate-envelope-model）’（Booth et al., 2014）。虽然它通常不如其他一些建模方法那么好（Elith et al., 2006）， 特别是在气候变化的背景下（Hijmans and Graham, 2006），尽管仍被一些其他的问题所包围，它仍然被使用着。因为算法易于理解，从而有助于教学物种分布建模。BIOCLIM算法通过将任何地点的环境变量的值与在已知发生地点（‘训练场所（training sites’））处的值的百分位数分布进行比较来计算地点的相似度。越接近第五百分位数（中位数），地点越适合。分布的结尾没有区别，也就是说，10百分位数被视为相当于90百分位数。在‘dismo’的实施中，上尾值被转换为下尾，并且使用所有环境变量的最小百分位数得分（也就是说，BIOCLIM使用像Liebig定律的最低限度的方法）。该值从1中被减去，然后乘以2，使结果在0和1之间。这种缩放方法的原因是为了使结果更像其他分布建模方法，从而更容易去解释。值1将很少被观察到，因为需要具有所考虑的所有变量的训练数据（training data）的中值的位置。值0是非常普遍的，因为它被分配给对于至少一个超出百分位数分布（训练数据的范围）的环境变量的值的所有单元格。
早些时候，我们使用data.frame拟合了一个Bioclim模型，将每个行代表已知物种存在位置的环境数据。在这里，我们只使用预测因子和出现点（函数会为我们提取）拟合一个bioclim模型。
我们以类似的方式评估模型，通过提供存在和背景（不存在）点，模和一个RasterStack。
e <- evaluate(pres_test, backg_test, bc, pred_nf) e ## class : ModelEvaluation ## n presences : 23 ## n absences : 200 ## AUC : 0.6957609 ## cor : 0.17802 ## max TPR+TNR at : 0.08592151
找到阈。
tr <- threshold(e, 'spec_sens') tr ## [1] 0.08592151
然后我们使用具有预测变量的RasterStack对RasterLayer作出预测。
pb <- predict(pred_nf, bc, ext=ext, progress='') pb ## class : RasterLayer ## dimensions : 112, 116, 12992 (nrow, ncol, ncell) ## resolution : 0.5, 0.5 (x, y) ## extent : -90, -32, -33, 23 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 ## data source : in memory ## names : layer ## values : 0, 0.7096774 (min, max) par(mfrow=c(1,2)) plot(pb, main='Bioclim, raw values') plot(wrld_simpl, add=TRUE, border='dark grey') plot(pb > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+')
请注意预测函数中参数的顺序。在上方的例子中，我们使用了predict(pred\_nf, bc)（首先是RasterStack，然后是模型对象），这比predict(bc,pred_nf)（首先是模型，然后是RasterStack）效率低一些。我们使用这种顺序的原因是它适用于所有模型，而其他选择只适用于在dismo包中定义好的模型，比如Bioclim、Domain和Maxent，而不是定义在其他包中的模型（random forest、boosted regression trees、glm等）。
6.2.2 Domain
