增长回归树模型（Boosted Regression Trees）

增长回归树模型（Boosted Regression Trees, BRT）是由 Elith et al. 2008 提出的，其用于生态学统计模型中的解释和预测，对某些典型特征如非线性的变量和变量之间的相互关系有很好的解释和预测。

BRT 是一种拟合统计模型，与传统的拟合吝啬模型有很大不同，其将回归树（regression tree）和增长（boosting）两种方法结合了起来。回归树（regression tree）是一种用递归二进制拆分（recursive binary splits）将因变量与其预测因子联合起来的模型。增长（boosting），则可以结合多个简单模型并提升预测能力。二者结合而成的 BRT 模型，则可以看做是一种加持的回归模型，将一个一个的简单树（simple tree）向前逐步拟合。

BRT 包含基于树模型（tree-based model）的各种优势，能处理不同类型的预测变量，并允许缺失值数据。该方法无需先验数据转换或异常值去除，能够拟合复杂非线性关系，并能自动处理预测因子之间的相互影响。BRT 拟合多个树的过程，能最大程度地弥补单一树模型预测能力弱的缺点，比大多传统模型方法有更强的预测能力，可在模型拟合中处理大量的实际问题。

关于 BRT 的原理和使用方法，在 Elith et al. 2008. A working guide to boosted regression trees 中有详细论述。其使用的数据集是新西兰一种淡水鱼类短鳍鳗 ( Anguilla australis ) 在多个地点上存在或不存在的二进制数据，并包含各个地点的环境因子变量和捕鱼方式等数据。在 R 中使用 dismo 和 gbm 程辑包已可实现 BRT 模型的计算，参见 Elith & Leathwick 2017. Boosted Regression Trees for ecological modeling。该教程实现了模型建立、测试数据集的拟合、空间栅格数据集的拟合等步骤，本文即在 R 中测试该方法。

library(dismo)
library(gbm)

模型建立

data(Anguilla_train)
# 新西兰短鳍鳗在许多地点上的存在(1)与不存在(0)、以及每个地点的环境变量数据
# 环境变量数据是研究区域的栅格化数据
head(Anguilla_train)
#   Site Angaus SegSumT SegTSeas SegLowFlow DSDist DSMaxSlope USAvgT USRainDays USSlope USNative DSDam   Method LocSed
# 1    1      0    16.0    -0.10      1.036  50.20       0.57   0.09      2.470     9.8     0.81     0 electric    4.8
# 2    2      1    18.7     1.51      1.003 132.53       1.15   0.20      1.153     8.3     0.34     0 electric    2.0
# 3    3      0    18.3     0.37      1.001 107.44       0.57   0.49      0.847     0.4     0.00     0      spo    1.0
# 4    4      0    16.7    -3.80      1.000 166.82       1.72   0.90      0.210     0.4     0.22     1 electric    4.0
# 5    5      1    17.2     0.33      1.005   3.95       1.15  -1.20      1.980    21.9     0.96     0 electric    4.7
# 6    6      0    15.1     1.83      1.015  11.17       1.72  -0.20      3.300    25.7     1.00     0 electric    4.5
angaus.tc5.lr01 <- gbm.step(data = Anguilla_train, 
    gbm.x = 3:13, gbm.y = 2, family = "bernoulli", 
    tree.complexity = 5, learning.rate = 0.01, bag.fraction = 0.5)
# tree.complexity 为个体树的复杂程度
# learning.rate 为应用于个体树的权重，权重越小，树的数量会越多
angaus.tc5.lr01$gbm.call$best.trees  # 树的最佳数量
# [1] 900
# 即下图模型中显示的纵向绿色线

查看各变量的影响程度：

summary(angaus.tc5.lr01)
#                   var     rel.inf
# SegSumT       SegSumT 23.94229012
# USNative     USNative 11.76602688
# DSDist         DSDist 11.67527241
# Method         Method 10.37636294
# DSMaxSlope DSMaxSlope  8.54970751
# USSlope       USSlope  8.25540475
# USRainDays USRainDays  8.17452647
# USAvgT         USAvgT  6.63247915
# SegTSeas     SegTSeas  6.56238329
# SegLowFlow SegLowFlow  4.00769505
# DSDam           DSDam  0.05785143

选择树数量

上述模型 angaus.tc5.lr01 的最佳树数量是 900 ，树数量越多，意味着模型预测能力相对越强。如果想要更多的树数量，则要将 gbm_step() 中的 learning.rate 改到更小：

angaus.tc5.lr001 <- gbm.step(data=Anguilla_train, gbm.x = 3:13, gbm.y = 2,
    family = "bernoulli", tree.complexity = 5,
    learning.rate = 0.001, bag.fraction = 0.5)
angaus.tc5.lr001$gbm.call$best.trees
# [1] 5450

模型简化

在变量较多的情况下，如果想提前简化（减少）模型使用的变量，则可以使用 gbm_simplify() 函数查看可被简化的变量的数量：

angaus.simp <- gbm.simplify(angaus.tc5.lr001, n.drops = 5)
# gbm.simplify - version 2.9 
# simplifying gbm.step model for Angaus with 11 predictors and 1000 observations 
# original deviance = 0.6846(0.0267)
# a fixed number of 5 drops will be tested
# creating initial models...
# dropping predictor: 1 2 3 4 5processing final dropping of variables with full data
# 1-DSDam
# 2-SegLowFlow
# 3-SegTSeas
# 4-USAvgT
# 5-DSMaxSlope

# 如下图所示，可以去除的变量有 2 个，即纵向红色线所示
angaus.tc5.lr001.simp <- gbm.step(Anguilla_train, 
    gbm.x=angaus.simp$pred.list[[2]], gbm.y=2,
    tree.complexity=5, learning.rate=0.001)
# 在实际应用中据情况而定

绘制模型的函数及拟合值

# 绘制每个变量的拟合函数的曲线
gbm.plot(angaus.tc5.lr001)

# 绘制 gbm 对象的拟合值
# 可更可靠地反映拟合状况，尤其是当预测变量与相关参数
# 和/或样本在环境空间中分布不均匀的时候
gbm.plot.fits(angaus.tc5.lr001, v = 0, mask.presence = FALSE) 
# 默认 v = 0 为输出全部变量，可选 v 的值以专门输出指定的变量

# 结果中，图上的 “wtm” 表示：相对于每个非因子预测变量的拟合值的加权平均值

变量之间的相互影响

使用函数 gbm.interactions() 查看环境变量之间的相互作用。使用 gbm.perspec() 函数绘制变量之间相互作用 3D 图。

find.int <- gbm.interactions(angaus.tc5.lr001)
find.int$interactions   # 所有变量之间的相互作用，矩阵式的表
find.int$rank.list      # 相互作用最高的 6 组值
#   var1.index var1.names var2.index var2.names int.size
# 1          7 USRainDays          1    SegSumT    21.52
# 2         11     Method          1    SegSumT     8.79
# 3          4     DSDist          1    SegSumT     8.06
# 4          6     USAvgT          1    SegSumT     7.85
# 5          9   USNative          1    SegSumT     4.72
# 6          4     DSDist          2   SegTSeas     4.58
# 最显著的相互作用发生在变量 `USRainDays` (7) 和 `SegSumT` (1) 之间
gbm.perspec(angaus.tc5.lr001, x = 7, y = 1, z.range = c(0, 0.7))

使用模型预测新数据集

在使用建立的模型模拟新数据集时，必须确保新数据集与训练数据集的表结构、列名称等一致。

data(Anguilla_test)
head(Anguilla_test)
#   Angaus_obs SegSumT SegTSeas SegLowFlow DSDist DSMaxSlope USAvgT USRainDays USSlope USNative DSDam Method LocSed
# 1          0    16.6     1.01      1.017   5.23       0.29  -1.40      1.980    10.0     1.00     0   <NA>    4.9
# 2          1    16.8    -0.51      1.002   2.24       0.00   0.27      0.460     0.7     0.00     0   <NA>    2.3
# 3          0    16.3     0.76      1.023 162.28       5.14  -0.60      0.806    21.4     0.66     0   <NA>    4.3
# 4          0    15.6     1.56      1.003   4.05       0.57   1.14      3.300     0.9     0.75     0   <NA>    1.0
# 5          0    14.6    -0.20      1.023 127.03       1.72  -1.90      1.940    28.9     0.97     0   <NA>     NA
# 6          0    16.7     0.80      1.003   2.48      14.57  -1.50      1.980    22.0     0.99     0   <NA>    2.8
preds <- predict.gbm(angaus.tc5.lr001, Anguilla_test, 
    n.trees = angaus.tc5.lr001$gbm.call$best.trees, 
    type = "response")

# 计算观测值与预测值之间的偏差
calc.deviance(obs = Anguilla_test$Angaus_obs, 
    weights = rep(1,length(obs)), family="binomial", 
    pred = preds, calc.mean = TRUE)
# [1] 0.7415488

d <- cbind(Anguilla_test$Angaus_obs, preds)  # 将观测值和预测值放在同一表中
pres <- d[d[,1] == 1, 2]                     # 观测值 = 1，存在
abs <- d[d[,1] == 0, 2]                      # 观测值 = 0，不存在

# 存在与不存在数据之间的交叉验证
e <- evaluate(p = pres, a = abs)
e
# class          : ModelEvaluation 
# n presences    : 107 
# n absences     : 393 
# AUC            : 0.8616918   # 特征曲线下的面积
# cor            : 0.5265351   # 相关系数
# max TPR+TNR at : 0.1129443   # TPR 真实正比率、TNR 真实负比率

使用模型预测空间栅格数据

这里使用的空间栅格数据，应当是一些栅格图层的集合，这些栅格图层的名称和值分别对应于模型中的各种环境变量。本例用模型预测 dismo 程辑包中的 Anguilla_grids 空间栅格数据，即预测栅格中各地点存在短鳍鳗的概率大小。

data(Anguilla_grids)
Anguilla_grids
# class      : RasterBrick 
# dimensions : 196, 250, 49000, 11  (nrow, ncol, ncell, nlayers)
# resolution : 100, 100  (x, y)
# extent     : 0, 25000, 0, 19600  (xmin, xmax, ymin, ymax)
# crs        : NA 
# source     : memory
# names      :    DSDam,   DSDist, DSMaxSlope,   LocSed, SegLowFlow,  SegSumT, SegTSeas,   USAvgT, USNative, USRainDays,  USSlope 
# min values :  0.00000,  0.76455,    0.06000,  1.00000,    1.00000, 16.80000,  0.47000, -2.07000,  0.00000,    0.84710,  0.05000 
# max values :   0.0000, 149.8895,    12.9500,   5.5000,     3.1640,  18.4000,   1.3300,   0.7200,   1.0000,     1.2200,  22.2700

该栅格数据包含了 11 个环境变量栅格组成的栅格 brick ，但是没有 “Method” 的栅格，因此需要用特定值（Method 的一种）定义一个数据框架，应用于预测函数 raster::predict 中的 const = 中：

elec <- data.frame(factor(c("electric")))
colnames(elec) <- "Method"

sp.pred.elec <- predict(Anguilla_grids, angaus.tc5.lr001, const = elec, 
    n.trees = angaus.tc5.lr001$gbm.call$best.trees, type ="response")
sp.pred.elec <- mask(sp.pred.elec, raster(Anguilla_grids, 1))
plot(sp.pred.elec, main = "Angaus - BRT prediction - electric")

此外，如果将训练集中 “Method” 的所有方法都列举出来，则会出现不同的结果：

par(mfrow = c(2, 3))

elec <- data.frame(factor(c("electric")))
colnames(elec) <- "Method"
sp.pred.elec <- predict(Anguilla_grids, angaus.tc5.lr001, const=elec, n.trees=angaus.tc5.lr001$gbm.call$best.trees, type="response")
sp.pred.elec <- mask(sp.pred.elec, raster(Anguilla_grids, 1))
plot(sp.pred.elec, main = "Angaus - BRT prediction - electric")

net <- data.frame(factor(c("net")))
colnames(net) <- "Method"
sp.pred.net <- predict(Anguilla_grids, angaus.tc5.lr001, const=net, n.trees=angaus.tc5.lr001$gbm.call$best.trees, type="response")
sp.pred.net <- mask(sp.pred.net, raster(Anguilla_grids, 1))
plot(sp.pred.net, main = "Angaus - BRT prediction - net")

mix <- data.frame(factor(c("mixture")))
colnames(mix) <- "Method"
sp.pred.mix <- predict(Anguilla_grids, angaus.tc5.lr001, const=mix, n.trees=angaus.tc5.lr001$gbm.call$best.trees, type="response")
sp.pred.mix <- mask(sp.pred.mix, raster(Anguilla_grids, 1))
plot(sp.pred.mix, main = "Angaus - BRT prediction - mixture")

spo <- data.frame(factor(c("spo")))
colnames(spo) <- "Method"
sp.pred.spo <- predict(Anguilla_grids, angaus.tc5.lr001, const= spo, n.trees=angaus.tc5.lr001$gbm.call$best.trees, type="response")
sp.pred.spo <- mask(sp.pred.spo, raster(Anguilla_grids, 1))
plot(sp.pred.spo, main = "Angaus - BRT prediction - spo")

trap <- data.frame(factor(c("trap")))
colnames(trap) <- "Method"
sp.pred.trap <- predict(Anguilla_grids, angaus.tc5.lr001, const=trap, n.trees=angaus.tc5.lr001$gbm.call$best.trees, type="response")
sp.pred.trap <- mask(sp.pred.trap, raster(Anguilla_grids, 1))
plot(sp.pred.trap, main = "Angaus - BRT prediction - trap")


all <- data.frame(factor(c("electric", "mixture", "net", "spo", "trap")))
colnames(all) <- "Method"
sp.pred.all <- predict(Anguilla_grids, angaus.tc5.lr001, const=all, n.trees=angaus.tc5.lr001$gbm.call$best.trees, type="response")
sp.pred.all <- mask(sp.pred.all, raster(Anguilla_grids, 1))
plot(sp.pred.all, main = "Angaus - BRT prediction - all")

至于这几种不同捕鱼方式的预测结果有什么意义，作者并没有提及，而是只输出了 “electric” 方法的预测结果。有待进一步说明。