R语言机器学习——BRT提升回归树

R语言BRT介绍与模型解释

模型解释

集成（Ensemble）是合并多个机器学习模型来构建更强大模型的方法。在机器学习文献中有许多模型都属于这一类，但已证明有两种集成模型对大量分类和回归的数据集都是有效的，二者都以决策树为基础，分别是随机森林（Random Forest, RF）和提升树 （Boosted Tree, BT）。梯度提升树（Gradient Boosting Tree, GBT）或称梯度提升决策树（Gradient Boosted Decision Tree, GBDT）是一种优化的提升回归树，在损失函数拟合方面，提升回归树损失函数拟合用的是平方损失，而梯度提升树则是使用损失函数的负梯度来拟合本轮损失的近似值，进而拟合一个回归树。梯度提升树在工业上用途广泛，属于最流行、最实用的算法之一。
提升回归树（Boosted Regression Tree, BRT）特指梯度提升树，这是最主流的方法，在一些R函数中GBM就是指BRT。
BRT理解与区别*

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

在gbm包中，采用的是决策树作为基学习器，重要的参数设置如下：
损失函数的形式(distribution)
迭代次数(n.trees)
学习速率(shrinkage)
再抽样比率(bag.fraction)
决策树的深度(interaction.depth)

损失函数的形式容易设定，分类问题一般选择bernoulli分布，而回归问题可以选择gaussian分布。学习速率方面，我们都知道步子迈得太大容易扯着，所以学习速率是越小越好，但是步子太小的话，步数就得增加，也就是训练的迭代次数需要加大才能使模型达到最优，这样训练所需时间和计算资源也相应加大了。gbm作者的经验法则是设置shrinkage参数在0.01-0.001之间，而n.trees参数在3000-10000之间。
BRT

2. BRT在R语言中的应用

1. gbm：gbm的原始R实现
2. xgboost：一个快速高效的梯度提升框架（C++后端）。
3.h2o：一个强大的基于 Java 的接口，提供并行分布式算法和高效的生产化。
4. dismo：一个用于物种分布建模的包，拓展了gbm包的基础功能。
5. caret：堪称机器学习的百宝箱，可以使用gbm方法进行建模。

2.1 GBM

gbm包
gbm实例

基本构建

gbm(formula = formula(data),
distribution = "bernoulli", #这里需要改成gaussian
data = list(), #我们的数据
weights, #不用管
var.monotone = NULL, #不用管
n.trees = 100, #一般来说先越大越好，然后选择合适的数目
interaction.depth = 1,
n.minobsinnode = 10,
shrinkage = 0.001, #0.001可能太慢了，可以先试试0.01
bag.fraction = 0.5,#从训练样本选择建立下一个树的比例，不用管
train.fraction = 1.0, #从训练样本选择样本数量建立模型的比例，不用管
cv.folds=0, #交叉验证的则数，可以用来提取最适的回归树数目
keep.data = TRUE, #不用管
verbose = FALSE, #如果为TRUE，将输出进度和性能指标
class.stratify.cv=NULL, #不用管
n.cores = NULL) #使用CPU核心的数量

重要参数理解
**（1）distribution：**模型计算损失函数时，需要对输出变量的数据分布做出假设。一般来说，对于分类问题，选择bernoulli或者adaboost，前者更为推荐；对于连续因变量，选择gaussian或者laplace。此外，gbm包还为一些具体问题提供了不少其它选择。

**（2）shrinkage：**学习速率，即每一步迭代中向梯度下降方向前进的速率。一般来说学习速率越小，模型表现越好。令shrinkage=0.001得出的模型几乎一定比shrinkage=0.01的模型好，然而代价是前者运算所需的时间和所耗内存将是后者的10倍。所以选择的准则是在可以承受的时间和内存范围内，shrinkage越小越好。

（3）n.trees：即number of iteration—迭代次数。迭代次数的选择与学习速率密切相关，下图展示了模型表现、学习速率和迭代次数之间的关系
迭代次数可以设得稍微大一点，因为模型训练完后，gbm中的gbm.perf可以估计出最佳迭代次数以供预测阶段使用。在模型训练阶段，gbm作者的经验法则是：3000-10000之间的迭代次数搭配0.01-0.001之间的学习速率。

（4）interaction.depth和n.minobsinnode：子决策树即基础学习器的深度和决策树叶节点包含的最小观测树，若基础学习器训练得过于复杂，将提升模型对于样本的拟合能力而导致过拟合问题，因此子决策树深度不宜过大，叶节点可包含的最小观测书不宜过小。

2.3 h2o.gbm

包介绍

2.4 dismo

gbm.step参数包

gbm.step(data, gbm.x, gbm.y, offset = NULL, fold.vector = NULL, tree.complexity = 1,
learning.rate = 0.01, bag.fraction = 0.75, site.weights = rep(1, nrow(data)),
var.monotone = rep(0, length(gbm.x)), n.folds = 10, prev.stratify = TRUE,
family = "bernoulli", n.trees = 50, step.size = n.trees, max.trees = 10000,
tolerance.method = "auto", tolerance = 0.001, plot.main = TRUE, plot.folds = FALSE,
verbose = TRUE, silent = FALSE, keep.fold.models = FALSE, keep.fold.vector = FALSE,
keep.fold.fit = FALSE, ...)

参数说明：
data : 输入数据框
gbm.x : 数据中预测变量的索引或名称
gbm.y : 数据中响应变量的索引或名称
offset : 抵消
fold.vector : 要读入以进行偏移交叉验证的折叠向量
tree.complexity : 设置单个树的复杂性
learning.rate : 设置应用于子树的权重，权重越小，说明需要树的数量就会越多，模型也越准确
bag.fraction : 设置用于选择变量的观察值的比例，默认为0.75
site.weights : 允许站点的不同权重
var.monotone : 将对单个预测值的响应限制为单调
n.folds : 折叠次数
prev.stratify : 患病率分层褶皱-仅用于存在/不存在数据
family : 族-伯努利（=二项式）、泊松、拉普拉斯或高斯
对于分类问题，选择bernoulli或者adaboost，前者更为推荐；对于连续因变量，选择gaussian或者laplace
n.trees : 要拟合的初始树数
step.size : 每个循环中要添加的树数
max.trees : 停止前要适应的最大树数
实例参考

#加载包
library(dismo)
#导入数据
data(Anguilla_train) #建模数据集
data(Anguilla_test)  #测试数据集
LST_dismo <- gbm.step(data = Anguilla_train, 
   gbm.x = 3:13, gbm.y = 2, family = "gaussian", 
   tree.complexity = 5, learning.rate = 0.005, bag.fraction = 0.5)
LST_dismo$gbm.call$best.trees  # 树的最佳数量950
#树数量越多，意味着模型预测能力相对越强。如果想要更多的树数量，则要将 gbm_step() 中的 learning.rate 改到更小：
summary(LST_dismo)#查看各变量的影响程度

#使用pROC包计算AUC值
library(pROC)
auc(test0$Spring_LST,preds)#Area under the 
curve: 

贡献因子

contributions<-LST_dismo$contributions
#查看各变量的影响程度$rel.inf
a<-summary(LST_dismo)
#然后计算不同类型的因子的贡献比例，画图
library(ggplot2)
ggplot(contributions)+
  geom_col(aes(x=rel.inf,y=var,fill=var))+
  xlab("Ralative Influence")+
  ylab("")+
  guides(fill=F)
#或者直接summary(angaus.tc5.lr01)

barplot(a$rel.inf, 
        names.arg = a$var,  # x轴标签
        main = "Relative influence",   # 主标题
        xlab = "factor",         # x轴标签
        ylab = "importance",            # y轴标签
        col = "skyblue",            # 柱子颜色
        ylim = c(0, 70)             # y轴范围
)

ggplot(a, aes(x = a$rel.inf, y = a$var)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  labs(title = "Relative influence",
       x = "importance",  # 将x轴标签设为Values
       y = "factor") +  # 将y轴标签设为Categories
  coord_flip()+ # 使用coord_flip()函数进行轴互换
  theme_classic()# 去除背景

3. 评价

1.求R方

使用模型预测值和实际观测值之间的平方和来计算总变异（Total Sum of Squares，TSS）和残差平方和（Residual Sum of Squares，RSS）。然后，R方可以通过以下公式计算：
R方 = 1 - (RSS / TSS)
其中，TSS表示总变异，RSS表示残差平方和。你可以使用模型的预测值和实际观测值来计算这些值。
#你还可以考虑使用caret包中的R2函数来计算R方值。R2函数可以计算R方值，它通常适用于各种回归模型，包括Gradient Boosting

library(caret)
r_squared <- R2(pred = predict(LST_dismo, newdata = test_data), obs = test_data$your_response_variable)

实例

library(ggplot2)
library(lattice)
library(caret)
install.packages("caret")
r_squared <- R2(pred=preds, obs = test0$Spring_LST)#0.3559662

2.RMSE

rmse <- sqrt(mean((preds-test0$Spring_LST)^2))#1.456581