一分钟读完全文
GBM是集成树模型的一种，具有高精度、鲁棒性强以及一定的可解释性。本文介绍了GBM模型的使用全过程。包括调参、训练及最终的feature importance 和 Partial dependence 的绘制，以此说明了如何对集成树模型进行解释。本文基于sklearn实现。相关文章已发表: Hu, Songhua, Kun Xie*, Xiaonian Shan, Hangfei Lin, and Xiaohong Chen. “Modeling Usage Frequencies and Vehicle Preferences in a Large-scale Electric Vehicle Sharing System.” IEEE Intelligent Transportation Systems Magazine 。欢迎引用。
写在前面
所有机器学习的模型在面对统计模型时，最不得不认怂的就是模型的可解释性。做回归时，统计模型有系数有P-value，机器学习却只能谈精度——但当然，若论精度，统计模型折腾来折腾去，还是玩不过机器学习的。
——The tension between interpretability and accuracy is a fundamental trade-off when using models for prediction。
但人们对机器学习模型的可解析性并没有放弃，并一直在努力进行改进。树模型便是其中的佼佼者。集成树中，最出名的当属Random Forest（RF）和Gradient boosting trees（GBM），后者也是近年来大火的XGB的根基。而解释集成树模型的两大利器：Feature importance和Partial dependence，则成了树模型的炫耀资本——高精度+快速+可解析性，是不是觉得很完美。
下面将对GBM进行一个简单的操作介绍。
准备工作
准备工作的第一步当然是掉一大堆包，本文主要借助sklearn下的相关包来完成建模：
import numpy as np
import pandas as pd
import math
import glob
import matplotlib.pyplot as plt
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_squared_error, make_scorer
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.ensemble.partial_dependence import plot_partial_dependence
from sklearn.ensemble.partial_dependence import partial_dependence
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import FormatStrFormatter然后就是把数据导入到DF里面，并对数据进行特征与标签、训练集和测试集的拆分：
os.chdir('/Users/husonghua/Desktop/data')
rawdata = pd.read_csv('rawdata_gbm_new.csv', index_col=0)
rawdata.info()
rawdata = rawdata.dropna()
rawxx = rawdd.drop('Return number', axis=1)
rawyy = rawdd['Return number']
# 拆分数据
X_train, X_test, y_train, y_test = train_test_split(rawxx, rawyy, test_size=0.2, random_state=123)
names = X_train.columns调参
老规矩，借助GRIDSEARCH，比较耗时间，我一般都是把实验室所有电脑打开一起跑，有钱还是买个服务器慢慢训练：
# 调参
cv_params = {'learning_rate': [0.1, 0.05,0.01], 'max_depth': [1,3,5,7,10], 'n_estimators': [100,200,300]}
ind_params = {'random_state': 10}
optimized_GBM = GridSearchCV(GradientBoostingRegressor(**ind_params),
cv_params,
scoring='neg_mean_squared_error', cv=5, n_jobs=-1, verbose=10)
optimized_GBM.fit(X_train, y_train)介绍几个参数：
1）对于GBM，一般训练最多的三个参数就是：
- learning_rate：学习速率，又叫shrinkage，是对每一棵新加上去的树设置一个收缩权重，用于防止过拟合。——Lower learning rate is generally recommended as it makes the model more robust and prevents over-fitting; however, it would require more iterations to get accurate outcome and thus will be computationally expensive.
- max_depth：最大树深度，又叫tree complexity，用于表征单棵树的复杂程度。——A tree with higher depth would be able to capture more specific relations; however, over-fitting would happen since high-depth tree may focus on the relationship of a particular sample.
- n_estimators：迭代次数，这个比较好理解，一般和learning_rate搭配使用：Lowering the learning rate and adding more trees can improve the generalization ability of the model and avoid over-fitting.
2）在GridSearchCV时，主要需要设置的参数有：
- 需要训练的参数写在cv_params中，不需要的写在ind_params中。
- scoring用于指定你想要的评分机制，我这里用的是MSE，注意sklearn内置的MSE是neg_mean_squared_error，也就是说是正常的MSE前加个负数，原因在于对于分数的默认设置是greater_is_better=False。
- cv用于指定交叉验证的折数(cv=5: 5-fold)，不知道交叉验证的请自行补充基础知识。
- n_jobs=-1表示全线程训练，基本n_jobs=-1你电脑就卡得啥也干不了了。
- verbose用于设置交互输出，我把verbose设置为10就会把每一次训练的结果实时输出。
输出最优参数并利用最优参数进行最终模型训练
最终的调参结果都存储在optimized_GBM.grid_scores_，包括每一种参数组合下所对应的5折CV的MSE的均值与方差。可以利用grid_scores_中的数据进一步的绘制调参图：
# 最优组合评分
best_parameters, score, _ = max(optimized_GBM.grid_scores_, key=lambda x: x[1])
print('score:', score)
optimized_GBM.grid_scores_
# 最优参数
for param_name in sorted(best_parameters.keys()):
print("%s: %r" % (param_name, best_parameters[param_name]))
# 训练模型
params = {'n_estimators': 100, 'max_depth': 4, 'learning_rate': 0.1, 'loss': 'huber', 'random_state': 1}
clf = GradientBoostingRegressor(**params, verbose=10)
clf.fit(X_train, y_train)
mse = mean_squared_error(y_test, clf.predict(X_test))
print("MSE: %.4f" % mse)提取重要度并绘制
重要度是树模型最主要的解释工具，顾名思义，重要度越大说明变量越重要。但这个指标也一直被诟病——不稳定性大、重要度量化存在争议（即用什么来表示重要度）、以及无法反应正相关还是负相关。目前有很多人都在对这些工作进行改进，有兴趣可以参考：
https://towardsdatascience.com/interpretable-machine-learning-with-xgboost-9ec80d148d27
[1802.03888] Consistent Individualized Feature Attribution for Tree Ensembles
并且也做出了一些包（比如SHAP）来实现，之后有空我会专门来写写。
feature_importances中存有所有特征的重要度，提取出来绘制一下即可。sklearn本身也有集成好的plot_importances函数。
importance_frame = pd.DataFrame({'Importance': list(clf.feature_importances_), 'Feature': list(names)})
importance_frame.sort_values(by='Importance', inplace=True)
#importance_frame['rela_imp'] = importance_frame['Importance'] / sum(importance_frame['Importance'])
importance_frame.plot(kind='barh', x='Feature', figsize=(8, 8), color='orange')
Plot_partial_dependence
最后便是对Partial_dependence的介绍了。Partial_dependence思路类似于统计中的边际效应，就是在控制其他变量不变的情况下，改变target feature的值，来看模型的fit结果如何变化。具体理论可以参考 Friedman, J.H., 2001. Greedy function approximation: a gradient boosting machine. Annals of statistics, 1189-1232. 顺便说一下，GBM就是这个大佬做的。
Partial_dependence可以说在一定程度上解决了重要度这个指标无法反应正负关系的问题。sklearn中的GBM已经集成了这个功能：Partial Dependence Plots，而早期的XGB是只能绘制重要度的，但近期也可以借助外包如SHAP来实现Partial_dependence的绘制：slundberg/shap
因为变量往往比较多，这里只拎出最主要的几个来画，每画一个都保存一下，axx.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))用于设置Y轴小数点只显示两位。
sorted_idx = np.argsort(feature_importance)
for jj in range(-10, -1):
print(jj)
features = sorted_idx[jj:(jj + 1)]
fig, axx = plot_partial_dependence(clf, X_train, features, feature_names=names, n_jobs=3, grid_resolution=50)
plt.subplots_adjust(left=0.2, top=0.9)
for axx in axx:
for item in ([axx.title, axx.xaxis.label, axx.yaxis.label] + axx.get_xticklabels() + axx.get_yticklabels()):
item.set_fontsize(20)
axx.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
fig.savefig('return month' + str(jj) + '.png')
当然也是可以对双变量进行绘制的，这样看起来会更加的高逼格：
fig = plt.figure()
target_feature = (1, 3)
pdp, axes = partial_dependence(clf, target_feature, X=X_train, grid_resolution=50)
XX, YY = np.meshgrid(axes[0], axes[1])
Z = pdp[0].reshape(list(map(np.size, axes))).T
ax = Axes3D(fig)
surf = ax.plot_surface(XX, YY, Z, rstride=1, cstride=1, cmap=plt.cm.BuPu, edgecolor='k')
ax.set_xlabel(names[target_feature[0]])
ax.set_ylabel(names[target_feature[1]])
ax.set_zlabel('Partial dependence')
# pretty init view
ax.view_init(elev=22, azim=122)
plt.colorbar(surf)
plt.subplots_adjust(top=0.9)
plt.show()
通过Partial_dependence我们可以明显看到相关变量对Y值的影响——relationship is positive or negative。对于写论文的同学，要做一些因果推断的时候，借助Partial_dependence便有理有据了。
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
