一、原理

GBDT + LR 是什么,用在哪

二、说明

GBDT + LR 的结构

RF + LR ? Xgb + LR?

GBDT + LR 模型提升

三、实践

1 如何获得样本落在哪个叶子节点

2 举例

2.2.1 训练集准备

2.2.2 RF+LR

2.2.3 GBDT+LR

2.2.4 Xgboost+LR

2.2.5 单独使用RF, GBDT和Xgboost

2.2.6 结果对比

三、为什么Xgboost+LR的融合效果没有想象中那么好

参考资料

一、原理

为什么要使用LR模型进行融合呢？这是因为LR (逻辑回归) 算法简单有效，成为工业界最常用的算法之一。但 LR 算法是线性模型，不能捕捉到非线性信息，需要大量特征工程找到特征组合。为了发现有效的特征组合，Facebook 在 2014年的论文ractical Lessons from Predicting Clicks on Ads at Facebook介绍了通过 GBDT （Gradient Boost Decision Tree）+ LR 的方案（XGBoost 是 GBDT 的后续发展）。在这篇论文中他们提出了一种将Xgboost作为feature transform的方法。随后在多个Kaggle 竞赛实践中，均证明了此思路的有效性。

GBDT构造组合特征的方式

利用GBDT进行特征构造依据其模型组合方式一共有两种方式：

GBDT + LR

与原论文的实现方式一样，利用GBDT构造组合特征，再将组合特征进行one-hot编码（本实践代码也属此类）；
GBDT + FFM 或者 GBDT + 树模型

此时，使用利用GBDT构造的组合特征不再进行one-hot编码，而是直接利用输出叶节点的索引信息，如果将GBDT组合特征输出到其他树模型，则可直接利用节点索引信息；若是将GBDT信息输出到FFM中，依旧是利用索引信息，但是需要将索引信息组织成FFM数据输入形式。

大概的思想可以描述为如下：先用已有特征训练Xgboost模型，然后利用Xgboost模型学习到的树来构造新特征，最后把这些新特征加入原有特征一起训练模型。构造的新特征向量是取值0/1的，向量的每个元素对应于Xgboost模型中树的叶子结点。当一个样本点通过某棵树最终落在这棵树的一个叶子结点上，那么在新特征向量中这个叶子结点对应的元素值为1，而这棵树的其他叶子结点对应的元素值为0。新特征向量的长度等于XGBoost模型里所有树包含的叶子结点数之和。最后将新的特征扔到LR模型进行训练。

GBDT + LR 是什么,用在哪

本质上GBDT+LR是一种具有stacking思想的二分类器模型，所以可以用来解决二分类问题。这个方法出自于Facebook 2014年的论文 Practical Lessons from Predicting Clicks on Ads at Facebook 。

GBDT+LR 使用最广泛的场景是CTR点击率预估，即预测当给用户推送的广告会不会被用户点击。

点击率预估模型涉及的训练样本一般是上亿级别，样本量大，模型常采用速度较快的LR。但LR是线性模型，学习能力有限，此时特征工程尤其重要。现有的特征工程实验，主要集中在寻找到有区分度的特征、特征组合，折腾一圈未必会带来效果提升。GBDT算法的特点正好可以用来发掘有区分度的特征、特征组合，减少特征工程中人力成本。

二、说明

GBDT + LR 的结构

1）为什么要使用集成的决策树模型，而不是单棵的决策树模型：一棵树的表达能力很弱，不足以表达多个有区分性的特征组合，多棵树的表达能力更强一些。可以更好的发现有效的特征和特征组合
2）为什么建树采用GBDT而非RF：RF也是多棵树，但从效果上有实践证明不如GBDT。且GBDT前面的树，特征分裂主要体现对多数样本有区分度的特征；后面的树，主要体现的是经过前N颗树，残差仍然较大的少数样本。优先选用在整体上有区分度的特征，再选用针对少数样本有区分度的特征，思路更加合理，这应该也是用GBDT的原因。

GBDT+LR 由两部分组成，其中GBDT用来对训练集提取特征作为新的训练输入数据，LR作为新训练输入数据的分类器。

具体来讲，有以下几个步骤：

1 GBDT首先对原始训练数据做训练，得到一个二分类器，当然这里也需要利用网格搜索寻找最佳参数组合。

2 与通常做法不同的是，当GBDT训练好做预测的时候，输出的并不是最终的二分类概率值，而是要把模型中的每棵树计算得到的预测概率值所属的叶子结点位置记为1，这样，就构造出了新的训练数据。

举个例子，下图是一个GBDT+LR 模型结构，设GBDT有两个弱分类器，分别以蓝色和红色部分表示，其中蓝色弱分类器的叶子结点个数为3，红色弱分类器的叶子结点个数为2，并且蓝色弱分类器中对0-1 的预测结果落到了第二个叶子结点上，红色弱分类器中对0-1 的预测结果也落到了第二个叶子结点上。那么我们就记蓝色弱分类器的预测结果为[0 1 0]，红色弱分类器的预测结果为[0 1]，综合起来看，GBDT的输出为这些弱分类器的组合[0 1 0 0 1] ，或者一个稀疏向量（数组）。

这里的思想与One-hot独热编码类似，事实上，在用GBDT构造新的训练数据时，采用的也正是One-hot方法。并且由于每一弱分类器有且只有一个叶子节点输出预测结果，所以在一个具有n个弱分类器、共计m个叶子结点的GBDT中，每一条训练数据都会被转换为1*m维稀疏向量，且有n个元素为1，其余m-n 个元素全为0。

3.3 新的训练数据构造完成后，下一步就要与原始的训练数据中的label(输出)数据一并输入到Logistic Regression分类器中进行最终分类器的训练。思考一下，在对原始数据进行GBDT提取为新的数据这一操作之后，数据不仅变得稀疏，而且由于弱分类器个数，叶子结点个数的影响，可能会导致新的训练数据特征维度过大的问题，因此，在Logistic Regression这一层中，可使用正则化来减少过拟合的风险，在Facebook的论文中采用的是L1正则化。

RF + LR ? Xgb + LR?

有心的同学应该会思考一个问题，既然GBDT可以做新训练样本的构造，那么其它基于树的模型，例如Random Forest以及Xgboost等是并不是也可以按类似的方式来构造新的训练样本呢？没错，所有这些基于树的模型都可以和Logistic Regression分类器组合。至于效果孰优孰劣，我个人觉得效果都还可以，但是之间没有可比性，因为超参数的不同会对模型评估产生较大的影响。下图是RF+LR、GBT+LR、Xgb、LR、Xgb+LR 模型效果对比图，然而这只能做个参考，因为模型超参数的值的选择这一前提条件都各不相同。

顺便来讲，RF也是多棵树，但从效果上有实践证明不如GBDT。且GBDT前面的树，特征分裂主要体现对多数样本有区分度的特征；后面的树，主要体现的是经过前N颗树，残差仍然较大的少数样本。优先选用在整体上有区分度的特征，再选用针对少数样本有区分度的特征，思路更加合理，这应该也是用GBDT的原因。

GBDT + LR 模型提升

现在，我们思考这样一个问题，Logistic Regression是一个线性分类器，也就是说会忽略掉特征与特征之间的关联信息，那么是否可以采用构建新的交叉特征这一特征组合方式从而提高模型的效果？

其次，我们已经在2.3小节中了解到GBDT很有可能构造出的新训练数据是高维的稀疏矩阵，而Logistic Regression使用高维稀疏矩阵进行训练，会直接导致计算量过大，特征权值更新缓慢的问题。

针对上面可能出现的问题，可以翻看我之前的文章：FM算法解析及Python实现，使用FM算法代替LR，这样就解决了Logistic Regression的模型表达效果及高维稀疏矩阵的训练开销较大的问题。然而，这样就意味着可以高枕无忧了吗？当然不是，因为采用FM对本来已经是高维稀疏矩阵做完特征交叉后，新的特征维度会更加多，并且由于元素非0即1，新的特征数据可能也会更加稀疏，那么怎么办？

所以，我们需要再次回到GBDT构造新训练数据这里。当GBDT构造完新的训练样本后，我们要做的是对每一个特征做与输出之间的特征重要度评估并筛选出重要程度较高的部分特征，这样，GBDT构造的高维的稀疏矩阵就会减少一部分特征，也就是说得到的稀疏矩阵不再那么高维了。之后，对这些筛选后得到的重要度较高的特征再做FM算法构造交叉项，进而引入非线性特征，继而完成最终分类器的训练数据的构造及模型的训练。

三、实践

1 如何获得样本落在哪个叶子节点

在实践中的关键点是如何获得每个样本落在训练后的每棵树的哪个叶子结点上。

A、对于xgb/lgbm来说，因为其有sklearn接口和自带接口，因此有两种方法可以获得：

①、sklearn接口。可以设置pre_leaf=True获得每个样本在每颗树上的leaf_Index。（lgbm使用）

②、自带接口。利用apply()方法可以获得leaf indices。（xgb使用）

此过程需注意：无论是设置pre_leaf=True还是利用apply()方法，获得的都是叶子节点的 index，也就是说落在了具体哪颗树的哪个叶子节点上，并非是0/1变量，因此需要自己动手去做 onehot 编码。onehot 可以在 sklearn 的预处理包中调用即可。

B、对于其它的树模型，如随机森林和GBDT，我们只能使用apply()方法获得leaf indices。

2 举例

接下来，我们举例来说明如何利用树模型，尤其是Xgboost来构建新特征，并且是如何与LR模型进行融合。

    params = {"tree_method": "hist", "n_jobs": 4, "n_estimators":50}xgb = XGBClassifier(**params)xgb.fit(X_train, y_train, eval_set=[(X_test, y_test)], eval_metric='auc')print("With XGB Pred")y_train_pred_proba = xgb.predict_proba(X_train)print_eval_metric(y_train, y_train_pred_proba[:, 1])y_test_pred_proba = xgb.predict_proba(X_test)print_eval_metric(y_test, y_test_pred_proba[:, 1])print("---------------")from sklearn.linear_model import LogisticRegressiononehot = OneHotEncoder()X_train_leaves = onehot.fit_transform(xgb.apply(X_train))xgb_lr = LogisticRegression()xgb_lr.fit(X_train_leaves, y_train)print("With XGB + LR Pred")X_train_leaves = onehot.transform(xgb.apply(X_train))y_train_pred_xgb_lr = xgb_lr.predict_proba((X_train_leaves))[:, 1]print_eval_metric(y_train, y_train_pred_xgb_lr)X_test_leaves = onehot.transform(xgb.apply(X_test))y_test_pred_xgb_org_lr = xgb_lr.predict_proba(X_test_leaves)[:, 1]print_eval_metric(y_test, y_test_pred_xgb_org_lr)

XGB (train/test)	XGB+ LR (train/test)

观察leaves的形状

In[181]:  y_pred_leaves.shape
Out[181]: (8000, 100)

共有8000个样本，100棵树（在上面的参数中 num_trees=100),观察第 1 个样本y_pred_train[0]的前10个值：

In[182]:  y_pred_leaves[0][:10]
Out[182]: array([31, 29, 29, 32, 38, 46, 35, 36, 36, 42])

其中第一个数 31 表示这个样本落到了第一颗树的 31 叶子节点，29 表示落到了第二棵树的 29 叶子节点，注意31 、29表示节点编号，从0开始到63

LGBM 例子代码

    params = {"n_jobs": 4, "n_estimators": 50,  "max_depth": 3}lgbm = LGBMClassifier(**params)lgbm.fit(X_train, y_train, eval_set=[(X_test, y_test)], eval_metric='auc')y_test_pred_proba = lgbm.predict_proba(X_test)print_eval_metric(y_test, y_test_pred_proba[:, 1])print("---------------")from sklearn.linear_model import LogisticRegressiononehot = OneHotEncoder()X_train_leaves = onehot.fit_transform(lgbm.predict(X_train, pred_leaf=True))lgbm_lr = LogisticRegression()lgbm_lr.fit(X_train_leaves, y_train)X_test_leaves = onehot.transform(lgbm.predict(X_test, pred_leaf=True))y_test_pred_lgbm_org_lr = lgbm_lr.predict_proba(X_test_leaves)[:, 1]

2.2.1 训练集准备

从代码中，可以看到，我们对训练集X_train又进行了一次切分，生成了训练集X_train, X_train_lr和测试集y_train, y_train_lr。特别要注意数据集的大小，后续我们会进行分析。

注意：我们设定了n_estimator = 10，这意味着树模型中只有10颗树。

import numpy as np
np.random.seed(10)import matplotlib.pyplot as plt
import xgboost as xgbfrom sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import ( RandomForestClassifier,GradientBoostingClassifier)
from sklearn.preprocessing import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_curve,  roc_auc_score
from scipy.sparse import hstackX, y = make_classification(n_samples=80000)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5)# It is important to train the ensemble of trees on a different subset
# of the training data than the linear regression model to avoid
# overfitting, in particular if the total number of leaves is
# similar to the number of training samples
X_train, X_train_lr, y_train, y_train_lr = train_test_split(X_train, y_train, test_size=0.5)n_estimator = 10'''
X_train为20000*20
X_train_lr为20000*20
y_train为20000*1
y_train_lr为20000*1
y_test 为40000*1
'''

2.2.2 RF+LR

我们首先使用随机森林进行实验。这里我们需要对代码进行一个解读。

    # Supervised transformation based on random forestsrf = RandomForestClassifier(max_depth=3, n_estimators=n_estimator)rf_enc = OneHotEncoder()rf_lm = LogisticRegression(solver='lbfgs', max_iter=1000)# 使用随机森林进行训练rf.fit(X_train, y_train)'''rf.apply(X_train)的大小为20000*10，10就是10颗树。rf.apply(X_train)的元素代表了哪一个样本落在了哪一个树的第几个叶子节点上，如32，就代表了落在了一棵树的第32个叶子节点上。rf_enc.fit(rf.apply(X_train))对rf.apply(X_train)进行了onehot编码'''rf_enc.fit(rf.apply(X_train))'''rf_enc.transform(rf.apply(X_train_lr))按照rf.apply(X_train)的编码方式对rf.apply(X_train_lr)进行了onehot编码。这里需要注意得是：我们并没有像常规方式那样对rf.apply(X_train)和rf.apply(X_train_lr)先进行合并再进行onehot，这是因为在训练完成之后随机森林的树模型已经固定，即叶子节点的架构已经确定。因此rf.apply(X_train_lr)的值的范围和rf.apply(X_train)的值的范围必然一样！rf_lm.fit(rf_enc.transform(rf.apply(X_train_lr)), y_train_lr)使用随机森林构造的特征来训练LR'''rf_lm.fit(rf_enc.transform(rf.apply(X_train_lr)), y_train_lr)y_pred_rf_lm = rf_lm.predict_proba(rf_enc.transform(rf.apply(X_test)))[:, 1]fpr_rf_lm, tpr_rf_lm, _ = roc_curve(y_test, y_pred_rf_lm)print("RF+LR的AUC为：", roc_auc_score(y_test, y_pred_rf_lm))

2.2.3 GBDT+LR

我们其次使用随机森林进行实验。其代码与随机森林几乎一样。不再进行代码解读。

# Supervised transformation based on gradient boosted trees
grd = GradientBoostingClassifier(n_estimators=n_estimator)
grd_enc = OneHotEncoder()
grd_lm = LogisticRegression(solver='lbfgs', max_iter=1000)grd.fit(X_train, y_train)
grd_enc.fit(grd.apply(X_train)[:, :, 0])
grd_lm.fit(grd_enc.transform(grd.apply(X_train_lr)[:, :, 0]), y_train_lr)y_pred_grd_lm = grd_lm.predict_proba(grd_enc.transform(grd.apply(X_test)[:, :, 0]))[:, 1]
fpr_grd_lm, tpr_grd_lm, _ = roc_curve(y_test, y_pred_grd_lm)
print("GBT+LR的AUC为：", roc_auc_score(y_test, y_pred_grd_lm))

2.2.4 Xgboost+LR

我们使用Xgboost进行实验。

# Supervised transformation based on xgboost
xgb = xgb.XGBClassifier(nthread=4,     #含义：nthread=-1时，使用全部CPU进行并行运算（默认）, nthread=1时，使用1个CPU进行运算。learning_rate=0.08,    #含义：学习率，控制每次迭代更新权重时的步长，默认0.3。调参：值越小，训练越慢。典型值为0.01-0.2。n_estimators=50,       #含义：总共迭代的次数，即决策树的个数max_depth=5,           #含义：树的深度，默认值为6，典型值3-10。调参：值越大，越容易过拟合；值越小，越容易欠拟合gamma=0,               #含义：惩罚项系数，指定节点分裂所需的最小损失函数下降值。subsample=0.9,       #含义：训练每棵树时，使用的数据占全部训练集的比例。默认值为1，典型值为0.5-1。调参：防止overfitting。colsample_bytree=0.5) #训练每棵树时，使用的特征占全部特征的比例。默认值为1，典型值为0.5-1。调参：防止overfitting。xgb_enc = OneHotEncoder()
xgb_lm = LogisticRegression(solver='lbfgs', max_iter=1000)xgb.fit(X_train, y_train)
xgb_enc.fit(xgb.apply(X_train))
xgb_lm.fit(xgb_enc.transform(xgb.apply(X_train_lr)), y_train_lr)y_pred_xgb_lm = xgb_lm.predict_proba(xgb_enc.transform(xgb.apply(X_test)))[:, 1]
fpr_xgb_lm, tpr_xgb_lm, _ = roc_curve(y_test, y_pred_xgb_lm)
print("xgboost+LR的AUC为：", roc_auc_score(y_test, y_pred_xgb_lm))

在之前的代码中，我们只是用树模型构造的新特征来训练LR。

接下来，我们更近一步，将新特征与原始的20个特征进行拼接形成新的数据集来训练LR。

X_train_ext = hstack([xgb_enc.transform(xgb.apply(X_train_lr)), X_train_lr])
X_test_ext = hstack([xgb_enc.transform(xgb.apply(X_test)), X_test])
xgb_lm.fit(X_train_ext, y_train_lr)y_pred_xgb_originalfeature_lm = xgb_lm.predict_proba(X_test_ext)[:, 1]
fpr_xgb_originalfeature_lm, tpr_xgb_originalfeature_lm, _ = roc_curve(y_test, y_pred_xgb_originalfeature_lm)
print("xgboost新特征与原始特征+LR的AUC为：", roc_auc_score(y_test, y_pred_xgb_originalfeature_lm))

2.2.5 单独使用RF, GBDT和Xgboost

为了进行对比，我们也输出单独使用RF,GBDT和Xgboost的结果。

# The gradient boosted model by itself
y_pred_grd = grd.predict_proba(X_test)[:, 1]
fpr_grd, tpr_grd, _ = roc_curve(y_test, y_pred_grd)
print("GBT的AUC为：", roc_auc_score(y_test, y_pred_grd))# The random forest model by itself
y_pred_rf = rf.predict_proba(X_test)[:, 1]
fpr_rf, tpr_rf, _ = roc_curve(y_test, y_pred_rf)
print("RF的AUC为：", roc_auc_score(y_test, y_pred_rf))# The xgboost model by itself
xgb.fit(X_train, y_train)
y_pred_xgb = xgb.predict_proba(X_test)[:, 1]
fpr_xgb, tpr_xgb, _ = roc_curve(y_test, y_pred_xgb)
print('xgboost的AUC为：' , roc_auc_score(y_test, y_pred_xgb))

2.2.6 结果对比

我们运行整个代码，结果为：

   RF+LR的AUC为： 0.972532755993GBT+LR的AUC为： 0.984711442675xgboost+LR的AUC为： 0.992587688381xgboost新特征与原始特征+LR的AUC为： 0.992632312284GBT的AUC为： 0.98220013158RF的AUC为： 0.965762807823xgboost的AUC为： 0.99284427301

对于RF和GBT，与LR进行融合后的结果要比单独使用RF和GBT要好。

而对于Xgboost，单独使用Xgboost效果最好，其次是xgboost新特征与原始特征+LR，最后才是xgboost+LR。

这与我们预期不符。为什么会出现这样的结果，值得我们讨论。

画图来进一步看下ROC曲线：

plt.figure(1)
plt.plot([0, 1], [0, 1], 'k--')
plt.plot(fpr_rf, tpr_rf, label='RF')
plt.plot(fpr_rf_lm, tpr_rf_lm, label='RF + LR')
plt.plot(fpr_grd, tpr_grd, label='GBT')
plt.plot(fpr_grd_lm, tpr_grd_lm, label='GBT + LR')
plt.plot(fpr_xgb, tpr_xgb, label='XGB')
plt.plot(fpr_xgb_lm, tpr_xgb_lm, label='XGB + LR')
plt.plot(fpr_xgb_originalfeature_lm, tpr_xgb_originalfeature_lm, label='XGB + ori_fea+ LR')
plt.xlabel('False positive rate')
plt.ylabel('True positive rate')
plt.title('ROC curve')
plt.legend(loc='best')
plt.show()plt.figure(2)
plt.xlim(0, 0.2)
plt.ylim(0.8, 1)
plt.plot([0, 1], [0, 1], 'k--')
plt.plot(fpr_rf, tpr_rf, label='RF')
plt.plot(fpr_rf_lm, tpr_rf_lm, label='RF + LR')
plt.plot(fpr_grd, tpr_grd, label='GBT')
plt.plot(fpr_grd_lm, tpr_grd_lm, label='GBT + LR')
plt.plot(fpr_xgb, tpr_xgb, label='XGB')
plt.plot(fpr_xgb_lm, tpr_xgb_lm, label='XGB + LR')
plt.plot(fpr_xgb_originalfeature_lm, tpr_xgb_originalfeature_lm, label='XGB + ori_fea + LR')
plt.xlabel('False positive rate')
plt.ylabel('True positive rate')
plt.title('ROC curve (zoomed in at top left)')
plt.legend(loc='best')
plt.show()