通过Aggregated boosted tree（ABT）评估解释变量的重要性
几天前一同学咨询了一个问题，如何通过Aggregated boosted tree（ABT）评估变量的相对重要性。周末抽空了解了一下，顺便进行了简单的整理在这里作个分享。
先来看一下文献中通常使用ABT回答哪些问题。就以同学当时提供的参考文献为例，正文中的图4部分，作者通过ABT分析了地上植物以及多种土壤理化属性对土壤真菌群落的影响，该图展示了这些指标的相对重要性，数值越大代表该指标的重要程度越高。
（Wang et al, 2020)
Aggregated boosted tree（ABT）的简单描述
ABT建立在boosted tree的基础上，是boosted tree的延伸类型。
首先来看一下什么是boosted tree。类似随机森林（random forest），boosted tree也是决策树（decision tree）的拓展类型，同样是有监督的机器学习方法，旨在实现准确的预测和解释。决策树、随机森林以及boosted tree三者的特点可简单概括为：
（1）决策树是一种简单的决策生成图，非常容易理解，但单个决策树的预测和解释能力很弱。相比之下，随机森林或boosted tree是大量决策树的融合形式，尽管解释复杂，但结果更加稳健。
（2）随机森林在过程中独立生成大量的决策树，并在结束阶段，基于平均或者多数原则实现树的合并获得最终模型。
（3）boosted tree一次构建一棵树，工作方式类似加性模型，采用前向选择的方式，不断引入一个弱学习者（weak learner）来改善现有弱学习者的不足。树合并过程起始于开始阶段，并在整个过程中逐步组合每一棵树，最终实现所有树的合并获得最终模型。
ABT作为boosted tree的某种延伸类型，因此具有boosted tree的特点：对于非平衡的数据而言，boosted tree通常比随机森林的效果很好；但如果数据中的噪声很多，则boosted tree通常不是一个好的选择，因为它更容易导致过拟合，并且也往往比随机森林更难调整。相比传统的线性相关或多元线性回归等方法，ABT模型更擅长处理变量间的非线性以及相互作用，并同时定量评估各解释变量对响应变量的相对影响。
ABT的具体计算过程，同学们若有兴趣可参考De'ath（2007）的文章。白鱼同学也就粗略地知道个大概情况，但细节公式啥的直接就看晕了……
R包gbmplus执行ABT评估变量的重要性
接下来的内容，展示如何通过R语言执行ABT的方法过程。考虑到大多数文献中使用ABT的重点都是解释变量的效应，很少用作预测模型来使用，因此下文只展示如何通过ABT评估变量的重要性。
ABT中，响应变量可以是连续变量，也可以为类别变量。
当响应变量为连续变量时，ABT执行了回归的功能，返回的结果代表了解释变量对响应变量数值改变程度的相对贡献。
当响应变量为类别变量时，ABT执行了分类的功能，返回的结果代表了解释变量对区分已知分类差异的相对重要性。
下文所使用的示例数据来自De'ath（2007）文章中的补充材料1。
白鱼同学已经将它们下载了下来，并和R代码一起上传到网盘中（提取码，7c7r）：
https://pan.baidu.com/s/1Ausna2z1erkVtJbhrYgHKw
若百度盘失效，也可在GitHub的备份中获取：
https://github.com/lyao222lll/sheng-xin-xiao-bai-yu
关于gbmplus包的安装
白鱼同学查阅了好几篇使用到ABT方法的文献，均提到通过R语言gbmplus包执行ABT。
由于未再查阅更多的文献，所以不清楚还有没有其它的包可以执行ABT。同学们也可以自行查找下，毕竟gbmplus包的体验非常不佳……gbmplus包好多年没更新了，高版本的R不支持，要找个低于3.0版本的R，然后在这个链接中手动下载gbmplus包进行安装：http://www.esapubs.org/archive/ecol/E088/015/suppl-1.htm
白鱼同学安装了一个非常古老的R2.7版本后，gbmplus包确实测试成功了
。好了也不纠结了，继续看下文使用。
响应变量是连续变量时的ABT
示例数据“softcorals.csv”，调查了澳大利亚大堡礁（GBR）的150个珊瑚礁中软珊瑚属的丰富度（Richness）；并记录了5个环境因素，包括沉积厚度（Sediment）、可见度（Visibility）、波浪冲击（Wave）、坡度角（Slope）和周围水流速（Flow）；以及2个空间变量，跨GBR（Across）和沿GBR（Along）的相对距离。
读取数据并加载R包，建立ABT以评估各类环境或空间因素对软珊瑚属丰富度的效应，解释软珊瑚属丰富度变化的成因。
代码语言：javascript复制##响应变量是连续变量时，用于回归
library(gbmplus)
#读取软珊瑚属丰富度及环境和空间因素数据
softcorals <- read.csv('softcorals.csv', check.names = FALSE)
#将软珊瑚属丰富度作为响应变量，环境和空间因素作为自变量，构建 ABT 执行回归
#详情 ?gbm，该示例生成 500 个随机树用作 boosting，并使用 5 折交叉验证估计误差，其余使用默认值
set.seed(123)
fit_softcorals <- gbm(Richness~Across Along Visibility Slope Flow Wave Sediment,
data = softcorals, n.trees = 100, cv.folds = 5)
par(las = 1)
summary(fit_softcorals)
结果中，给出了各类环境或空间因素对软珊瑚属丰富度的效应，Relative influence表示了相对影响的百分比。
可以看到，跨GBR和沿GBR的相对距离（Along 和Across）、波浪冲击（Wave）是解释软珊瑚属丰富度变化的比较重要的因素。其它环境因素的效应则相对微弱。
响应变量是类别变量时的ABT
示例数据“barramundi.csv”，在淡水和河口栖息地收集了澳洲肺鱼的鳞片，共计270个样本中，141个来自淡水（Fresh，记录为1），129个来自河口（Fresh，记录为0）。测量了每个鳞片样本中的元素组成，包括锶（Sr）、钡（Ba）、钙（Ca）、铁（Fe）、钾（K）、镁（Mg）、锰（Mn）、磷（P）和硫（S），测量数据中所有变量均经过对数转换。
期望寻找鳞片中重要的元素，它们的含量能够较好地区分鱼类的来源（淡水或河口）。这里，通过建立ABT模型实现该目的。
其实过程和上文的连续型响应变量的是一样的。
代码语言：javascript复制##响应变量是类别变量时，用于分类
library(gbmplus)
#读取鱼鳞数据及其元素组成含量数据
barramundi <- read.csv('barramundi.csv', check.names = FALSE)
barramundi$Fresh <- as.factor(barramundi$Fresh)
#将鱼鳞样本来源（淡水或河口）作为响应变量，鱼鳞元素组成含量为解释变量，构建 ABT 执行分类
#详情 ?gbm，该示例生成 500 个随机树用作 boosting，并使用 5 折交叉验证估计误差，其余使用默认值
set.seed(123)
fit_barramundi <- gbm(Fresh~Sr Ba Mn Mg Fe Zn K P S,
data = barramundi, n.trees = 100, cv.folds = 5)
par(las = 1)
summary(fit_barramundi)
结果中，给出了鱼鳞中各类化学元素含量对区分鱼类已知来源（淡水或河口）的相对重要性，Relative influence表示了相对影响的百分比。
可以看到，锶（Sr）是最重要的元素；铁（Fe）、钡（Ba）、硫（S）等次之，但效应已经相对微弱了；其余元素的重要性则更低。
* 关于变量重要性的柱形图的美化
如果觉得默认出图不好看，不妨通过ggplot2重新绘制，这里以上文中软珊瑚属丰富度的结果为例展示一个简单的作图例子。
为了能正常运行gbmplus包，上文使用了一个早期版本的R（v2.7），然而若在2.7版本的R中安装新包（如ggplot2）就特别难受……例如，白鱼同学选择将上述结果导出，然后打开平时最常使用的3.6版本的R，加载ggplot2后作图……
#将变量重要性的结果信息提取后输出
softcorals_var_influence <- summary(fit_softcorals)
write.csv(softcorals_var_influence, 'softcorals.var_influence.csv', row.names = FALSE, quote = FALSE)
#然后打开一个已经安装 ggplot2 的较新的 R 版本
#加载 ggplot2，读取数据后重新绘制变量重要性的柱形图
library(ggplot2)
softcorals_var_influence <- read.csv('softcorals.var_influence.csv', stringsAsFactors = FALSE)
softcorals_var_influence <- softcorals_var_influence[order(softcorals_var_influence$rel.influence), ]
softcorals_var_influence$var <- factor(softcorals_var_influence$var, levels = softcorals_var_influence$var)
#颜色代表不同的变量
p <- ggplot(softcorals_var_influence, aes(var, rel.influence))
coord_flip()
theme(panel.grid = element_blank(), panel.background = element_blank(),
axis.line = element_line(colour = 'black'))
scale_y_continuous(expand = c(0, 0))
labs(x = '', y = 'Relative influence(%)', title = '')
p geom_col(aes(fill = var), width = 0.7, show.legend = FALSE)
scale_fill_manual(values = c('#8DD3C7', '#FFFFB3', '#BEBADA',
'#FB8072', '#80B1D3', '#FDB462', '#B3DE69'))
#或者颜色以渐变色表示变量的重要性
p geom_col(aes(fill = rel.influence), width = 0.7, show.legend = FALSE)
scale_fill_gradientn(colours = colorRampPalette(c('#86CEF7', '#0000F6'))(10))
参考资料
https://www.datasciencecentral.com/profiles/blogs/decision-tree-vs-random-forest-vs-boosted-trees-explained
De'ath G. Boosted trees for ecological modeling and prediction. Ecology, 2007, 88(1): 243-251.
Wang J, Liu G, Zhang C, et al. Effect of long-term destocking on soil fungal functional groups and interactions with plants. Plant and Soil, 2020: 1-14.
