(15)生态

发表于 更新于 分类于 ， 阅读次数： Waline： 本文字数： 7.2k 阅读时长 ≈ 26 分钟

本文由SCY翻译自《Geocomputation with R》第十五章

本章通过模拟雾绿洲（也叫lomas）的植被分布，揭示了明显受到水源供应控制的不同植被区域。这个案例研究不仅整合了前几章中的核心观点，还拓展了这些观点，帮助你更全面地掌握使用R进行地理计算的相关技能。

前提条件

本章假设你已经掌握了地理数据分析和处理的基础知识，这在空间数据章至几何操作章中进行了涵盖。本章利用了 GIS 软件的桥梁，以及在GIS桥梁和统计学习章中涉及到的空间交叉验证。本章使用了下列包：

1
2
3
4
5
6
7
8
9
10
11
12
13
library(sf)
library(terra)
library(dplyr)
library(data.table)        # fast data.frame manipulation (used by mlr3)
library(mlr3)              # machine learning (see Chapter 12)
library(mlr3spatiotempcv)  # spatio-temporal resampling 
library(mlr3tuning)        # hyperparameter tuning package
library(mlr3learners)      # interface to most important machine learning packages
library(paradox)           # defining hyperparameter spaces
library(ranger)            # random forest package
library(qgisprocess)       # bridge to QGIS (Chapter 10)
library(tree)              # decision tree package
library(vegan)             # community ecology package

引言

本章通过模拟雾绿洲（也叫lomas）的植被分布，揭示了明显受到水源供应控制的不同植被区域。这个案例研究不仅整合了前几章中的核心观点，还拓展了这些观点，帮助你更全面地掌握使用R进行地理计算的相关技能。

雾绿洲或lomas，特指那些分布在秘鲁和智利沿海沙漠山脉的植被区。在世界的其他角落，如纳米比亚的沙漠和也门与阿曼的海岸，也存在类似的生态系统[@galletti_land_2016]。尽管这些地区年均降水量极低，仅为30-50毫米，但是在南半球的冬季，雾气的沉积显著增加了植物可用的水分。这导致了秘鲁沿海一带南坡的山地呈现出一片翠绿。这种由于南半球冬季冷洋流引发的温度逆转层下形成的雾气，成为了这一生境名称的由来。

每隔数年，厄尔尼诺现象会给这片烈日炙烤下的土地带来一场暴雨，为树苗提供了一个难得的机会：它们有足够的时间发展出能够在随后的干旱环境中生存下去的根系[@dillon_lomas_2003]。

不幸的是，雾绿洲正面临严重的威胁，主要原因是农业活动和人为的气候变化。有关本土植被成分和空间分布的证据能够支持保护剩余雾绿洲碎片的努力[@muenchow_predictive_2013; @muenchow_soil_2013]。本章还演示了如何将前几章介绍的技术应用到一个重要的应用领域：生态学。在这一章中，你将分析秘鲁中北部海岸附近卡斯马的蒙贡山（Mt. Mongón）南坡上的维管植物（主要是指开花植物）的组成和空间分布。


The Mt. Mongón study area, from Muenchow, Schratz, and Brenning (2017).

在2011年的澳大利亚冬季，进行了一次前往蒙贡山（Mt. Mongón）的实地研究，记录了100个随机抽样的4x4平方米区域内所有生存的维管植物[@muenchow_predictive_2013]。
该抽样活动恰好与当年强烈的拉尼娜现象（La Niña）事件相吻合，这一点在美国国家海洋和大气管理局（NOAA）发布的数据中有所体现。这导致了沿海沙漠通常干燥程度的进一步加剧，以及秘鲁lomas山脉南坡上雾气活动的增加。

排序（Ordinations）是降维技术，允许从一个（有噪声的）数据集中提取主要梯度，对我们来说，即沿着南侧山坡发展的植物群落梯度（详见下一节）。在本章中，我们将对第一个排序轴进行建模，即植物群落梯度，作为海拔、坡度、流域面积和归一化植被指数（NDVI）等环境预测因子的函数。为此，我们将使用随机森林模型——一种非常流行的机器学习算法[@breiman_random_2001]。该模型将使我们能够在研究区域内的任何地方制作植物群落组成的空间分布图。为了保证最佳预测，建议提前使用空间交叉验证。

数据和数据准备

所有后续分析所需的数据都可以通过 spDataLarge 包获得。

1
2
3
data("study_area", "random_points", "comm", package = "spDataLarge")
dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))

study_area 是一个多边形，表示研究区域的轮廓，而 random_points 是一个包含100个随机选取的地点的 sf 对象。comm 是一个采用宽数据格式的群落矩阵[@wickham_tidy_2014]，其中行代表实地考察过的地点，列代表观察到的物种。^[在统计学中，这也被称为列联表或交叉表。]

1
2
3
4
5
6
7
8
9
10
# sites 35 to 40 and corresponding occurrences of the first five species in the
# community matrix
comm[35:40, 1:5]
#>    Alon_meri Alst_line Alte_hali Alte_porr Anth_eccr
#> 35         0         0         0       0.0     1.000
#> 36         0         0         1       0.0     0.500
#> 37         0         0         0       0.0     0.125
#> 38         0         0         0       0.0     3.000
#> 39         0         0         0       0.0     2.000
#> 40         0         0         0       0.2     0.125

值代表每个地点的物种覆盖面积，记录为物种覆盖面积与地点面积的比例（%；请注意，由于单个植物之间的覆盖面积重叠，一个地点的覆盖面积可能>100%）。comm 的行名与 random_points 的 id 列相对应。dem 是研究区域的数字高程模型（DEM），ndvi 是从Landsat场景的红色和近红外通道计算出的归一化植被指数（NDVI）。可视化数据有助于更熟悉它，如下图所示，其中 dem 被 random_points 和 study_area 叠加绘制。


Study mask (polygon), location of the sampling sites (black points) and DEM in the background.

下一步是计算不仅用于建模和预测性制图但也用于将非度量多维缩放（NMDS）轴与研究区域中的主要梯度，即海拔和湿度，对齐的变量。

具体而言，我们使用 R-GIS 桥从数字高程模型计算流域坡度和流域面积。曲率也可能是有价值的预测因子，在练习部分你可以了解它们如何影响建模结果。

为了计算流域面积和流域坡度，我们可以使用 sagang:sagawetnessindex 函数。^[诚然，知道 sagawetnessindex 计算所需的地形属性的唯一方法是熟悉 SAGA。]
qgis_show_help() 返回特定地理算法的所有函数参数和默认值。在这里，我们仅展示完整输出的部分内容。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
# if not already done, enable the saga next generation plugin
qgisprocess::qgis_enable_plugins("processing_saga_nextgen")
# show help
qgisprocess::qgis_show_help("sagang:sagawetnessindex")
#> Saga wetness index (sagang:sagawetnessindex)
#> ...
#> ----------------
#> Arguments
#> ----------------
#> 
#> DEM: Elevation
#> 	Argument type:	raster
#> 	Acceptable values:
#> 		- Path to a raster layer
#> ...
#> SLOPE_TYPE: Type of Slope
#> 	Argument type:	enum
#> 	Available values:
#> 		- 0: [0] local slope
#> 		- 1: [1] catchment slope
#> ...
#> AREA: Catchment area
#> 	Argument type:	rasterDestination
#> 	Acceptable values:
#> 		- Path for new raster layer
#>... 
#> ----------------
#> Outputs
#> ----------------
#> 
#> AREA: <outputRaster>
#> 	Catchment area
#> SLOPE: <outputRaster>
#> 	Catchment slope
#> ...

接下来，我们可以使用 R 的命名参数来指定所需的参数。记住，我们可以使用磁盘上文件的路径或者存储在 R’s全局环境中的 SpatRaster 来指定输入栅格 DEM。指定 SLOPE_TYPE 为1确保算法将返回流域坡度。生成的栅格将保存为具有.sdat扩展名的临时文件，这是 SAGA 的原生栅格格式。

1
2
3
4
5
6
7
8
# environmental predictors: catchment slope and catchment area
ep = qgisprocess::qgis_run_algorithm(
  alg = "sagang:sagawetnessindex",
  DEM = dem,
  SLOPE_TYPE = 1, 
  SLOPE = tempfile(fileext = ".sdat"),
  AREA = tempfile(fileext = ".sdat"),
  .quiet = TRUE)

这将返回一个名为ep的列表，其中包含计算出的输出栅格的路径。让我们将流域面积和流域坡度读入一个多层SpatRaster对象中。另外，我们还将向其中添加两个更多的栅格对象，即dem和ndvi。

1
2
3
4
5
6
7
# read in catchment area and catchment slope
ep = ep[c("AREA", "SLOPE")] |>
  unlist() |>
  terra::rast()
names(ep) = c("carea", "cslope") # assign proper names 
terra::origin(ep) = terra::origin(dem) # make sure rasters have the same origin
ep = c(dem, ndvi, ep) # add dem and ndvi to the multilayer SpatRaster object

此外，集水区域的值高度向右倾斜(‘hist (ep $carea)’)。Log10转换使分布更加正态。

1
ep$carea = log10(ep$carea)

为了方便读者，我们在spDataLarge中添加了ep:

1
ep = terra::rast(system.file("raster/ep.tif", package = "spDataLarge"))

最后，我们可以将地形属性提取到野外观测数据中。

1
2
3
4
# terra::extract adds automatically a for our purposes unnecessary ID column
ep_rp = terra::extract(ep, random_points) |>
  select(-ID)
random_points = cbind(random_points, ep_rp)

数据降维

排序法（Ordinations）在植被科学中是一种流行的工具，主要用于从大多数元素为0的大型物种——样地矩阵中提取主要信息，这些信息通常对应于生态梯度。然而，排序法也广泛应用于遥感、土壤科学、地理营销和许多其他领域。如果你对排序技术不熟悉或需要复习，可以查看Michael W. Palmer的网页，以获得生态学中常用排序技术的简要介绍，或参考@borcard_numerical_2011以深入了解如何在R中应用这些技术。vegan包的文档（通过运行vignette(package = "vegan")查看）也是一个非常有用的资源。

主成分分析（PCA）可能是最著名的排序技术。它是一个优秀的降维工具，尤其适用于变量之间有线性关系的情况，以及两个样地（观察）中变量共同缺失可以被视为相似性的情况。然而，这在植被数据中几乎不是问题。

首先，植物的出现通常沿着一个梯度（如湿度、温度或盐度）呈单峰分布，即非线性关系，最适宜的条件处达到峰值，然后在不适宜的条件下逐渐下降。

其次，两个样地中一个物种的共同缺失很难被视为相似性的指标。假设一个植物物种在我们样本中最干燥（例如极端沙漠）和最潮湿（例如树木草原）的地点都不存在。那么我们确实应该避免把这一点算作是相似性，因为在植被组成方面，这两个完全不同的环境设置很可能唯一的共同点就是物种的共同缺失（除了稀有的无处不在的物种）。

非度量多维尺度（NMDS）是生态学中常用的一种降维技术[@vonwehrden_pluralism_2009]。NMDS通过减少原始矩阵中对象间距离和排序后对象间距离的基于秩的差异来实现降维。这种差异用“应力值”来表示。应力值越低，说明降维效果越好，也就是说，这种低维表示更准确地反映了原始矩阵。应力值低于10表示拟合优秀，约15的值仍然不错，大于20的值则表示拟合差[@mccune_analysis_2002]。

在R中，vegan包的metaMDS()函数可以执行NMDS。作为输入，它需要一个以样地为行、以物种为列的社群矩阵。使用存在-不存在数据进行排序通常会得到更好的结果（就解释方差而言），尽管代价当然是输入矩阵信息较少。decostand() 函数用1和0分别表示物种的存在和缺失，将数值观测转化为存在和缺失数据。像NMDS这样的排序技术至少需要每个样地有一个观察。因此，我们需要去除所有没有发现物种的样地。

1
2
3
4
# presence-absence matrix
pa = vegan::decostand(comm, "pa")  # 100 rows (sites), 69 columns (species)
# keep only sites in which at least one species was found
pa = pa[rowSums(pa) != 0, ]  # 84 rows, 69 columns

所生成的矩阵将作为NMDS（非度量多维标度法）的输入。在这个例子中，我们设定k为4，代表输出轴的数量。^[选择最佳的k值的一种方法是尝试在1到6之间运行NMDS，并选择产生最佳应力值的结果[@mccune_analysis_2002]。]

NMDS是一个迭代过程，每一步都试图让排序后的空间更接近输入矩阵。为确保算法收敛，我们设定迭代次数为500（通过try参数）。

1
2
3
4
5
6
7
8
9
10
11
12
set.seed(25072018)
nmds = vegan::metaMDS(comm = pa, k = 4, try = 500)
nmds$stress
#> ...
#> Run 498 stress 0.08834745 
#> ... Procrustes: rmse 0.004100446  max resid 0.03041186 
#> Run 499 stress 0.08874805 
#> ... Procrustes: rmse 0.01822361  max resid 0.08054538 
#> Run 500 stress 0.08863627 
#> ... Procrustes: rmse 0.01421176  max resid 0.04985418 
#> *** Solution reached
#> 0.08831395

应力值为9代表了非常好的结果，意味着降维后的排序空间能代表输入矩阵的大部分方差。总体来说，NMDS会将更相似（在物种组成方面）的对象在排序空间中放得更近。然而，与大多数其他排序技术不同，这些轴是任意的，并不一定按重要性排序[@borcard_numerical_2011]。然而，我们已经知道湿度是研究区域内主要的梯度[@muenchow_predictive_2013;@muenchow_rqgis:_2017]。由于湿度与海拔高度高度相关，我们将NMDS轴与海拔进行旋转（有关旋转NMDS轴的更多细节，请参见?MDSrotate）。

1
2
3
4
5
6
7
8
9
elev = dplyr::filter(random_points, id %in% rownames(pa)) |> 
  dplyr::pull(dem)
# rotating NMDS in accordance with altitude (proxy for humidity)
rotnmds = vegan::MDSrotate(nmds, elev)
# extracting the first two axes
sc = vegan::scores(rotnmds, choices = 1:2, display = "sites")
# plotting the first axis against altitude
plot(y = sc[, 1], x = elev, xlab = "elevation in m", 
     ylab = "First NMDS axis", cex.lab = 0.8, cex.axis = 0.8)


Plotting the first NMDS axis against altitude.

绘制结果显示，第一个轴正如我们预期的，明显与海拔相关。

第一NMDS轴的得分代表了出现在Mongón山坡上的不同植被形态，即，沿坡度出现的植物群落梯度。为了空间可视化它们，我们可以使用先前创建的预测因子对NMDS得分进行建模，并使用生成的模型进行预测性映射（见下一节）。

植物群落梯度建模

为了在空间上预测植物群落梯度，我们使用随机森林模型[@hengl_random_2018]。随机森林模型在环境和生态建模中经常被应用，并且通常在预测性能方面提供最好的结果[@schratz_hyperparameter_2019]。这里，我们简短地介绍决策树和bagging，因为它们是随机森林的基础。更详细的随机森林及相关技术描述，读者可以参考@james_introduction_2013。

首先，为了通过示例介绍决策树，我们通过将旋转后的NMDS\得分与现场观察数据（random_points）合并，构建一个响应-预测因子矩阵。后续，我们还会用这个生成的数据框进行mlr3建模。

1
2
3
4
5
# construct response-predictor matrix
# id- and response variable
rp = data.frame(id = as.numeric(rownames(sc)), sc = sc[, 1])
# join the predictors (dem, ndvi and terrain attributes)
rp = inner_join(random_points, rp, by = "id")

决策树将预测因子空间划分成若干区域。为了说明这一点，我们使用第一NMDS轴\index{NMDS}的得分（sc）作为响应变量，海拔（dem）作为唯一的预测因子，对我们的数据应用一个决策树。

1
2
3
tree_mo = tree::tree(sc ~ dem, data = rp)
plot(tree_mo)
text(tree_mo, pretty = 0)


Simple example of a decision tree with three internal nodes and four terminal nodes.

生成的树由三个内部节点和四个终端节点组成。树顶部的第一个内部节点将所有观察值中海拔低于328.5米的分到左侧分支，其余的观察值分到右侧分支。落入左侧分支的观察值具有平均NMDS得分为-1.198。

总体而言，我们可以这样解释这棵树：海拔越高，NMDS得分也越高。这意味着这个简单的决策树已经揭示了四个不同的植物群落组合。

决策树有过拟合的倾向，也就是说它们过于贴近输入数据，包括其噪声，从而导致预测性能差[@james_introduction_2013]。Bootstrap集成（bagging）是一种可以帮助克服这个问题的集成技术。集成技术简单地将多个模型的预测结果组合起来。因此，bagging从相同的输入数据中重复抽样，并对预测进行平均。这减少了方差和过拟合，从而大大提高了与决策树相比的预测准确性。

最后，随机森林通过去相关化树来扩展和改进bagging，这是可取的，因为高度相关的树的预测平均值具有更高的方差，因此可靠性较低[@james_introduction_2013]。为了实现这一目标，随机森林使用了bagging，但与传统的bagging不同，其中每棵树都可以使用所有可用的预测因子，随机森林仅使用所有可用预测因子的一个随机样本。

mlr3 构建块

本节中的代码大体上遵循我们在第机器学习超参数调整节中介绍的步骤，唯一的不同之处如下：

  1. 响应变量是数值型的，因此将替换第机器学习超参数调整节中的分类任务为回归任务。
  2. 我们将使用均方根误差（RMSE）作为性能度量，而不是仅适用于分类响应变量的AUROC。
  3. 我们使用随机森林模型而不是支持向量机，这自然涉及到不同的超参数。
  4. 我们将评估偏差减少的性能度量留作读者的练习（参见练习题）。
    相反，我们展示如何为（空间）预测调整超参数。

请记住，当使用100次重复的5折空间交叉验证和50次迭代的随机搜索时，需要125,500个模型来获取偏差减少的性能估计。在超参数调优层面，我们找到了最佳的超参数组合，然后用于外层性能层对特定空间分区的测试数据进行预测。这样做总共产生了500个最佳超参数组合。那么，我们应该使用哪一个来进行空间分布图？答案很简单：完全不用。请记住，调优是为了获取偏差减少的性能估计，而不是为了做最佳可能的空间预测。对于后者，从完整的数据集中估计最佳超参数组合。这意味着，不再需要内部超参数调优层，这是非常合理的，因为我们将模型应用于新的数据（未访问的实地观察）, 这些数据的真实结果是不可用的，因此在任何情况下测试都是不可能的。因此，我们通过一次重复的5折空间CV在完整数据集上调整超参数，以实现良好的空间预测。

既然已经构建了输入变量（rp），我们就准备好指定mlr3的构建块（任务、学习器和重采样）了。为了指定一个空间任务，我们再次使用mlr3spatiotempcv包[@schratz_mlr3spatiotempcv_2021 & Section]，而由于我们的响应（sc）是数值型的，我们使用回归任务。

1
2
3
# create task
task = mlr3spatiotempcv::as_task_regr_st(select(rp, -id, -spri),
  id = "mongon", target = "sc")

使用sf对象作为后端会自动提供之后进行空间分区所需的几何信息。此外，我们去掉了id和spri这两列，因为这些变量不应该用作模型的预测因子。接下来，我们将使用ranger包[@wright_ranger_2017]来构建一个随机森林学习器。

1
lrn_rf = lrn("regr.ranger", predict_type = "response")

与支持向量机不同，随机森林通常在使用其超参数的默认值时就已经表现出良好的性能（这也许是它们受欢迎的一个原因）。然而，调优通常会适度地改善模型结果，因此值得付出努力 [@probst_hyperparameters_2018]。在随机森林中，超参数mtry、min.node.size 和 sample.fraction 决定了随机性的程度，并应进行调优[@probst_hyperparameters_2018]。mtry 指示每棵树应使用多少预测变量。如果使用所有预测变量，那么这实际上就是装袋法。sample.fraction 参数指定应在每棵树中使用的观察样本的比例。较小的比例会导致更大的多样性，从而产生相互间较少相关的树，这通常是可取的（如上文所述）。min.node.size 参数表示一个终端节点至少应具有多少观察值。自然地，随着树和计算时间的增加，min.node.size 会变得更小。

超参数组合将随机选择，但应落在特定的调优限制内（使用 paradox::ps() 创建）。mtry 应在 1 和预测变量的数量（4） 之间变化，sample.fraction 应在 0.2 和 0.9 之间变化，min.node.size 应在 1 和 10 之间变化 [@probst_hyperparameters_2018]。

1
2
3
4
5
6
# specifying the search space
search_space = paradox::ps(
  mtry = paradox::p_int(lower = 1, upper = ncol(task$data()) - 1),
  sample.fraction = paradox::p_dbl(lower = 0.2, upper = 0.9),
  min.node.size = paradox::p_int(lower = 1, upper = 10)
)

确定了搜索空间后，我们可以通过AutoTuner()函数来指定我们的调优。由于我们处理的是地理数据，因此我们将再次利用空间交叉验证来调整超参数。具体来说，我们将使用五折空间划分，且仅进行一次重复（rsmp()）。在这些空间划分中，我们运行50个模型（trm()），同时在预定义的限制（search_space）内使用随机选择的超参数配置（tnr()）来找到最优的超参数\index{hyperparameter}组合[参见机器学习超参数调整章节和 https://mlr3book.mlr-org.com/optimization.html#autotuner, @becker_mlr3_2022]。性能度量是均方根误差（RMSE）。

1
2
3
4
5
6
7
8
autotuner_rf = mlr3tuning::AutoTuner$new(
  learner = lrn_rf,
  resampling = mlr3::rsmp("spcv_coords", folds = 5), # spatial partitioning
  measure = mlr3::msr("regr.rmse"), # performance measure
  terminator = mlr3tuning::trm("evals", n_evals = 50), # specify 50 iterations
  search_space = search_space, # predefined hyperparameter search space
  tuner = mlr3tuning::tnr("random_search") # specify random search
)

调用AutoTuner对象的train()方法最终将执行超参数调优，并会为指定的参数找到最优的超参数组合。

1
2
3
# hyperparameter tuning
set.seed(0412022)
autotuner_rf$train(task)
1
2
3
autotuner_rf$tuning_result
#>    mtry sample.fraction min.node.size learner_param_vals  x_domain regr.rmse
#> 1:    4             0.9             7          <list[4]> <list[3]>     0.375

预测地图

经过调优的超参数现在可以用于预测了。为此，我们只需运行已拟合的AutoTuner对象的predict方法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# predicting using the best hyperparameter combination
autotuner_rf$predict(task)
#> Warning: Detected version mismatch: Learner 'regr.ranger.tuned' has been
#> trained with mlr3 version '0.13.3', not matching currently installed version
#> '0.16.1'
#> Warning: Detected version mismatch: Learner 'regr.ranger' has been trained with
#> mlr3 version '0.13.3', not matching currently installed version '0.16.1'
#> <PredictionRegr> for 84 observations:
#>     row_ids  truth response
#>           1 -1.084   -1.073
#>           2 -0.975   -1.050
#>           3 -0.912   -1.012
#> ---                        
#>          82  0.814    0.646
#>          83  0.814    0.790
#>          84  0.808    0.845

该predict方法将把模型应用于建模时用到的所有观测值。如果有一个多层的SpatRaster，其中包含与建模时使用的预测变量同名的栅格，terra::predict()也将生成空间分布图，即，对新数据进行预测。

1
pred = terra::predict(ep, model = autotuner_rf, fun = predict)


Predictive mapping of the floristic gradient clearly revealing distinct vegetation belts.

如果情况是terra::predict()不支持某个模型算法，你仍然可以手动进行预测。

1
2
3
4
5
6
7
8
9
10
11
newdata = as.data.frame(as.matrix(ep))
colSums(is.na(newdata))  # 0 NAs
# but assuming there were 0s results in a more generic approach
ind = rowSums(is.na(newdata)) == 0
tmp = autotuner_rf$predict_newdata(newdata = newdata[ind, ], task = task)
newdata[ind, "pred"] = data.table::as.data.table(tmp)[["response"]]
pred_2 = ep$dem
# now fill the raster with the predicted values
pred_2[] = newdata$pred
# check if terra and our manual prediction is the same
all(values(pred - pred_2) == 0)

预测性绘图清晰地揭示了不同的植被带。请参考@muenchow_soil_2013以获取lomas山脉上植被带的详细描述。蓝色调表示所谓的Tillandsia-带。Tillandsia是一个高度适应性的属，特别是在lomas山脉的多沙和相当干旱的山脚下以高数量存在。黄色调指的是与Tillandsia-带相比植物覆盖率要高得多的草本植被带。橙色代表了含有最多物种丰富度和植物覆盖的凤梨带。它位于温度逆转下方（大约在海拔750-850米处），那里因雾气而湿度最高。在温度逆转以上，水的可用性自然会减少，景观再次变成沙漠，只有少数多汁植物（多汁带；红色）。

结论

在本章中，我们用NMDS的帮助对lomas蒙贡山的群落矩阵进行了排序。第一个轴，代表研究区域的主要植被梯度，被建模为环境预测变量的函数，其中部分是通过R-GIS桥接得到的。mlr3包提供了空间调整超参数mtry、sample.fraction和min.node.size的构建块。经过调优的超参数作为最终模型的输入，反过来又应用于环境预测变量，以空间地表示植被梯度。结果在沙漠中心展示了令人惊奇的生物多样性。由于lomas山脉正面临严重的威胁，预测图可作为划定保护区和提醒当地居民其所在地区独特性的依据。

在方法论方面，还有一些额外的点可以考虑：

  • 对第二个排序轴也进行建模，并随后找到一种创新方式，将两个轴的模型得分同时在一个预测图中可视化，可能会很有趣
  • 如果我们有兴趣以生态有意义的方式解释模型，我们可能应该使用（半）参数模型[@muenchow_predictive_2013;@zuur_mixed_2009;@zuur_beginners_2017]
    然而，至少有一些方法有助于解释如随机森林这样的机器学习模型（参见，例如，https://mlr-org.github.io/interpretable-machine-learning-iml-and-mlr/）
  • 使用顺序模型基础优化（SMBO）可能优于本章中用于超参数\index{hyperparameter}优化的随机搜索[@probst_hyperparameters_2018]。

最后，请注意，随机森林和其他机器学习模型通常用于观察量多、预测因子多的情境，远超本章所使用的，并且其中哪些变量以及变量之间的交互作用有助于解释响应是不明确的。此外，这些关系可能是高度非线性的。在我们的用例中，响应与预测因子之间的关系相对明确，只有少量的非线性，观测量和预测因子的数量也较低。因此，尝试线性模型可能是值得的。线性模型比随机森林模型更容易解释和理解，因此更应优先选择（简约原则），并且它在计算上也更不费力（参见练习）。如果线性模型无法应对数据中存在的非线性程度，也可以尝试使用广义可加模型（GAM）。关键在于，数据科学家的工具箱包含多种工具，选择最适合手头任务或目的的工具是你的责任。在这里，我们希望向读者介绍随机森林建模以及如何用相应的结果进行预测性地图绘制。为此，一个具有已知响应与预测因子之间关系的研究得当的数据集是合适的。然而，这并不意味着随机森林模型在预测性能方面返回了最佳结果（参见练习）。

练习

E1. 使用社区矩阵的百分比数据运行一个NMDS\index{NMDS}。报告压力值，并将其与使用存在-不存在数据从NMDS中检索到的压力值进行比较。可能是什么原因导致了观察到的差异？

E2. 计算我们在章节中使用过的所有预测栅格\index{raster}（集水区坡度、集水区面积），并将它们放入一个SpatRaster-对象中。添加dem和ndvi。接下来，计算轮廓和切向曲率，并将它们作为额外的预测栅格添加（提示：grass7:r.slope.aspect）。最后，构造一个响应-预测矩阵。第一个NMDS\index{NMDS}轴的分数（当使用存在-不存在社区矩阵时的结果）按照海拔旋转，代表响应变量，并应连接到random_points（使用内连接）。要完成响应-预测矩阵，请将环境预测栅格对象的值提取到random_points中。

E3. 使用空间交叉验证\index{cross-validation!spatial CV}检索随机森林\index{random forest}和线性模型的偏差减少的RMSE。随机森林建模应包括在内部调优循环中估算最优超参数\index{hyperparameter}组合（随机搜索50次迭代）。在调优级别上实现并行化\index{parallelization}。报告平均RMSE\index{RMSE}，并使用箱线图可视化所有检索到的RMSE。请注意，最好使用mlr3函数benchmark_grid()和benchmark()来解决此练习（有关更多信息，请参见https://mlr3book.mlr-org.com/perf-eval-cmp.html#benchmarking）。

-------------已经到底啦-------------