Python-数据科学项目第二版-三-

Python 数据科学项目第二版（三）

原文：annas-archive.org/md5/c89da1393a37db56d0a53ed5ccaa3f00

译者：飞龙

协议：CC BY-NC-SA 4.0

第六章：6. 梯度提升、XGBoost 与 SHAP 值

概述

阅读完本章后，您将能够描述梯度提升的概念，这是 XGBoost 包的基本理念。然后，您将通过合成数据训练 XGBoost 模型，并在此过程中学习早期停止以及几个 XGBoost 超参数。除了使用我们之前讨论的相似方法来生成树（通过设置max_depth），您还将发现 XGBoost 提供的另一种生成树的新方式：基于损失的树生成。在学习了 XGBoost 之后，您将接触到一种新的、强大的模型预测解释方法，称为SHAP（SHapley Additive exPlanations）。您将看到如何使用 SHAP 值为任何数据集的模型预测提供个性化的解释，而不仅仅是训练数据，同时也理解 SHAP 值的加法属性。

介绍

正如我们在上一章中看到的，基于决策树和集成模型提供了强大的方法来创建机器学习模型。尽管随机森林已经存在了几十年，但最近关于另一种树集成方法——梯度提升树的研究，已经产生了最先进的模型，这些模型在预测建模领域，尤其是在使用结构化表格数据（如案例研究数据）方面，占据了主导地位。如今，机器学习数据科学家使用的两个主要包是 XGBoost 和 LightGBM，用于创建最准确的预测模型。在本章中，我们将使用合成数据集熟悉 XGBoost，然后在活动中将其应用到案例研究数据中。

注意

也许使用 XGBoost 的最佳动机之一，来自于描述这个机器学习系统的论文，尤其是在 Kaggle 这个流行的在线机器学习竞赛论坛的背景下：

“在 2015 年 Kaggle 博客上发布的 29 个挑战获胜解决方案中，有 17 个解决方案使用了 XGBoost。在这些解决方案中，8 个仅使用 XGBoost 来训练模型，而大多数其他解决方案则将 XGBoost 与神经网络结合使用在集成中。相比之下，第二受欢迎的方法——深度神经网络，在 11 个解决方案中得到了使用。”（陈和郭斯特林，2016，dl.acm.org/doi/abs/10.1145/2939672.2939785）

正如我们将看到的，XGBoost 将我们迄今为止讨论的几种不同思想结合起来，包括决策树、集成建模以及梯度下降。

除了性能更强的模型，近年来的机器学习研究还提供了更详细的模型预测解释方式。我们不再依赖仅仅表示模型训练集整体的解释方法，比如逻辑回归系数或随机森林的特征重要性，而是通过一个名为 SHAP 的新包，能够为任何我们希望的 dataset（如验证集或测试集）逐个解释模型预测。这对于帮助我们数据科学家以及我们的业务合作伙伴深入理解模型的工作原理非常有帮助，即使是对于新数据也是如此。

梯度提升和 XGBoost

什么是提升（Boosting）？

提升（Boosting）是一种创建多个机器学习模型或估计器集成的方法，类似于随机森林模型背后的袋装（bagging）概念。像袋装方法一样，提升可以与任何类型的机器学习模型一起使用，但它通常用于构建决策树的集成模型。与袋装方法的一个主要区别是，在提升过程中，每个新加入集成的估计器都依赖于之前加入的所有估计器。由于提升过程是按顺序进行的，并且集成成员的预测结果会加总以计算整体集成预测结果，因此它也被称为逐步加法建模（stagewise additive modeling）。袋装和提升的区别可以如图 6.1所示：

](tos-cn-i-73owjymdk6/b3f8705ca1b74df2b23c9752c69561a5)

图 6.1：袋装与提升的对比

虽然袋装方法使用不同的训练数据随机样本训练多个估计器，但提升则通过利用上一轮估计器错误分类的样本信息来训练新估计器。通过将新的估计器集中在这些错误分类的样本上，目标是使得整体集成在整个训练数据集上的表现更好。AdaBoost，作为XGBoost的前身，通过为错误分类的样本赋予更多的权重，达到了这一目标，从而训练新的估计器加入集成。

梯度提升与 XGBoost

XGBoost 是一种建模过程和 Python 包，是当今最受欢迎的机器学习方法之一，因其在许多领域（从商业到自然科学）的卓越表现而广受欢迎。XGBoost 还证明是机器学习竞赛中最成功的工具之一。我们不会讨论 XGBoost 实现的所有细节，而是了解它如何工作的高层次概念，并查看一些最重要的超参数。有关详细信息，感兴趣的读者应参考 Tianqi Chen 和 Carlos Guestrin 的论文《XGBoost: A Scalable Tree Boosting System》（dl.acm.org/doi/abs/10.1145/2939672.2939785）。

XGBoost 实现的梯度提升过程是一个逐步加法模型，类似于 AdaBoost。然而，XGBoost 并不是直接在模型训练过程中对分类错误的样本赋予更多权重，而是使用一种类似于梯度下降的过程。回想一下第四章，偏差方差权衡，梯度下降优化利用了损失函数（成本函数的另一种称呼）导数的信息，在训练逻辑回归模型时更新估计的系数。损失函数的导数包含了关于每次迭代中如何调整系数估计的方向和幅度的信息，从而减少预测中的误差。

XGBoost 将梯度下降的思想应用于逐步加法建模，利用损失函数的梯度（导数的另一种说法）信息训练新的决策树，加入到集成模型中。实际上，XGBoost 在梯度下降的基础上更进一步，正如第四章，偏差方差权衡中所描述的，它同时利用了损失函数的一阶和二阶导数。使用误差梯度训练决策树的方法是第五章，决策树与随机森林中引入的节点不纯度思想的替代方案。从概念上讲，XGBoost 训练新树的目标是将集成预测朝着减少误差的方向移动。朝这个方向迈多大的步伐由 learning_rate 超参数控制，类似于第 4.01 章 使用梯度下降最小化成本函数中的 learning_rate，该内容也出现在第四章，偏差方差权衡中。

到此为止，我们应该已经掌握了足够的知识，了解 XGBoost 如何工作，可以开始动手实践了。为了说明 XGBoost，我们将使用 scikit-learn 的 make_classification 函数创建一个用于二分类的合成数据集。该数据集包含 5,000 个样本和 40 个特征。这里的其余选项控制了分类任务的难度，你应当查阅 scikit-learn 文档以更好地理解它们。特别值得注意的是，我们将使用多个簇（n_clusters_per_class），意味着在多维特征空间中会有几个属于某一类别的点的区域，类似于上一章图 5.3中展示的簇。基于树的模型应该能够识别这些簇。此外，我们还指定在 40 个特征中只有 3 个是有信息的（n_informative），另外有 2 个冗余特征（n_redundant），它们包含与有信息特征相同的信息。因此，总的来说，在 40 个特征中，只有 5 个特征在做出预测时是有用的，而且这些信息都编码在其中的 3 个特征中。

如果你想在自己的电脑上跟随本章的示例进行操作，请参考 Jupyter 笔记本：packt.link/L5oS7：

from sklearn.datasets import make_classification
X, y = make_classification(n_samples=5000, n_features=40,\
                           n_informative=3, n_redundant=2,\
                           n_repeated=0, n_classes=2,\
                           n_clusters_per_class=3,\
                           weights=None, flip_y=0.05,\
                           class_sep=0.1, hypercube=True,\
                           shift=0.0,scale=1.0, shuffle=True,\
                           random_state=2)

请注意，响应变量y的类别比例大约为 50%：

y.mean()

这应该输出以下内容：

0.4986

本章中，我们将一次性将此合成数据集划分为训练集和验证集，而不是使用交叉验证。然而，我们在这里介绍的概念可以扩展到交叉验证场景。我们将把合成数据集划分为 80%的训练集和 20%的验证集。在现实世界的数据问题中，我们还希望保留一个测试集，供以后评估最终模型时使用，但在这里我们不做此操作：

from sklearn.model_selection import train_test_split
X_train, X_val, y_train, y_val = \
train_test_split(X, y, test_size=0.2, random_state=24)

现在我们已经为建模准备好了数据，接下来需要实例化一个XGBClassifier类的对象。请注意，我们现在将使用 XGBoost 包，而不是 scikit-learn，来开发预测模型。然而，XGBoost 的 API（应用程序接口）设计得与 scikit-learn 类似，因此使用这个类应该是直观的。XGBClassifier类可以用来创建一个模型对象，并提供fit、predict方法以及其他常见功能，我们还可以在实例化类时指定模型的超参数。我们在此仅指定几个已经讨论过的超参数：n_estimators是用于模型的提升轮数（换句话说，是阶段性加法建模过程中的阶段数），objective是用来计算梯度的损失函数，learning_rate控制每个新估计器对集成模型的贡献，或者本质上控制减少预测误差的步长。剩余的超参数与在模型训练过程中希望看到的输出量（verbosity）以及即将被弃用的label_encoder选项有关，XGBoost 的开发者建议将其设置为False：

xgb_model_1 = xgb.XGBClassifier(n_estimators=1000,\
                                verbosity=1,\
                                use_label_encoder=False,\
                                objective='binary:logistic',\
                                learning_rate=0.3)

我们所指定的超参数值表示：

• 我们将使用 1,000 个估计器，或提升轮次。稍后我们会详细讨论需要多少轮次；默认值是 100。

• 我们熟悉二元逻辑回归的目标函数（也称为代价函数），这在第四章《偏差-方差权衡》中有所介绍。XGBoost 还提供了各种目标函数，适用于分类和回归等多种任务。

• 学习率设置为0.3，这是默认值。可以通过超参数搜索程序探索不同的值，我们将在后续演示。

  注意

  推荐使用 Anaconda 环境安装 XGBoost 和 SHAP，具体方法可参考前言。如果安装与此处所示版本不同，您的结果可能与此处展示的不同。

现在我们有了模型对象和一些训练数据，可以开始拟合模型了。这与在 scikit-learn 中的做法类似：

%%time
xgb_model_1.fit(X_train, y_train,\
                eval_metric="auc",\
                verbose=True)

在这里，我们正在使用 Jupyter notebook 中的 %%time “单元魔法”来跟踪拟合过程所需的时间。我们需要提供训练数据的特征 X_train 和响应变量 y_train。我们还需要提供 eval_metric 并设置详细程度，我们稍后将解释。执行此单元格后，应显示类似以下的输出：

CPU times: user 52.5 s, sys: 986 ms, total: 53.4 s
Wall time: 17.5 s
Out[7]:
XGBClassifier(base_score=0.5, booster='gbtree',\
              colsample_bylevel=1, colsample_bynode=1,\
              colsample_bytree=1, gamma=0, gpu_id=-1,\
              importance_type='gain',interaction_constraints='',\
              learning_rate=0.3, max_delta_step=0, max_depth=6,\
              min_child_weight=1, missing=nan,\
              monotone_constraints='()', n_estimators=1000,\
              n_jobs=4, num_parallel_tree=1, random_state=0,\
              reg_alpha=0, reg_lambda=1, scale_pos_weight=1,\
              subsample=1, tree_method='exact',\
              use_label_encoder=False, validate_parameters=1,\
              verbosity=1)

输出告诉我们，这个单元格执行了 17.5 秒，称为“墙时”，即你墙上的时钟上显示的经过时间。CPU 时间比这个更长，因为 XGBoost 高效地同时使用多个处理器。XGBoost 还打印出所有超参数，包括我们设置的和那些保留默认值的超参数。

现在，为了检查这个拟合模型的性能，我们将评估验证集上的 ROC 曲线下的面积。首先，我们需要获取预测的概率：

val_set_pred_proba = xgb_model_1.predict_proba(X_val)[:,1]
from sklearn.metrics import roc_auc_score
roc_auc_score(y_val, val_set_pred_proba)

此单元格的输出应如下所示：

0.7773798710782294

这表示 ROC AUC 约为 0.78。这将是我们模型性能的基准，使用 XGBoost 的几乎默认选项。

XGBoost 超参数

早期停止

在使用 XGBoost 训练决策树集成时，有许多选项可以用于减少过拟合并利用偏差-方差权衡。早期停止是其中一种简单的方式，可以帮助自动回答“需要多少次提升轮次？”的问题。值得注意的是，早期停止依赖于拥有一个独立的验证集数据，而不仅仅是训练集。然而，这个验证集实际上会在模型训练过程中使用，因此它不算作“未见过”的数据，这与我们在第四章《偏差-方差权衡》中使用验证集来选择模型超参数的方式类似。

当 XGBoost 训练连续的决策树以减少训练集上的误差时，可能会出现添加越来越多的树到集成中，能够提供越来越好的拟合训练数据，但开始导致在保留数据上的性能下降。为避免这种情况，我们可以使用一个验证集，也叫做评估集或 XGBoost 称之为 eval_set。评估集将作为特征及其对应响应变量的元组列表提供。该列表中最后的元组将用于早期停止。我们希望这个验证集，因为训练数据将用于拟合模型，无法提供对样本外泛化的估计：

eval_set = [(X_train, y_train), (X_val, y_val)]

现在我们可以再次拟合模型，但这次我们提供eval_set关键字参数，并传入我们刚刚创建的评估集。此时，eval_metric中的auc变得非常重要。这意味着在每轮提升之后，在训练另一个决策树之前，ROC 曲线下面积将在所有通过eval_set提供的数据集上进行评估。由于我们将verbosity=True，因此我们将在单元格下方看到输出，其中包含训练集和验证集的 ROC AUC。这提供了一个实时视图，展示随着更多提升轮次的训练，模型在训练和验证数据上的性能变化。

由于在预测建模中，我们主要关心模型在新数据和未见过的数据上的表现，当我们发现进一步的提升轮次不再改善模型在样本外数据上的表现时，我们希望停止训练额外的提升轮次。early_stopping_rounds=30参数表示一旦完成 30 轮提升，且在验证集上的 ROC AUC 没有任何进一步的提升，XGBoost 就会停止模型训练。一旦模型训练完成，最终拟合的模型将只包含获得最高验证集性能所需的集成成员。这意味着最后的 30 个集成成员将被丢弃，因为它们没有提供验证集性能的提升。现在我们来拟合这个模型，并观察进展：

%%time
xgb_model_1.fit(X_train, y_train, eval_set=eval_set,\
                eval_metric='auc',\
                verbose=True, early_stopping_rounds=30)

输出应该类似于这样：

[0]	validation_0-auc:0.80412	validation_1-auc:0.75223
[1]	validation_0-auc:0.84422	validation_1-auc:0.79207
[2]	validation_0-auc:0.85920	validation_1-auc:0.79278
[3]	validation_0-auc:0.86616	validation_1-auc:0.79517
[4]	validation_0-auc:0.88261	validation_1-auc:0.79659
[5]	validation_0-auc:0.88605	validation_1-auc:0.80061
[6]	validation_0-auc:0.89226	validation_1-auc:0.80224
[7]	validation_0-auc:0.89826	validation_1-auc:0.80305
[8]	validation_0-auc:0.90559	validation_1-auc:0.80095
[9]	validation_0-auc:0.91954	validation_1-auc:0.79685
[10]	validation_0-auc:0.92113	validation_1-auc:0.79608
…
[33]	validation_0-auc:0.99169	validation_1-auc:0.78323
[34]	validation_0-auc:0.99278	validation_1-auc:0.78261
[35]	validation_0-auc:0.99329	validation_1-auc:0.78139
[36]	validation_0-auc:0.99344	validation_1-auc:0.77994
CPU times: user 2.65 s, sys: 136 ms, total: 2.78 s
Wall time: 2.36 s
…

请注意，这比之前的拟合花费了更少的时间。这是因为通过提前停止，我们只训练了 37 轮提升（请注意提升轮次是零索引的）。这意味着提升过程只需 8 轮就能达到最佳验证得分，而不是我们之前尝试的 1,000 轮！你可以通过模型对象的booster属性访问达到最佳验证集得分所需的提升轮数以及该得分。这个属性提供了比我们一直在使用的 scikit-learn API 更底层的模型接口：

xgb_model_1.get_booster().attributes()

输出应该如下所示，确认迭代次数和最佳验证得分：

{'best_iteration': '7', 'best_score': '0.80305'}

从训练过程，我们还可以看到每一轮的 ROC AUC，包括训练数据的validation_0-auc和验证数据的validation_1-auc，这些都能提供有关随着提升过程进行，模型是否过拟合的洞见。在这里，我们可以看到验证得分在第 8 轮之前一直在上升，但之后开始下降，表明进一步提升可能会导致模型过拟合。然而，训练得分则持续上升，直到过程被终止，这展示了 XGBoost 在拟合训练数据方面的强大能力。

我们还可以进一步确认拟合后的模型对象仅表示七轮提升，并通过手动计算 ROC AUC（如前所述）来检查验证集的表现：

val_set_pred_proba_2 = xgb_model_1.predict_proba(X_val)[:,1]
roc_auc_score(y_val, val_set_pred_proba_2)

这应该输出如下内容：

0.8030501882609966

这与经过七轮提升后获得的最高验证分数一致。因此，通过对模型训练过程进行一次简单的调整——使用验证集和提前停止，我们成功地将模型在验证集上的表现从大约 0.78 提升到 0.80，取得了显著的提高。这展示了提前停止在提升中的重要性。

这里自然会有一个问题：“我们怎么知道 30 轮提前停止就足够了？”你可以像调整任何超参数一样尝试不同的轮数，对于不同的数据集，可能需要不同的值。你可以观察验证分数在每轮提升中的变化来大致判断。有时，验证分数在每轮之间可能会出现跳跃性波动，因此最好有足够的轮数，以确保找到最大值，并且越过任何暂时的下降。

调整学习率

学习率在 XGBoost 文档中也被称为eta，以及步长缩减。这个超参数控制每个新估算器对集成预测的贡献大小。如果你增加学习率，你可能会更快地达到最佳模型，即在验证集上表现最好的模型。然而，设置学习率过高的风险在于，这可能会导致提升步骤过大。在这种情况下，梯度提升过程可能无法收敛到最优模型，这与第四章 偏差方差权衡中的练习 4.01，使用梯度下降最小化成本函数中讨论的梯度下降中的大学习率问题类似。接下来，让我们探讨学习率如何影响我们合成数据上的模型表现。

学习率是一个介于零和 1 之间的数字（包含端点，尽管零学习率没有实际意义）。我们创建一个包含 25 个在 0.01 和 1 之间均匀分布的数字的数组，用来测试不同的学习率：

learning_rates = np.linspace(start=0.01, stop=1, num=25)

现在我们设置一个for循环，为每个学习率训练一个模型，并将验证分数保存在数组中。我们还会追踪达到最佳迭代所需的提升轮数。接下来的几个代码块应该作为一个单元格在 Jupyter notebook 中一起运行。我们首先通过测量所需的时间，创建空列表来存储结果，并开启for循环：

%%time
val_aucs = []
best_iters = []
for learning_rate in learning_rates:

在每次循环迭代中，learning_rate变量将保存learning_rates数组的连续元素。一旦进入循环，第一步是用新的学习率更新模型对象的超参数。这是通过set_params方法完成的，我们用双星号**和一个字典来传递超参数名称和值。Python 中的**函数调用语法允许我们传递任意数量的关键字参数，也称为kwargs，以字典的形式传递。在这种情况下，我们只改变一个关键字参数，所以字典中只有一个项：

    xgb_model_1.set_params(**{'learning_rate':learning_rate})

现在我们已经在模型对象上设置了新的学习率，我们按照之前的方法使用早期停止来训练模型：

    xgb_model_1.fit(X_train, y_train, eval_set=eval_set,\
                    eval_metric='auc',\
                    verbose=False, early_stopping_rounds=30)

经过拟合后，我们获得了验证集的预测概率，并使用这些概率来计算验证集的 ROC AUC。然后，我们使用append方法将其添加到结果列表中：

    val_set_pred_proba_2 = xgb_model_1.predict_proba(X_val)[:,1]
    val_aucs.append(roc_auc_score(y_val, val_set_pred_proba_2))

最后，我们还捕捉了每个学习率所需的轮次数量：

     best_iters.append(
        int(xgb_model_1.get_booster().\
                        attributes()['best_iteration']))

前面提到的五个代码片段应该一起在一个单元格中运行。输出结果应该类似于以下内容：

CPU times: user 1min 23s, sys: 526 ms, total: 1min 24s
Wall time: 22.2 s

现在我们已经得到了这个超参数搜索的结果，可以可视化验证集的性能和迭代次数。由于这两个指标的尺度不同，我们希望创建一个双* y *轴图。pandas 使这变得简单，因此首先我们将所有数据放入一个数据框：

learning_rate_df = \
pd.DataFrame({'Learning rate':learning_rates,\
              'Validation AUC':val_aucs,\
              'Best iteration':best_iters})

现在我们可以像这样可视化不同学习率下的性能和迭代次数，注意：

• 我们设置了索引（set_index），使得学习率绘制在* x 轴上，其他列绘制在 y *轴上。

• secondary_y关键字参数指示要绘制在右侧* y *轴上的列。

• style参数允许我们为每一列指定不同的线条样式。-o是带点的实线，而--o是带点的虚线：

  mpl.rcParams['figure.dpi'] = 400
  learning_rate_df.set_index('Learning rate')\
  .plot(secondary_y='Best iteration', style=['-o', '--o'])
  

结果图应该如下所示：

图 6.2：XGBoost 模型在验证集上的表现，以及不同学习率下直到最佳迭代的提升轮次

总的来说，似乎较小的学习率能够使这个合成数据集上的模型表现更好。通过使用小于默认值 0.3 的学习率，我们能得到的最佳表现如下所示：

max(val_aucs)

输出如下：

0.8115309360232714

通过调整学习率，我们将验证 AUC 从大约 0.80 提高到 0.81，表明使用适当的学习率是有益的。

通常，较小的学习率通常会导致更好的模型性能，尽管它们需要更多的提升轮次，因为每轮的贡献较小。这将导致模型训练所需的时间增加。在图 6.2的轮次与最佳迭代的关系图中，我们可以看到这一点。在这种情况下，看起来可以通过少于 50 轮来获得良好的性能，并且对于这组数据，模型训练时间不会很长。对于较大的数据集，训练时间可能会更长。根据你的计算时间，如果减少学习率并训练更多轮次，可能是提高模型性能的有效方法。

在探索较小学习率时，务必将n_estimators超参数设置得足够大，以便训练过程能够找到最佳模型，理想情况下与早停策略结合使用。

XGBoost 中的其他重要超参数

我们已经看到，XGBoost 中的过拟合可以通过使用不同的学习率以及早停策略来进行补偿。还有哪些其他可能相关的超参数？XGBoost 有很多超参数，我们不会在这里列出所有内容。建议你查阅文档（xgboost.readthedocs.io/en/latest/parameter.html）获取完整列表。

在接下来的练习中，我们将对六个超参数的范围进行网格搜索，其中包括学习率。我们还将包含max_depth，该参数在第五章《决策树与随机森林》中应该已经熟悉，它控制集成中树的生长深度。除了这些，我们还会考虑以下内容：

• gamma通过仅在损失函数值减少超过一定量时才允许节点分裂，从而限制集成中树的复杂性。

• min_child_weight也通过仅在节点的“样本权重”达到一定量时才分裂节点，从而控制树的复杂性。如果所有样本的权重相等（如我们练习中一样），这就相当于节点中的最小训练样本数。这类似于 scikit-learn 中决策树的min_weight_fraction_leaf和min_samples_leaf。

• colsample_bytree是一个随机选择的特征子集，用于在集成中生长每棵树。这类似于 scikit-learn 中的max_features参数（该参数在节点级别进行选择，而不是在树级别进行选择）。XGBoost 还提供了colsample_bylevel和colsample_bynode，分别在每棵树的每一层和每个节点进行特征采样。

• subsample控制从训练数据中随机选择的样本比例，在为集成模型生长新树之前进行选择。这类似于 scikit-learn 中随机森林的bootstrap选项。subsample和colsample参数都限制了模型训练期间可用的信息，增加了各个集成成员的偏差，但希望也能减少整体集成的方差，从而改善样本外模型的表现。

正如你所看到的，XGBoost 中的梯度提升树实现了许多来自决策树和随机森林的概念。现在，让我们探索这些超参数如何影响模型性能。

练习 6.01：用于调优 XGBoost 超参数的随机网格搜索

在这个练习中，我们将使用随机网格搜索来探索六个超参数的空间。当你有许多超参数值需要搜索时，随机网格搜索是一个不错的选择。在这里，我们将研究六个超参数。例如，如果每个超参数有五个值需要测试，我们就需要进行56 = 15,625*次搜索。即使每次模型拟合只需要一秒钟，我们也仍然需要几个小时来穷举所有可能的组合。通过只搜索这些组合的随机样本，随机网格搜索也能取得令人满意的结果。这里，我们将展示如何使用 scikit-learn 和 XGBoost 来实现这一点。

随机网格搜索的第一步是为每个超参数指定你希望从中采样的值的范围。你可以通过提供值列表或分布对象来实现这一点。对于离散超参数，如max_depth，它只有少数几个可能的值，因此可以将其指定为列表。而对于连续超参数，如subsample，它可以在区间(0, 1]上任意变化，因此我们不需要指定一个值列表。相反，我们可以要求网格搜索在这个区间内以均匀的方式随机采样值。我们将使用均匀分布来采样我们考虑的几个超参数：

注意

本练习的 Jupyter 笔记本可以在packt.link/TOXso找到。

1. 从scipy导入uniform分布类，并使用字典指定所有需要搜索的超参数的范围。uniform可以接受两个参数，loc和scale，分别指定采样区间的下界和区间的宽度：

   from scipy.stats import uniform
   param_grid = {'max_depth':[2,3,4,5,6,7],
                 'gamma':uniform(loc=0.0, scale=3),
                 'min_child_weight':list(range(1,151)),
                 'colsample_bytree':uniform(loc=0.1, scale=0.9),
                 'subsample':uniform(loc=0.5, scale=0.5),
                 'learning_rate':uniform(loc=0.01, scale=0.5)}
   

   在这里，我们根据实验和经验选择了参数范围。例如，对于subsample，XGBoost 文档建议选择至少 0.5 的值，因此我们表示为uniform(loc=0.5, scale=0.5)，这意味着从区间[0.5, 1]中采样。

2. 现在我们已经指明了要从哪些分布中采样，我们需要进行采样。scikit-learn 提供了ParameterSampler类，它将随机采样提供的param_grid参数，并返回请求的样本数量（n_iter）。我们还设置了RandomState，以便在不同的 notebook 运行中得到可重复的结果：

   from sklearn.model_selection import ParameterSampler
   rng = np.random.RandomState(0)
   n_iter=1000
   param_list = list(ParameterSampler(param_grid, n_iter=n_iter,
                                      random_state=rng))
   

   我们已经将结果以字典列表的形式返回，字典包含特定的参数值，对应于 6 维超参数空间中的位置。

   请注意，在本次练习中，我们正在遍历 1,000 个超参数组合，这可能需要超过 5 分钟。您可能希望减少这个数字，以便获得更快的结果。

3. 检查param_list的第一个项目：

   param_list[0]
   

   这应该返回一个包含六个参数值的组合，来自所指示的分布：

   {'colsample_bytree': 0.5939321535345923,
    'gamma': 2.1455680991172583,
    'learning_rate': 0.31138168803582195,
    'max_depth': 5,
    'min_child_weight': 104,
    'subsample': 0.7118273996694524}
   

4. 观察如何使用**语法，通过字典同时设置多个 XGBoost 超参数。首先，为本次练习创建一个新的 XGBoost 分类器对象。

   xgb_model_2 = xgb.XGBClassifier(
       n_estimators=1000,
       verbosity=1,
       use_label_encoder=False,
       objective='binary:logistic')
   xgb_model_2.set_params(**param_list[0])
   

   输出应该显示所设置的超参数：

   XGBClassifier(base_score=0.5, booster='gbtree',\
                 colsample_bylevel=1, colsample_bynode=1,\
                 colsample_bytree=0.5939321535345923,\
                 gamma=2.1455680991172583, gpu_id=-1,\
                 importance_type='gain',interaction_constraints='',\
                 learning_rate=0.31138168803582195,\
                 max_delta_step=0, max_depth=5,\
                 min_child_weight=104, missing=nan,\
                 monotone_constraints='()', n_estimators=1000,\
                 n_jobs=4, num_parallel_tree=1,\
                 random_state=0, reg_alpha=0, reg_lambda=1,\
                 scale_pos_weight=1, subsample=0.7118273996694524,\
                 tree_method='exact', use_label_encoder=False,\
                 validate_parameters=1, verbosity=1)
   

   我们将在循环中使用此过程，以查看所有超参数值。

5. 接下来的几个步骤将包含在for循环的一个单元格中。首先，测量执行此操作所需的时间，创建一个空列表以保存验证 AUC 值，然后启动计数器：

   %%time
   val_aucs = []
   counter = 1
   

6. 打开for循环，设置超参数，并拟合 XGBoost 模型，类似于调整学习率的前述示例：

   for params in param_list:
       #Set hyperparameters and fit model
       xgb_model_2.set_params(**params)
       xgb_model_2.fit(X_train, y_train, eval_set=eval_set,\
                       eval_metric='auc',\
                       verbose=False, early_stopping_rounds=30)
   

7. 在for循环内，获取预测概率和验证集 AUC 值：

       #Get predicted probabilities and save validation ROC AUC
       val_set_pred_proba = xgb_model_2.predict_proba(X_val)[:,1]
       val_aucs.append(roc_auc_score(y_val, val_set_pred_proba))
   

8. 由于这个过程需要几分钟时间，最好将进度打印到 Jupyter notebook 输出中。我们使用 Python 余数语法%，每 50 次迭代打印一次信息，换句话说，当counter除以 50 的余数为零时。最后，我们增加计数器：

       #Print progress
       if counter % 50 == 0:
           print('Done with {counter} of {n_iter}'.format(
               counter=counter, n_iter=n_iter))
       counter += 1
   

9. 将步骤 5-8 组合在一个单元格中并运行 for 循环，应该得到如下输出：

   Done with 50 of 1000
   Done with 100 of 1000
   …
   Done with 950 of 1000
   Done with 1000 of 1000
   CPU times: user 24min 20s, sys: 18.9 s, total: 24min 39s
   Wall time: 6min 27s
   

10. 现在我们已经得到所有超参数探索的结果，我们需要对其进行检查。我们可以很容易地将所有超参数组合放入数据框中，因为它们是作为字典列表组织的。做这件事并查看前几行：

    xgb_param_search_df = pd.DataFrame(param_list)
    xgb_param_search_df.head()
    

    输出应该像这样：

    

    图 6.3：来自随机网格搜索的超参数组合

11. 我们还可以将验证集的 ROC AUC 值添加到数据框中，查看最大值：

    xgb_param_search_df['Validation ROC AUC'] = val_aucs
    max_auc = xgb_param_search_df['Validation ROC AUC'].max()
    max_auc
    

    输出应该如下：

    0.8151220995602575
    

    在超参数空间搜索的结果是验证集 AUC 大约为 0.815。虽然这个值比我们通过提前停止和学习率搜索得到的 0.812 稍大（图 6.3），但差距不大。这意味着，对于这组数据，默认的超参数（除了学习率）已经足够实现相当好的性能。虽然我们没有通过超参数搜索大幅提升性能，但观察超参数变化如何影响模型表现依然具有启发性。在接下来的步骤中，我们将逐个检查 AUC 相对于每个参数的边际分布。这意味着我们将查看当一个超参数逐渐变化时，AUC 如何变化，并牢记在网格搜索结果中其他超参数也在变化。

12. 使用以下代码设置一个六个子图的网格，用于绘制每个超参数与性能之间的关系，同时调整图像分辨率并启动一个计数器，我们将用它来遍历子图：

    mpl.rcParams['figure.dpi'] = 400
    fig, axs = plt.subplots(3,2,figsize=(8,6))
    counter = 0
    

13. 打开一个 for 循环来遍历超参数名称，它们是数据框的列，不包括最后一列。通过将 subplot 返回的 3 x 2 数组展平，并使用 counter 索引来访问轴对象。对于每个超参数，使用数据框的 plot.scatter 方法在适当的坐标轴上绘制散点图。x 轴将显示超参数，y 轴显示验证 AUC，其他选项帮助我们获得具有白色内部的黑色圆形标记：

    for col in xgb_param_search_df.columns[:-1]:
        this_ax = axs.flatten()[counter]
        xgb_param_search_df.plot.scatter(x=col,\
                                         y='Validation ROC AUC',\
                                         ax=this_ax, marker='o',\
                                         color='w',\
                                         edgecolor='k',\
                                         linewidth=0.5)
    

14. 数据框的 plot 方法会自动创建 x 和 y 轴标签。然而，由于 y 轴标签对于所有这些图来说是相同的，我们只需要在第一个图中包含它。因此，我们将其他所有图的标签设置为空字符串 ''，并递增计数器：

        if counter > 0:
            this_ax.set_ylabel('')
        counter += 1
    

    由于我们将绘制边际分布图，在观察验证集 AUC 随着给定超参数变化时，其他所有超参数也会发生变化。这意味着关系可能会有噪声。为了了解整体趋势，我们还将创建折线图，显示每个超参数的分位数中验证 AUC 的平均值。分位数将数据按值是否落入最低 10%、接下来的 10% 等，直到最高 10% 来组织成不同的箱体。pandas 提供了一个名为qcut的函数，可以将一个 Series 切分成分位数（分位数是均匀大小的箱体中的一个，例如在 10 个箱体的情况下是一个十分位数），返回另一个 Series，包含这些分位数以及分位数的边界值，您可以把它们看作是直方图的边缘。

15. 使用 pandas 的qcut为每个超参数（除了max_depth）生成一个分位数序列（10 个分位数），返回区间边界（对于 10 个分位数会有 11 个边界），如果唯一值不足以分成 10 个分位数，则删除不需要的区间边界（duplicates='drop'）。创建一个列表，包含每对区间边界之间的中点，用于绘图：

        if col != 'max_depth':
            out, bins = pd.qcut(xgb_param_search_df[col], q=10,\
                                retbins=True, duplicates='drop')
            half_points = [(bins[ix] + bins[ix+1])/2
                           for ix in range(len(bins)-1)]
    

16. 对于max_depth，由于只有六个唯一值，我们可以像处理分位数一样直接使用这些值：

        else:
            out = xgb_param_search_df[col]
            half_points = np.sort(xgb_param_search_df[col].unique())
    

17. 通过复制超参数搜索数据框创建一个临时数据框，创建一个包含分位数序列的新列，并利用此列查找每个超参数分位数内的验证 AUC 平均值：

        tmp_df = xgb_param_search_df.copy()
        tmp_df['param_decile'] = out
        mean_df = tmp_df.groupby('param_decile').agg(
            {'Validation ROC AUC':'mean'})
    

18. 我们可以通过在每个散点图的同一坐标轴上，绘制表示每个分位数平均值的虚线图来可视化结果。关闭for循环并使用plt.tight_layout()清理子图格式：

        this_ax.plot(half_points,\
                     mean_df.values,\
                     color='k',\
                     linestyle='--')
    plt.tight_layout()
    

    运行for循环后，生成的图像应如下所示：

    

    图 6.4：验证集 AUC 与每个超参数的关系图，并显示每个超参数分位数内的平均值

    尽管我们注意到，本次练习中的超参数搜索并未显著提高验证 AUC，相较于本章之前的尝试，但图 6.4中的图表仍能展示 XGBoost 超参数如何影响该特定数据集的模型表现。XGBoost 通过在生成树时限制可用数据来对抗过拟合，一种方式是随机选择每棵树可用的特征的一部分（colsample_bytree），或者随机选择训练样本的一部分（subsample）。然而，在该合成数据中，似乎当每棵树使用 100% 的特征和样本时，模型表现最佳；低于此比例时，模型表现逐渐下降。另一种控制过拟合的方法是通过控制树的复杂度来限制集成中的树，方法包括限制max_depth、叶子节点中的最小训练样本数（min_child_weight）或分裂节点所需的最小损失函数减少值（gamma）。在我们这里的例子中，max_depth和gamma似乎对模型表现没有太大影响，而限制叶子节点中的样本数量似乎会带来负面效果。

    看起来在这个案例中，梯度提升过程本身足够稳健，能够在没有额外技巧的情况下实现良好的模型表现，以减少过拟合。然而，正如我们上面所观察到的，较小的learning_rate是有益的。

19. 我们可以显示最佳的超参数组合及其对应的验证集 AUC，如下所示：

    max_ix = xgb_param_search_df['Validation ROC AUC'] == max_auc
    xgb_param_search_df[max_ix]
    

    这应该返回类似于以下的数据框行：

](tos-cn-i-73owjymdk6/1d97a47eb33f489c8f8705e3a56c8093)

图 6.5：最佳超参数组合与验证集 AUC

验证集 AUC 与我们在上一步（步骤 10）通过仅调整学习率所达到的结果相似。

另一种生长树的方式：XGBoost 的 grow_policy

除了通过max_depth超参数限制树的最大深度外，还有一种控制树生长的范式：寻找一个节点，在这个节点上分裂将导致损失函数的最大减少，并进行分裂，而不考虑这会使树变得多深。这可能导致树有一个或两个非常深的分支，而其他分支可能不会生长得很远。XGBoost 提供了一个名为grow_policy的超参数，设置为lossguide时会产生这种树的生长方式，而depthwise选项是默认设置，会将树生长到指定的max_depth，正如我们在第五章《决策树与随机森林》中所做的，以及在本章至今为止的操作一样。lossguide生长策略是 XGBoost 中的一个较新选项，它模拟了 LightGBM（另一个流行的梯度提升包）的行为。

要使用lossguide策略，需要设置我们尚未讨论的另一个超参数tree_method，它必须设置为hist或gpu-hist。不详细讲解，hist方法将使用更快的方式来搜索分裂。它不会在节点中对训练样本的每一对排序特征值进行逐一比较，而是构建一个直方图，仅考虑直方图边缘的分裂。例如，如果一个节点中有 100 个样本，它们的特征值可能被分为 10 组，这意味着只考虑 9 个可能的分裂，而不是 99 个。

我们可以按照如下方式实例化一个使用lossguide生长策略的 XGBoost 模型，使用学习率0.1，这是根据我们在前一个练习中进行的超参数探索的直觉得出的：

xgb_model_3 = xgb.XGBClassifier(
    n_estimators=1000,
    max_depth=0,
    learning_rate=0.1,
    verbosity=1,
    objective='binary:logistic',
    use_label_encoder=False,
    n_jobs=-1,
    tree_method='hist',
    grow_policy='lossguide')

请注意，我们已经设置了max_depth=0，因为该超参数与lossguide策略无关。相反，我们将设置一个名为max_leaves的超参数，它简单地控制将要生长的树的最大叶子数。我们将进行一个超参数搜索，范围从 5 到 100 个叶子：

max_leaves_values = list(range(5,105,5))
print(max_leaves_values[:5])
print(max_leaves_values[-5:])

这将输出如下内容：

[5, 10, 15, 20, 25]
[80, 85, 90, 95, 100]

现在我们准备好在这一系列超参数值范围内反复进行模型拟合和验证，类似于我们之前做过的操作：

%%time
val_aucs = []
for max_leaves in max_leaves_values:
    #Set parameter and fit model
    xgb_model_3.set_params(**{'max_leaves':max_leaves})
    xgb_model_3.fit(X_train, y_train, eval_set=eval_set,\
                    eval_metric='auc', verbose=False,\
                    early_stopping_rounds=30)
    #Get validation score
    val_set_pred_proba = xgb_model_3.predict_proba(X_val)[:,1]
    val_aucs.append(roc_auc_score(y_val, val_set_pred_proba))

输出将包括所有这些拟合的壁钟时间，在测试中大约是 24 秒。现在，让我们把结果放入数据框中：

max_leaves_df = \
pd.DataFrame({'Max leaves':max_leaves_values,
              'Validation AUC':val_aucs})

我们可以可视化验证集 AUC 随最大叶子数的变化，类似于我们对学习率的可视化：

mpl.rcParams['figure.dpi'] = 400
max_leaves_df.set_index('Max leaves').plot()

这将产生如下图所示的图形：

](tos-cn-i-73owjymdk6/7721619aaf324c79a3178e1b244f61de)

图 6.6：验证集 AUC 与 max_leaves 超参数的关系

较小的max_leaves值会限制生成的树的复杂度，这将理想地增加偏差，但也会减少方差，从而提高样本外的表现。当树的叶子数限制在 15 或 20 时，我们可以看到验证集 AUC 有所提高。那么，最大验证集 AUC 是多少呢？

max_auc = max_leaves_df['Validation AUC'].max()
max_auc

这应该输出以下内容：

0.8151200989120475

让我们确认最大验证 AUC 出现在max_leaves=20时，如图 6.6所示：

max_ix = max_leaves_df['Validation AUC'] == max_auc
max_leaves_df[max_ix]

这应该返回数据框的一行：

图 6.7：最优 max_leaves

通过使用lossguide增长策略，我们可以实现至少与我们迄今为止尝试过的任何方法一样好的性能。lossguide策略的一个关键优势是，对于较大的数据集，它能够提供比depthwise策略更快的训练速度，特别是在max_leaves较小的情况下。尽管这里的数据集足够小，这一速度差异并没有实际重要性，但在其他应用中，这种速度可能是理想的。

使用 SHAP 值解释模型预测

随着像 XGBoost 这样的前沿建模技术的发展，解释模型预测的实践在近年来有了显著的进展。到目前为止，我们已经了解到，逻辑回归的系数，或者随机森林中的特征重要性，可以为模型预测的原因提供洞察。2017 年，Scott Lundberg 和 Su-In Lee 在论文《统一模型预测解释方法》中描述了一种更强大的模型预测解释技术（arxiv.org/abs/1705.07874）。该技术被称为SHAP（Shapley 加法解释），它基于数学家 Lloyd Shapley 的早期工作。Shapley 发展了博弈论中的一个领域，用以理解玩家联盟如何为博弈的总体结果做出贡献。近期的机器学习研究在模型解释方面借鉴了这一概念，考虑了预测模型中的特征组或联盟如何贡献于模型的最终预测输出。通过考虑不同特征组的贡献，SHAP 方法能够孤立出单一特征的影响。

注意

在撰写本文时，第六章《梯度提升、XGBoost 和 SHAP 值》中使用的 SHAP 库与 Python 3.9 不兼容。因此，如果您使用的是 Python 3.9 作为基础环境，我们建议您按照前言中所述，设置 Python 3.8 环境。

使用 SHAP 值解释模型预测的一些显著特点包括：

• SHAP 值可以用来对模型预测做个性化解释；换句话说，可以通过 SHAP 理解单个样本的预测结果，分析每个特征的贡献。这与我们之前见过的随机森林特征重要性解释方法不同，后者仅考虑特征在模型训练集中的平均重要性。

• SHAP 值是相对于背景数据集计算的。默认情况下，这是训练数据集，当然也可以提供其他数据集。

• SHAP 值是可加的，这意味着对于单个样本的预测，SHAP 值可以加起来恢复预测值，例如预测的概率。

SHAP 方法有不同的实现，适用于各种类型的模型，这里我们将专注于树模型的 SHAP（Lundberg 等，2019 年，arxiv.org/abs/1802.03888），以便了解我们在验证集上的 XGBoost 模型预测。首先，我们将使用最优的max_leaves（即 20）重新拟合上一节中的xgb_model_3：

%%time
xgb_model_3.set_params(**{'max_leaves':20})
xgb_model_3.fit(X_train, y_train,\
                eval_set=eval_set,\
                eval_metric='auc',
                verbose=False,\
                early_stopping_rounds=30)

现在我们准备开始计算验证数据集的 SHAP 值。这里有 40 个特征和 1000 个样本：

X_val.shape

这将输出以下内容：

(1000, 40)

为了自动标注我们可以使用shap包绘制的图，我们将把验证集特征放入一个带列名的数据框中。我们将使用列表推导来生成通用的特征名称，例如“Feature 0, Feature 1, …”，并按如下方式创建数据框：

feature_names = ['Feature {number}'.format(number=number)
                 for number in range(X_val.shape[1])]
X_val_df = pd.DataFrame(data=X_val, columns=feature_names)
X_val_df.head()

dataframe的头部应如下所示：

图 6.8：验证特征的数据框

使用训练好的模型xgb_model_3和验证特征的数据框，我们准备创建一个explainer接口。SHAP 包有多种类型的解释器，我们将使用专门针对树模型的解释器：

explainer = shap.explainers.Tree(xgb_model_3, data=X_val_df)

这创建了一个使用模型验证数据作为背景数据集的解释器。现在我们已经准备好使用该解释器来获取 SHAP 值。SHAP 包使这变得非常简单。我们需要做的就是传入我们想要解释的数据集：

shap_values = explainer(X_val_df)

就是这样！那么，创建的这个变量shap_values是什么呢？如果你直接检查shap_values变量的内容，你会看到它包含三个属性。第一个是values，它包含 SHAP 值。让我们查看它的形状：

shap_values.values.shape

这将返回以下结果：

(1000, 40)

因为 SHAP 提供了个性化的解释，每个验证集中的 1,000 个样本都有一行数据。总共有 40 列，因为我们有 40 个特征，SHAP 值告诉我们每个特征对每个样本预测的贡献。shap_values 还包含一个 base_values 属性，即在考虑任何特征贡献之前的初始预测值，也定义为整个数据集的平均预测值。每个样本（1,000 个）都有一个这样的值。最后，还有一个 data 属性，包含特征值。所有这些信息可以通过不同的方式结合起来解释模型预测。

幸运的是，shap 包不仅提供了快速便捷的计算 SHAP 值的方法，还提供了一整套丰富的可视化技术。其中最受欢迎的一种是 SHAP 汇总图，它可视化每个特征对每个样本的贡献。我们来创建这个图表，然后理解其中展示的内容。请注意，大多数有趣的 SHAP 可视化使用了颜色，因此如果你正在阅读的是黑白版本，请参考 GitHub 仓库中的彩色图形：

mpl.rcParams['figure.dpi'] = 75
shap.summary_plot(shap_values.values, X_val_df)

这应该会生成以下内容：

图 6.9：合成数据验证集的 SHAP 汇总图

注意

如果你正在阅读本书的印刷版，你可以通过访问以下链接下载并浏览本章中部分图像的彩色版本：packt.link/ZFiYH

图 6.9 包含了大量信息，帮助我们解释模型。汇总图可能包含多达 40,000 个绘制点，每个特征对应一个点，每个样本（1,000 个验证样本）对应一个点（尽管默认情况下只显示前 20 个特征）。我们先从理解 x 轴开始。SHAP 值表示每个特征值对样本预测的加性贡献。这里显示的 SHAP 值是相对于预期值的，预期值即前面提到的 base_values。因此，如果某个特征对某个样本的预测影响较小，它将不会使预测偏离预期值太远，SHAP 值接近于零。然而，如果某个特征的影响较大，对于我们的二分类问题而言，这意味着预测的概率将被推向 0 或 1，SHAP 值会远离 0。负的 SHAP 值表示某个特征使得预测值更接近 0，正的 SHAP 值则表示更接近 1。

注意，图 6.9中显示的 SHAP 值不能直接解释为预测概率。默认情况下，XGBoost 二分类模型的 SHAP 值，使用 binary:logistic 目标函数计算并绘制，使用的是概率的对数赔率表示法，这在第三章中的逻辑回归细节和特征探索部分的*为什么逻辑回归被认为是线性模型？*小节中介绍过。这意味着 SHAP 值可以进行加减，换句话说，我们可以对其进行线性变换。

那么，图 6.9中点的颜色如何解释呢？这些颜色代表了每个样本的特征值，红色表示特征值较高，蓝色表示特征值较低。因此，例如，在图的第四行，我们可以看到特征 29 的高特征值（红点）对应的是最低的 SHAP 值。

点的垂直排列，换句话说，每个特征的点带宽度，表示在该位置上 x 轴上有多少个点。如果样本多，点带的宽度就会更大。

图中特征的垂直排列是根据特征重要性进行的。最重要的特征，也就是那些对模型预测具有最大平均影响（均值绝对 SHAP 值）的特征，排在列表的顶部。

虽然图 6.9中的总结图是查看所有最重要特征及其 SHAP 值的好方法，但它可能无法揭示一些有趣的关系。例如，最重要的特征——特征 3，似乎在特征值范围的中间部分有一大簇紫色点，这些点的 SHAP 值为正，而该特征的负 SHAP 值可能来自于特征值过高或过低。

这是什么情况呢？通常，当从 SHAP 总结图中看不出特征的影响时，我们使用的基于树的模型正在捕捉特征之间的交互效应。为了进一步了解单个特征及其与其他特征的交互作用，我们可以使用 SHAP 散点图。首先，我们绘制特征 3 的 SHAP 值散点图。注意，我们可以像索引数据框一样索引 shap_values 对象：

shap.plots.scatter(shap_values[:,'Feature 3'])

这将生成如下图所示：

](tos-cn-i-73owjymdk6/9ad96a0be7c24d828c2bb1e524ae5ab0)

图 6.10：特征 3 的 SHAP 值散点图

从图 6.10中，我们几乎可以得出与图 6.9总结图相同的信息：特征值在范围中间的部分对应较高的 SHAP 值，而极值部分的 SHAP 值较低。然而，scatter方法还允许我们通过另一个特征值来为散点图上的点上色，这样我们就可以看到特征之间是否存在交互作用。我们将使用第二重要特征——特征 5，为点上色：

shap.plots.scatter(shap_values[:,'Feature 3'],
                   color=shap_values[:,'Feature 5'])

生成的图应该像这样：

图 6.11：特征 3 的 SHAP 值散点图，按特征 5 的特征值着色。箭头 A 和 B 指示这些特征之间有趣的交互效应

图 6.11 显示了特征 3 和特征 5 之间的有趣交互。当样本位于特征 3 的特征值范围的中间时，也就是说，位于 图 6.11 中山丘形状的顶部时，从底部到顶部的点的颜色似乎变得越来越红（箭头 A）。这意味着，对于特征 3 范围中间的特征值，随着特征 5 的值增加，特征 3 的 SHAP 值也会增加。我们还可以看到，随着特征 3 的特征值沿 x 轴从中间向上升高，这种关系发生反转，特征 5 的较高特征值开始对应于特征 3 较低的 SHAP 值（箭头 B）。因此，与特征 5 的交互似乎对特征 3 的 SHAP 值有显著影响。

图 6.11 展示了复杂的关系，说明当存在交互效应时，增加特征值可能导致 SHAP 值增加或减少。图 6.11 中模式的具体原因与我们建模的合成数据集的创建有关，在此数据集中，我们在特征空间中指定了多个簇。如 第五章 中讨论的，决策树与随机森林，在 使用决策树：优势与预测概率 部分中提到，基于树的模型，如 XGBoost，能够有效地在多维特征空间中建模属于某一类别的点簇。SHAP 解释有助于我们理解模型如何进行这些表示。

在这里，我们使用了合成数据，特征没有现实世界的解释，因此我们无法为观察到的交互分配任何含义。然而，对于现实世界的数据，结合 SHAP 值和交互作用的详细探索可以提供有关模型如何表示客户或用户属性之间复杂关系的见解。例如，SHAP 值也很有用，因为它们可以提供相对于任何背景数据集的解释。尽管逻辑回归系数和随机森林的特征重要性完全由模型训练数据决定，SHAP 值可以为任何背景数据集计算；到目前为止，在本章中，我们一直使用的是验证数据。这为当预测模型部署到生产环境中时，提供了一个了解新预测是如何做出的机会。如果新预测的 SHAP 值与模型训练和测试数据的 SHAP 值非常不同，这可能表明传入数据的性质发生了变化，可能是时候考虑开发一个新模型了。在最后一章中，我们将讨论这些在实际应用中使用模型的实际问题。

练习 6.02：绘制 SHAP 交互作用、特征重要性，并从 SHAP 值中重构预测概率

在本练习中，你将更熟悉使用 SHAP 值来提供模型工作原理的可见性。首先，我们将再次查看特征 3 和特征 5 之间的交互作用，然后使用 SHAP 值计算特征重要性，类似于我们在第五章 决策树与随机森林中使用随机森林模型所做的那样。最后，我们将看到如何从 SHAP 值中获取模型输出，利用它们的加法性质：

注意

本练习的 Jupyter notebook 可以在packt.link/JcMoA找到。

1. 由于本节已经完成了初步步骤，我们可以再次查看特征 3 和特征 5 之间的交互作用，它们是合成数据集中的两个最重要的特征。使用以下代码制作图 6.11的另一个版本，不过这次我们关注的是特征 5 的 SHAP 值，并按照特征 3 的值进行着色：

   shap.plots.scatter(shap_values[:,'Feature 5'],
                      color=shap_values[:,'Feature 3'])
   

   结果图应如下所示：

   

   图 6.12：特征 5 的 SHAP 值散点图，按特征 3 的特征值着色

   与图 6.11不同，这里我们看到的是特征 5 的 SHAP 值。一般来说，从散点图来看，我们可以看到，随着特征 5 的特征值增加，SHAP 值也趋向增加。然而，确实有一些反例与这一普遍趋势相悖，并且与特征 3 有一个有趣的交互作用：对于特征 5 的给定值，可以将其视为图像中的一个垂直切片，点的颜色可以从下到上变得更红，表示负特征值，或者对于正特征值，点的颜色变得不那么红。这意味着对于特征 5 的给定值，它的 SHAP 值依赖于特征 3 的值。这进一步说明了特征 3 和 5 之间有趣的交互作用。在实际项目中，您选择展示哪个图表取决于您希望通过数据讲述什么样的故事，涉及特征 3 和 5 可能代表的实际世界中的量。

2. 使用以下代码创建特征重要性条形图：

   mpl.rcParams['figure.dpi'] = 75
   shap.summary_plot(shap_values.values, X_val, plot_type='bar')
   

   

   图 6.13: 使用 SHAP 值的特征重要性条形图

   特征重要性条形图提供了类似于第五章《决策树与随机森林》中的练习 5.03，拟合随机森林所获得的信息的可视化呈现：这是每个特征的一个数字，表示它对数据集的总体重要性。

   这些结果合理吗？回想一下，我们是通过创建包含三个有信息的特征和两个冗余特征的合成数据来实现的。在图 6.13中，似乎有四个特征的显著性比其他所有特征重要，因此，可能其中一个冗余特征的创建方式导致 XGBoost 经常选择它来进行节点分裂，而另一个冗余特征则使用较少。

   与我们在第五章《决策树与随机森林》中找到的特征重要性相比，这里的特征重要性有些不同。我们从 scikit-learn 中获得的随机森林模型的特征重要性是通过特征导致的节点杂质下降以及特征划分的训练样本的比例来计算的。相比之下，使用 SHAP 值计算的特征重要性是这样得到的：首先，取所有 SHAP 值（shap_values.values）的绝对值，然后对每个特征取所有样本的平均值，正如 x 轴标签所示。感兴趣的读者可以通过直接从 shap_values 计算这些指标来确认这一点。

   现在我们已经熟悉了 SHAP 值的多种用法，让我们看看它们的加法性质是如何允许重建预测概率的。

3. SHAP 值是相对于模型的期望值或基准值计算的。这可以解释为背景数据集中所有样本的平均预测值。然而，正如前面提到的，为了支持可加性，预测将以对数几率单位表示，而不是概率。可以通过以下方式从解释器对象中访问模型的期望值：

   explainer.expected_value
   

   输出应如下所示：

   -0.30949621941894295
   

   这条信息本身并不特别有用。然而，它为我们提供了一个基准，我们可以用它来重构预测概率。

4. 回顾一下，SHAP 值矩阵的形状是样本数和特征数。在我们对验证数据进行的练习中，形状将是 1,000 x 40。为了将每个样本的所有 SHAP 值相加，我们需要对列轴（axis=1）进行求和。这将把所有特征的贡献加在一起，实际上提供了与期望值的偏移量。如果我们将期望值加到其中，就得到了以下的预测值：

   shap_sum = shap_values.values.sum(axis=1) + explainer.expected_value
   shap_sum.shape
   

   这应该返回以下结果：

   (1000,)
   

   这意味着我们现在为每个样本有了一个单一的数值。然而，这些预测是在对数几率空间中。为了将它们转换为概率空间，我们需要应用在第三章，逻辑回归细节与特征探索中介绍的逻辑函数。

5. 如此应用逻辑转换到对数几率预测：

   shap_sum_prob = 1 / (1 + np.exp(-1 * shap_sum))
   

   现在，我们希望将通过 SHAP 值获得的预测概率与直接模型输出进行比较以确认。

6. 获取模型验证集的预测概率，并使用以下代码检查其形状：

   y_pred_proba = xgb_model_3.predict_proba(X_val)[:,1]
   y_pred_proba.shape
   

   输出应如下所示：

   (1000,)
   

   这与我们从 SHAP 导出的预测结果形状相同，如预期的那样。

7. 将模型输出和 SHAP 值的总和放在数据框中进行并排比较，并随机检查五行数据：

   df_check = pd.DataFrame(
       {'SHAP sum':shap_sum_prob,
        'Predicted probability':y_pred_proba})
   df_check.sample(5, random_state=1)
   

   输出应该确认这两种方法的结果是相同的：

   

   图 6.14：SHAP 导出的预测概率与直接从 XGBoost 获得的预测概率比较

   随机检查表明这五个样本具有相同的值。虽然由于机器算术的四舍五入误差，这些值可能不完全相等，但你可以使用 NumPy 的 allclose 函数，确保它们在用户配置的四舍五入误差范围内相同。

8. 确保 SHAP 导出的概率和模型输出的概率非常接近，如下所示：

   np.allclose(df_check['SHAP sum'],\
               df_check['Predicted probability'])
   

   输出应如下所示：

   True
   

   这表明这两列中的所有元素在四舍五入误差范围内是相等的。allclose 在出现四舍五入误差时非常有用，而精确相等（可以通过 np.array_equal 进行测试）通常不成立。

到目前为止，你应该对 SHAP 值在帮助理解机器学习模型方面的强大功能有了一些印象。SHAP 值的样本特定、个性化特性开启了非常详细的分析可能性，这可以帮助回答来自业务利益相关者的各种潜在问题，比如“模型会如何对像这样的人的数据做出预测？”或“为什么模型会对这个特定人做出这样的预测？”现在，我们已经熟悉了 XGBoost 和 SHAP 值，这两种先进的机器学习技术，接下来我们将回到案例研究数据上应用它们。

缺失数据

关于同时使用 XGBoost 和 SHAP 的最后一点说明，这两个包的一个宝贵特性是它们能够处理缺失值。回想一下在第一章，数据探索与清洗中，我们发现案例研究数据中一些样本在PAY_1特征上存在缺失值。到目前为止，我们的方法是，在构建模型时，简单地将这些样本从数据集中移除。这是因为，如果不以某种方式专门处理缺失值，scikit-learn 实现的机器学习模型将无法处理这些数据。忽略缺失值是一种方法，尽管这可能不令人满意，因为它涉及丢弃数据。如果缺失的数据仅占很小的比例，这可能没问题；然而，一般来说，知道如何处理缺失值是很重要的。

有几种方法可以填充缺失值，例如使用该特征非缺失值的均值或众数，或者随机选择一个非缺失值。你还可以构建一个模型，将该特征作为响应变量，所有其他特征作为该新模型的特征，然后预测缺失的特征值。这些方法在本书的第一版中有所探讨（packt.link/oLb6C）。然而，由于 XGBoost 通常在使用我们这里的表格数据进行二分类任务时，表现至少与其他机器学习模型一样好，并且能够处理缺失值，因此我们将不再深入探讨填充缺失值的问题，而是让 XGBoost 来为我们完成这项工作。

XGBoost 如何处理缺失数据？在每次有机会分裂节点时，XGBoost 仅考虑非缺失的特征值。如果一个特征包含缺失值且被选择用于分裂，缺失该特征值的样本会被发送到一个子节点，在最小化损失函数的基础上选择最优路径。

将 Python 变量保存到文件

在本章的活动中，为了进行文件读写，我们将使用新的 Python 语句（with）和pickle包。with语句使得处理文件更加方便，因为它们不仅打开文件，还会在使用完成后自动关闭文件，而无需用户单独进行这些操作。你可以使用如下代码片段将变量保存到文件中：

with open('filename.pkl', 'wb') as f:
    pickle.dump([var_1, var_2], f)

其中filename.pkl是你选择的文件路径，'wb'表示文件以二进制格式打开以供写入，pickle.dump将变量var_1和var_2保存到该文件中。要打开此文件并加载这些变量，可能需要在另一个 Jupyter Notebook 中，代码类似，但现在需要以二进制格式（'rb'）打开文件：

with open('filename.pkl', 'rb') as f:
    var_1, var_2 = pickle.load(f)

活动 6.01：使用 XGBoost 建模案例研究数据并使用 SHAP 解释模型

在本活动中，我们将利用本章所学的内容，使用一个合成数据集，并将其应用于案例研究数据。我们将观察 XGBoost 模型在验证集上的表现，并使用 SHAP 值解释模型预测。我们已通过替换先前忽略的、PAY_1特征缺失值的样本来准备数据集，同时保持没有缺失值的样本的训练/测试划分。你可以在本活动的笔记本附录中查看数据是如何准备的。

注意

包含解决方案以及附录的 Jupyter 笔记本可以在这里找到：packt.link/YFb4r。

1. 加载已为此练习准备好的案例研究数据。文件路径为../../Data/Activity_6_01_data.pkl，变量包括：features_response, X_train_all, y_train_all, X_test_all, y_test_all。

2. 定义一个验证集，用于训练 XGBoost 并进行早期停止。

3. 实例化一个 XGBoost 模型。使用lossguide增长策略，以便检查验证集在多个max_leaves值下的表现。

4. 创建一个max_leaves值的列表，范围从 5 到 200，步长为 5。

5. 创建用于早期停止的评估集。

6. 遍历超参数值并创建一个验证 ROC AUC 的列表，使用与练习 6.01：随机网格搜索调优 XGBoost 超参数相同的技术。

7. 创建一个超参数搜索结果的数据框，并绘制验证 AUC 与max_leaves的关系图。

8. 观察对应于验证集上最高 ROC AUC 的max_leaves数量。

9. 使用最佳超参数重新拟合 XGBoost 模型。这样我们就可以检查验证集的 SHAP 值，制作该数据的数据框。

10. 使用验证数据作为背景数据集，为我们的新模型创建 SHAP 解释器，获取 SHAP 值，并绘制总结图。

11. 绘制LIMIT_BAL SHAP 值的散点图，按与最强交互特征相关的颜色进行着色。

12. 将训练好的模型以及训练数据和测试数据保存到一个文件中。

    注意

    该活动的解决方案可以通过此链接找到。

总结

在本章中，我们学习了一些构建机器学习模型的前沿技术，尤其是针对表格数据。虽然其他类型的数据，如图像或文本数据，需要使用不同类型的模型（如神经网络）进行探索，但许多标准的商业应用仍然依赖于表格数据。XGBoost 和 SHAP 是一些最先进且流行的工具，你可以用它们来构建和理解这类数据的模型。在通过这些工具与合成数据进行实践并积累经验后，在接下来的活动中，我们将回到案例研究的数据集，看看如何使用 XGBoost 来对其建模，包括包含缺失特征值的样本，并利用 SHAP 值来理解模型。

第七章：7. 测试集分析、财务洞察和交付给客户

概述

本章介绍了几种分析模型测试集的技术，用于推导未来可能的模型表现。这些技术包括我们已经计算过的相同模型性能指标，例如 ROC AUC，以及一些新的可视化方法，如按预测概率分组的违约风险变化和预测概率的校准。阅读本章后，您将能够弥合机器学习的理论指标与商业世界财务指标之间的差距。您将能够识别关键洞察，并估算模型的财务影响，向客户提供如何实现这一影响的指导。最后，我们讨论了交付和部署模型时需要考虑的关键因素，如交付格式和监控模型使用情况的方法。

引言

在上一章中，我们使用了 XGBoost 进一步提高了模型的性能，超越了我们之前的所有努力，并学会了如何通过 SHAP 值解释模型预测。现在，我们将认为模型构建已完成，并解决交付给客户之前需要关注的剩余问题。本章的关键内容是对测试集的分析，包括财务分析，以及交付模型给希望在现实世界中使用它的客户时需要考虑的事项。

我们通过查看测试集来了解模型未来的表现。通过计算我们已知的指标，如 ROC AUC，但这次是针对测试集的，我们可以增强对模型在新数据上有用性的信心。我们还将学习一些直观的方法来可视化模型将客户分为不同违约风险等级的能力，比如十分位图。

您的客户可能会欣赏您在创建更准确的模型或更高 ROC AUC 模型方面所做的努力。然而，他们一定会更重视了解模型能够帮助他们赚取或节省多少钱，并且可能会乐于接受关于如何最大化模型潜力的具体指导。对测试集的财务分析可以模拟基于模型的不同策略场景，帮助客户选择适合他们的策略。

在完成财务分析后，我们将总结如何交付模型给客户使用，并讨论如何随着时间推移监控其表现。

模型结果回顾

为了开发满足客户商业需求的二分类模型，我们已经尝试了几种建模技术，取得了不同程度的成功。最终，我们希望选择性能最好的模型进行进一步分析，并呈现给客户。然而，向客户传达我们探索的其他选项也是很重要的，这样可以展示我们进行了一项彻底的研究项目。

在这里，我们回顾了针对案例研究问题所尝试的不同模型、需要调整的超参数以及交叉验证结果，或者在 XGBoost 的情况下使用的验证集。我们仅包括使用所有可能特征所做的工作，而不包括早期只使用一两个特征的探索性模型：

图 7.1：使用案例研究数据的建模活动总结

在向客户展示结果时，你应该准备好为各个技术背景层次的业务伙伴解释这些结果，包括那些技术背景很少的人。例如，业务伙伴可能不理解 ROC AUC 指标的推导过程；然而，这是一个重要的概念，因为它是我们用来评估模型的主要性能指标。你可能需要解释它是一个介于 0.5 和 1 之间变化的指标，并给出这些界限的直观解释：0.5 就像是抛硬币一样，而 1 是完美的，实际上几乎是无法达到的。

我们的结果介于两者之间，最佳模型接近 0.78。虽然一个给定模型的 ROC AUC 可能单独来看没有太大意义，图 7.1 显示我们尝试了几种方法，并且在最初的尝试基础上取得了性能的提升。最终，对于像案例研究这样的商业应用，像 ROC AUC 这样的抽象模型性能指标如果能配合财务分析，效果会更好。我们将在本章后面深入探讨这一点。

注：关于解释 ROC AUC

ROC AUC 分数的一个有趣解释是，对于两个样本，其中一个是正面结果，另一个是负面结果，正面样本的预测概率会比负面样本高。换句话说，对于评估数据集中所有可能的正负样本对，正面样本的模型预测高于负面样本的比例等于 ROC AUC。

从图 7.1中，我们可以看到，对于案例研究，通过工程化新特征来增强简单的逻辑回归模型或创建决策树集成模型来构建更复杂模型的努力，取得了更好的模型性能。特别是随机森林和 XGBoost 模型表现相似，尽管这些验证得分在技术上并不直接可比，因为在随机森林的情况下我们排除了缺失值，并且使用了 4 倍交叉验证，而在 XGBoost 中，缺失值被包含在内，且只有一个验证集用于提前停止。然而，图 7.1 表明，XGBoost 或随机森林可能是最好的选择。我们将在这里继续使用 XGBoost 模型。

现在我们已经决定了要交付的模型，考虑一下在模型开发过程中我们可能尝试过的其他方法也是好的。这些概念本书中不会展开，但你可能希望自己进行实验。

特征工程

提高模型性能的另一种方式，我们简要提到过的是LIMIT_BAL。与此特征最强交互的特征是两个月前的账单金额。尽管 XGBoost 可以找到类似的交互并在一定程度上对其进行建模，我们还可以构造一个新特征：过去每月账单金额与信用额度的比率，假设账单金额是账户的余额。这样计算得到的信用利用率可能是一个更强的特征，并且在这种方式下计算可能会比将信用额度和每月账单金额分别提供给模型的效果更好。

特征工程可能表现为通过操作现有特征来创建新特征，如之前的示例，或者它可能涉及引入完全新的数据源并用它们来创建特征。

新特征的灵感可能来自领域知识：与业务伙伴讨论他们认为可能是良好特征的内容会非常有帮助，尤其是在你对所从事的应用领域了解不如他们时。检查现有特征的交互作用也是假设新特征的一种方式，就像我们在活动 6.01中看到的与信用利用率相关的交互作用那样，使用 XGBoost 建模案例数据并用 SHAP 解释模型。

集成多个模型

在选择最终交付的案例研究项目模型时，交付随机森林或 XGBoost 中的任何一个都可能是可以接受的。机器学习中另一种常用的方法是集成多个模型。这意味着将不同模型的预测结果结合起来，类似于随机森林和 XGBoost 将多个决策树结合的方式。但在这种情况下，如何结合模型预测由数据科学家决定。创建模型集成的一种简单方式是取它们预测结果的平均值。

集成通常是在存在多个模型时进行的，可能是不同种类的模型或使用不同特征训练的模型，这些模型都具有良好的表现。在我们的案例中，可能使用随机森林和 XGBoost 的平均预测会比单独使用任何一个模型的表现更好。为了探索这一点，我们可以在验证集上比较表现，例如，用于 XGBoost 早停的验证集。

不同的建模技术

根据你为项目分配的时间和在不同建模技术方面的专业知识，你可能希望尝试尽可能多的方法。更高级的方法，例如用于分类的神经网络，可能在这个问题上提供更好的表现。我们鼓励你继续学习并掌握如何使用这些模型。然而，对于像我们这个案例研究中使用的表格数据，XGBoost 是一个非常好的默认选择，并且很可能提供优秀的表现，甚至可能是所有方法中最好的表现。

平衡类别

请注意，我们并没有处理响应变量中的类别不平衡。我们鼓励你尝试使用 scikit-learn 中的 class_weight='balanced' 选项或在 XGBoost 中使用 scale_pos_weight 超参数来拟合模型，以观察其效果。

尽管这些是进一步模型开发的有趣方向，但就本书而言，我们在此时已经完成了模型构建。我们将继续向前，检查 XGBoost 模型在测试集上的表现。

测试集上的模型表现

我们已经通过验证集对 XGBoost 模型的外部样本表现有了一些了解。然而，验证集在模型拟合过程中通过提前停止被使用。因此，我们能够做出的最严格的预期未来表现估计应该使用完全没有用于模型拟合的数据。这也是我们从模型构建过程中预留一个测试数据集的原因。

你可能会注意到，我们已经在一定程度上检查了测试集，例如，在第一章中评估数据质量和进行数据清理时。预测建模的金标准是在项目开始时预留出一个测试集，并且在模型完成之前不要对其进行任何检查。这是确保测试集中的任何知识没有在模型开发过程中“泄漏”到训练集中的最简单方法。当发生这种情况时，就有可能测试集不再能真实地代表未来的未知数据。然而，有时将所有数据一起探索和清理是很方便的，就像我们做的那样。如果测试数据与其余数据有相同的质量问题，那么就不会发生泄漏。最重要的是确保在决定使用哪些特征、拟合不同的模型并比较它们的表现时，不要查看测试集。

我们通过加载来自 活动 6.01、使用 XGBoost 构建案例研究数据模型并使用 SHAP 解释模型 中训练好的模型，连同训练数据、测试数据和特征名称，使用 Python 的 pickle 开始对测试集进行检查：

with open('../../Data/xgb_model_w_data.pkl', 'rb') as f:
    features_response, X_train_all, y_train_all, X_test_all,\
    y_test_all, xgb_model_4 = pickle.load(f)

在将这些变量加载到笔记本后，我们可以对测试集进行预测并进行分析。首先，获取测试集的预测概率：

test_set_pred_proba = xgb_model_4.predict_proba(X_test_all)[:,1]

现在，从 scikit-learn 导入 ROC AUC 计算函数，使用它来计算测试集的这一指标并显示出来：

from sklearn.metrics import roc_auc_score
test_auc = roc_auc_score(y_test_all, test_set_pred_proba)
test_auc

结果应如下所示：

0.7735528979671706

测试集上的 ROC AUC 为 0.774，略低于我们在验证集上看到的 XGBoost 模型的 0.779；然而，它们差异不大。由于模型拟合过程是针对验证集的性能进行了优化，因此在新数据上看到稍微低一点的表现并不完全令人意外。总体来说，测试性能符合预期，我们可以认为该模型在 ROC AUC 指标上已经成功测试。

虽然我们这里不会做这个，但在交付训练好的模型之前，最后一步可能是使用所有可用的数据，包括未见过的测试集，来拟合模型。这可以通过将训练数据和测试数据的特征（X_train_all，X_test_all）以及标签（y_train_all，y_test_all）进行拼接，使用它们来拟合一个新的模型，可能通过定义一个新的验证集来进行早停，或者使用当前的测试集来进行早停。这个方法的动机是机器学习模型通常在更多数据上训练时表现更好。缺点是，由于在这种情况下没有未见过的测试集，最终的模型可能被认为是未经测试的。

数据科学家对于采用哪种方法有不同的看法：是仅使用未见过的测试集进行模型评估，还是在完成所有前期步骤后，利用尽可能多的数据（包括测试集）来训练最终模型。一个考虑因素是模型是否会从更多数据中受益。这可以通过构建学习曲线来确定。虽然我们在这里不做说明，但学习曲线的概念是训练一个模型，使用不断增加的数据量，并计算在同一验证集上的验证分数。例如，如果你有 10,000 个训练样本，你可以留出 500 个作为验证集，然后先用前 1,000 个样本训练模型，再用前 2,000 个样本，依此类推，直到用完所有 9,500 个不在验证集中的样本。如果在更多数据上训练时，验证分数持续提高，甚至在使用所有可用数据时仍然有效，这表明使用更多数据训练模型会带来好处。然而，如果模型性能在某个点开始趋于平稳，并且似乎额外的数据并不能提升模型表现，那么你可能不需要这么做。学习曲线可以为如何使用测试集以及项目中是否需要更多数据提供指导。

就本案例研究而言，我们假设重新使用测试集来重新拟合模型不会带来任何益处。所以，现在我们主要关注的是向客户展示模型，帮助他们设计使用模型以实现业务目标的策略，并提供如何监控模型表现随时间变化的指导。

预测概率分布和十分位图

ROC AUC 指标很有帮助，因为它提供了一个总结模型在数据集上表现的单一数值。然而，观察模型在不同子集上的表现也是很有洞察力的。将数据集划分为不同子集的一种方式是使用模型预测结果。使用测试集，我们可以通过直方图可视化预测的概率：

mpl.rcParams['figure.dpi'] = 400
plt.hist(test_set_pred_proba, bins=50)
plt.xlabel('Predicted probability')
plt.ylabel('Number of samples')

这段代码应当生成如下图：

图 7.2：测试集的预测概率分布

测试集预测概率的直方图显示大多数预测集中在 [0, 0.2] 范围内。换句话说，根据模型，大多数借款人违约的概率在 0 到 20% 之间。然而，似乎有一小部分借款人存在更高的风险，集中在 0.7 附近。

检查模型在不同预测违约风险区域表现的直观方式是创建一个十分位图表，该图表根据预测概率的十分位将借款人分组。在每个十分位中，我们可以计算真实的违约率。我们预计会看到从最低预测十分位到最高预测十分位违约率的稳步上升。

我们可以像在 练习 6.01 中那样，使用 pandas 的 qcut 计算十分位数，随机网格搜索调优 XGBoost 超参数：

deciles, decile_bin_edges = pd.qcut(x=test_set_pred_proba,\
                                    q=10,\
                                    retbins=True)

在这里，我们正在拆分测试集的预测概率，这些预测概率通过 x 关键字参数提供。我们希望将它们分成 10 个大小相等的区间，预测概率最低的 10% 在第一个区间，以此类推，因此我们设置 q=10 分位数。然而，你可以将其分成任何数量的区间，比如 20（百分位）或 5（五分位）。由于我们设置了 retbins=True，区间的边界会保存在 decile_bin_edges 变量中，而十分位标签的序列则保存在 deciles 中。我们可以查看创建 10 个区间所需的 11 个区间边界：

decile_bin_edges

这应当产生如下图所示：

array([0.02213463, 0.06000734, 0.08155108, 0.10424594, 0.12708404,
       0.15019046, 0.18111563, 0.23032923, 0.32210371, 0.52585585,
       0.89491451])

为了使用 decile 序列，我们可以将其与测试集的真实标签和预测概率合并成一个 DataFrame：

test_set_df = pd.DataFrame({'Predicted probability':test_set_pred_proba,\
                            'Prediction decile':deciles,\
                            'Outcome':y_test_all})
test_set_df.head()

DataFrame 的前几行应如下所示：

图 7.3：包含预测概率和十分位数的 DataFrame

在 DataFrame 中，我们可以看到每个样本都被标注上了一个十分位区间，这个区间通过包含预测概率的区间边界来指示。结果显示了真实标签。我们希望在十分位图表中展示的是真实的违约率。为此，我们可以使用 pandas 的 groupby 功能。首先，我们通过对 decile 列进行分组来创建一个 groupby 对象：

test_set_gr = test_set_df.groupby('Prediction decile')

groupby对象可以通过其他列进行聚合。特别是在这里，我们关注的是每个十分位箱中的违约率，它是outcome变量的平均值。我们还计算了每个箱中数据的计数。由于分位数（如十分位）将总体分成相等大小的箱，因此我们预期计数应该相同或相似：

gr_df = test_set_gr.agg({'Outcome':['count', 'mean']})

查看我们的分组 DataFrame，gr_df：

图 7.4：测试集上预测概率的十个分位中的违约率

在图 7.4中，我们可以看到所有箱中的计数几乎相等。我们还可以看出，真实的违约率随着十分位的增加而增加，这是我们所期望的，因为我们知道我们的模型表现良好。在可视化数据之前，值得注意的是，这个 DataFrame 有一种特殊的列索引，叫做Outcome，以及第二级索引，标签分别为count和mean。访问具有多重索引的 DataFrame 的数据比我们之前使用的 DataFrame 要稍微复杂一些。我们可以通过以下方式显示列索引：

gr_df.columns

这将生成以下结果：

MultiIndex([('Outcome', 'count'),
            ('Outcome',  'mean')],
           )

在这里我们可以看到，要访问多重索引中的某一列，我们需要使用元组来指定索引的每一层级，例如，gr_df[('Outcome','count')]。虽然在这里 MultiIndex 并不是必需的，因为我们只对一列（Outcome）进行了聚合，但在对多个列进行聚合时，它会变得非常有用。

现在我们想要创建一个可视化图表，展示模型的预测如何将借款人准确地分组，并且这些组的违约风险不断上升。我们将展示每个箱中的计数，以及每个箱中的违约风险。由于这些列的量级不同，计数在数百之间，风险在 0 到 1 之间，因此我们应该使用双y轴图。为了对图表外观有更多的控制，我们将使用 Matplotlib 函数来创建这个图，而不是通过 pandas 来做。首先，我们绘制每个箱中样本数量的图，并用与图表相同的颜色标注y轴刻度，确保清晰。请参阅 GitHub 上的笔记本，如果你是以黑白阅读的，因为颜色对该图很重要。这段代码应该与接下来的代码段在同一个单元格中运行。这里我们创建了一组坐标轴，然后添加了图表并进行了格式化和标注：

ax_1 = plt.axes()
color_1 = 'tab:blue'
gr_df[('Outcome', 'count')].plot.bar(ax=ax_1, color=color_1)
ax_1.set_ylabel('Count of observations', color=color_1)
ax_1.tick_params(axis='y', labelcolor=color_1)
ax_1.tick_params(axis='x', labelrotation = 45)

请注意，我们正在为样本大小创建一个bar图。我们希望在此基础上添加一条线形图，显示每个箱中的违约率，且该线图使用右侧的y轴，但与现有图表共享相同的x轴。Matplotlib 为此目的提供了一个叫做twinx的方法，该方法可以在axes对象上调用，返回一个共享相同x轴的新坐标轴对象。我们采取类似的步骤，然后绘制违约率并进行标注：

ax_2 = ax_1.twinx()
color_2 = 'tab:red'
gr_df[('Outcome', 'mean')].plot(ax=ax_2, color=color_2)
ax_2.set_ylabel('Default rate', color=color_2)
ax_2.tick_params(axis='y', labelcolor=color_2)

在运行前两个代码片段后，应该会出现以下图形：

图 7.5：根据模型预测十分位数的违约率

图 7.5 包含与 图 7.4 中显示的数据框相同的信息，但呈现方式更为美观。很明显，违约风险随着每个十分位数的增加而增加，风险最大的 10%借款人的违约率接近 70%，而风险最小的借款人违约率低于 10%。当一个模型能够有效地区分出违约风险持续增加的借款人群体时，称该模型倾斜了所研究的总体。还要注意，违约率在最低的 5 到 7 个十分位数之间相对平坦，这可能是因为这些观察值大多集中在预测风险范围[0, 0.2]中，如图 7.2中的直方图所示。

将测试集分成等人口十分位数是评估模型性能的一种方式，特别是在违约风险倾斜方面。然而，客户可能希望查看按不同群体划分的违约率，例如按等间隔箱（例如，将所有预测范围为 0, 0.2)、[0.2, 0.4)等的观察值放在一起，而不管每个箱中的样本大小），或以其他方式。你将在接下来的练习中探索如何在 pandas 中轻松实现这一点。

在接下来的练习中，我们将使用一些统计学概念来帮助创建误差条，包括我们之前学到的均值的标准误差和二项分布的正态近似。

我们从第五章，决策树和随机森林中知道，我们可以估计样本均值的方差为 ![1，其中 n 是样本大小， 是一个理论上更大总体的未观察到的方差。虽然我们不知道 ，但可以通过我们观察到的样本方差来估计。对于二元变量，样本方差可以计算为 p(1-p)，其中 p 是成功的比例，或者说是案例研究中的违约比例。根据上面的样本均值方差公式，我们可以代入观察到的方差，然后取平方根得到均值的标准误差：。在某些情况下，这个公式也被称为二项分布的正态近似。我们将在下面使用它来为不同模型预测区间的违约率的等间隔图表创建误差条。有关这些概念的更多细节，建议查阅统计学教材。

练习 7.01：等间隔图表

在这个练习中，你将制作一个类似于图 7.5所示的图表；然而，与你将测试集拆分为预测概率的等人口十分位不同，你将使用预测概率的等间隔区间。如果业务伙伴希望使用某些得分范围来思考潜在的基于模型的策略，指定这些区间可能会有帮助。你可以使用 pandas 的 cut 来创建等间隔区间，或者使用一个区间边界数组创建自定义区间，这类似于你使用 qcut 来创建分位标签的方式：

注意

你可以在 packt.link/4Ev3n 找到这个练习的 Jupyter notebook。

1. 使用以下代码创建等间隔标签系列，针对 5 个区间：

   equal_intervals, equal_interval_bin_edges = \
       pd.cut(x=test_set_pred_proba,\
              bins=5,\
              retbins=True)
   

   请注意，这与调用 qcut 类似，不同的是在这里使用 cut 时，我们可以通过向 bins 参数传递一个整数来指定我们想要的等间隔区间数。你也可以为此参数提供一个数组来指定自定义的区间边界。

2. 使用以下代码检查等间隔区间边界：

   equal_interval_bin_edges
   

   结果应如下所示：

   array([0.02126185, 0.1966906 , 0.37124658, 0.54580256, 0.72035853,
          0.89491451])
   

   你可以通过从第一个元素到倒数第二个元素的子数组与从第二个元素开始到最后一个元素的子数组相减，来确认这些区间边界之间是等间隔的。

3. 使用以下代码检查区间边界之间的间隔：

   equal_interval_bin_edges[1:] - equal_interval_bin_edges[:-1]
   

   结果应如下所示：

   array([0.17542876, 0.17455598, 0.17455598, 0.17455598, 0.17455598])
   

   你可以看到区间边界之间的距离大致相等。第一个区间边界比最小的预测概率稍小，你可以自行确认这一点。

   为了创建一个类似于图 7.5的图表，首先我们需要将区间标签与响应变量放入一个 DataFrame，就像我们之前使用十分位标签时做的那样。我们还将预测概率放入 DataFrame 以供参考。

4. 创建一个包含预测概率、区间标签和测试集响应变量的 DataFrame，如下所示：

   test_set_bins_df =\
   pd.DataFrame({'Predicted probability':test_set_pred_proba,\
                 'Prediction bin':equal_intervals,\
                 'Outcome':y_test_all})
   test_set_bins_df.head()
   

   结果应如下所示：

   

   图 7.6：包含等间隔区间的 DataFrame

   我们可以使用这个 DataFrame 按照区间标签进行分组，然后获取我们感兴趣的指标：表示默认率和每个区间内样本数量的聚合数据。

5. 使用以下代码按区间标签分组并计算各区间内的默认率和样本数：

   test_set_equal_gr = test_set_bins_df.groupby('Prediction bin')
   gr_eq_df = test_set_equal_gr.agg({'Outcome':['count', 'mean']})
   gr_eq_df
   

   结果的 DataFrame 应该如下所示：

   

   图 7.7：五个等间隔区间的分组数据

   请注意，与分位数不同，这里每个区间内的样本数是不同的。默认率在各区间之间呈现一致的增长趋势。让我们绘制这个 DataFrame，创建一个类似于图 7.5的可视化。

   在创建此可视化之前，为了考虑到由于这些范围内的样本量减少，高预测概率的违约率估计可能不够稳健，我们将计算违约率的标准误差。

6. 使用以下代码计算箱体内违约率的标准误差：

   p = gr_eq_df[('Outcome', 'mean')].values
   n = gr_eq_df[('Outcome', 'count')].values
   std_err = np.sqrt(p * (1-p) / n)
   std_err
   

   结果应该如下所示：

   array([0.00506582, 0.01258848, 0.02528987, 0.02762643, 0.02683029])
   

   注意，对于那些高分范围且样本较少的箱体，标准误差较大。将这些标准误差与违约率一起可视化会非常有帮助。

7. 使用以下代码创建一个违约率与样本大小的等间隔图。该代码与 图 7.5 所需的代码非常相似，不同之处在于这里我们使用 yerr 关键字和上一阶段的结果，在违约率图上加入了误差条：

   ax_1 = plt.axes()
   color_1 = 'tab:blue'
   gr_eq_df[('Outcome', 'count')].plot.bar(ax=ax_1, color=color_1)
   ax_1.set_ylabel('Count of observations', color=color_1)
   ax_1.tick_params(axis='y', labelcolor=color_1)
   ax_1.tick_params(axis='x', labelrotation = 45)
   ax_2 = ax_1.twinx()
   color_2 = 'tab:red'
   gr_eq_df[('Outcome', 'mean')].plot(ax=ax_2, color=color_2,
                                      yerr=std_err)
   ax_2.set_ylabel('Default rate', color=color_2)
   ax_2.tick_params(axis='y', labelcolor=color_2)
   

   结果应该如下所示：

   

图 7.8：等间隔箱体中的违约率与样本数量的图

我们可以在 图 7.8 中看到，不同箱体中的样本数量差异较大，与分位数方法形成对比。尽管高分数箱体中的样本相对较少，导致较大的标准误差，但违约率图上的误差条仍然较小，相比于从低分到高分箱体的总体趋势，我们可以对这一趋势有信心。

预测概率的校准

图 7.8 的一个有趣特点是，违约率的折线图从一个箱体到下一个箱体大致增加相同的量。与 图 7.5 中的十分位图对比，违约率最初增加缓慢，然后迅速增加。还要注意，违约率大致是每个箱体内预测概率边缘的中点。这意味着违约率与每个箱体的平均模型预测相似。换句话说，我们的模型不仅能有效地按照从低到高的违约风险对借款人进行排序（通过 ROC AUC 量化），而且还似乎能够准确预测违约概率。

衡量预测概率与实际概率的匹配程度是校准****概率的目标。概率校准的标准度量遵循上述讨论的概念，称为期望校准误差（ECE），定义为：

图 7.9：期望校准误差

其中，索引 i 从 1 到箱体数目（N）范围，Fi 是所有样本中落入箱体 i 的比例，oi 是箱体 i 中为正样本（即在案例研究中为违约者）的比例，ei 是箱体 i 内预测概率的平均值。

我们可以使用一个与图 7.4中非常相似的 DataFrame 来计算测试集中分位区间内预测概率的 ECE，这个 DataFrame 用于创建分位图。我们唯一需要添加的是每个区间内的平均预测概率。按如下方式创建这样的 DataFrame：

cal_df = test_set_gr.agg({'Outcome':['count', 'mean'],\
                          'Predicted probability':'mean'})
cal_df

输出的 DataFrame 应如下所示：

图 7.10：计算 ECE 指标的 DataFrame

为了方便起见，我们定义一个变量F，表示每个区间中样本的比例。这是从上面的 DataFrame 中每个区间的计数除以样本总数，后者是从测试集响应变量的形状中获得的：

F = cal_df[('Outcome', 'count')].values/y_test_all.shape[0]
F

输出应为：

array([0.10003368, 0.10003368, 0.10003368, 0.09986527, 0.10003368,
       0.10003368, 0.09986527, 0.10003368, 0.10003368, 0.10003368])

因此，每个区间大约包含 10%的样本。这是预期的，当然，因为区间是使用分位数方法创建的。然而，对于其他分箱方法，区间中的样本数量可能不相等。现在让我们在代码中实现 ECE 的公式来计算这个指标：

ECE = np.sum(
    F
    * np.abs(
             cal_df[('Outcome', 'mean')]
             - cal_df[('Predicted probability', 'mean')]))
ECE

输出应为：

0.008144502190176022

这个数字代表我们最终模型在测试集上的 ECE。单独来看，这个数字并没有太多意义。然而，像这样的指标可以在时间上进行监控，在模型投入生产并在实际应用中使用之后。如果 ECE 开始增加，这表示模型的校准性变差，可能需要重新训练，或者对输出应用校准过程。

检查我们预测概率的校准的一个更直观方法是绘制 ECE 所需的各个成分，特别是响应变量的真实违约率，与每个区间内模型预测的平均值进行对比。在此基础上，我们加上一条 1-1 线，表示完美校准，作为参考点：

ax = plt.axes()
ax.plot([0, 0.8], [0, 0.8], 'k--', linewidth=1,
        label='Perfect calibration')
ax.plot(cal_df[('Outcome', 'mean')],\
        cal_df[('Predicted probability', 'mean')],\
        marker='x',\
        label='Model calibration on test set')
ax.set_xlabel('True default rate in bin')
ax.set_ylabel('Average model prediction in bin')
ax.legend()

结果图应如下所示：

图 7.11：预测概率的校准图

图 7.11显示了模型预测的概率与真实的违约率非常接近，因此模型似乎被良好校准。为了获得更多见解，你可以尝试自己向这个图添加误差条作为练习。还要注意，scikit-learn 提供了一个函数来计算创建图 7.11所需的信息：sklearn.calibration.calibration_curve。然而，这个函数并不会返回每个区间的样本大小。

对于概率校准，另一个需要注意的点是，一些处理类别不平衡的方法，如过采样或欠采样，会改变训练数据集中类别的比例，这会影响预测的概率，并可能使其不那么准确。不过，这一点可能不那么重要，特别是与模型根据借款人违约风险进行排名的能力（通过 ROC AUC 衡量）相比，这取决于客户的需求。

财务分析

我们目前计算的模型性能指标基于一些抽象的度量，可以应用于分析任何分类模型：模型的准确度、模型在不同阈值下识别真阳性与假阳性能力的准确性（ROC AUC）、正向预测的正确性（精确度），或诸如斜坡风险等直观的度量。这些指标对于理解模型的基本工作原理非常重要，并且在机器学习领域广泛使用，因此理解它们非常关键。然而，在将模型应用于商业案例时，我们并不能总是直接使用这些性能指标来制定如何使用模型来指导商业决策或评估模型可能创造的价值的策略。为了更进一步地将预测概率和阈值的数学世界与成本和收益的商业世界联系起来，通常需要进行某种财务分析。

为了帮助客户进行此分析，数据科学家需要了解根据模型预测可能采取的决策和行动。这应该是与客户的对话内容，最好在项目生命周期的早期进行。我们将其留到本书的最后，以便建立对预测建模是什么以及如何工作的基础理解。然而，在项目开始时了解模型使用的商业背景，可以让你为模型性能设定目标，并在整个项目过程中追踪这些目标，就像我们追踪不同模型的 ROC AUC 一样。将模型性能指标转化为财务术语是本节的主题。

对于像案例研究中的二分类模型，数据科学家需要了解以下几个问题的答案，以帮助客户弄清楚如何使用该模型：

• 客户希望通过模型帮助做出哪些决策？

• 如何利用二分类模型的预测概率来帮助做出这些决策？

• 它们是是/否决策吗？如果是，那么选择一个单一的预测概率阈值就足够了。

• 是否会根据模型结果做出超过两级的活动决策？如果是，那么选择两个或更多的阈值，例如将预测结果分为低、中、高风险，可能是解决方案。例如，预测概率低于 0.5 可以视为低风险，0.5 到 0.75 之间为中风险，0.75 以上为高风险。

• 根据模型指导，采取不同行动方案的成本是什么？

• 从模型指导下采取成功行动可能获得的潜在好处是什么？

与客户的财务对话

我们向案例研究客户询问了上述要点，得知以下信息：对于那些高风险违约的信贷账户，客户正在设计一个新项目，为账户持有人提供个性化的咨询，鼓励他们按时支付账单，或在无法按时支付的情况下提供替代付款选项。信贷咨询由经过培训的客户服务代表在呼叫中心进行。每次咨询的费用为新台币 7,500 元，每次咨询的预期成功率为 70%，意味着平均来说，70% 接到提供咨询电话的客户会按时支付账单，或者采取对债权人可接受的替代安排。成功咨询的潜在收益在于，如果账户本应违约但在咨询后未违约，则账户的每月账单金额将视为节省。目前，违约账户的每月账单会被报告为损失。

在与客户进行前述对话后，我们已经获得了进行财务分析所需的材料。客户希望我们帮助他们决定应该联系哪些成员并提供信贷咨询。如果我们能帮助他们缩小需要联系进行咨询的人员名单，就能通过避免不必要和昂贵的联系方式为他们节省费用。客户在信贷咨询上的有限资源将更合适地用于高风险违约账户。这应该能够通过预防违约带来更大的节省。此外，客户告诉我们，如果我们能给他们一个关于值得提供多少次咨询的概念，我们的分析可以帮助他们申请咨询项目的预算。

在进行财务分析时，我们看到模型将帮助客户做出的决策是逐个账户进行的“是/否”决策：是否为给定账户的持有人提供咨询。因此，我们的分析应该专注于找到一个合适的预测概率阈值，通过该阈值我们可以将账户分为两组：高风险账户将接受咨询，低风险账户则不接受。

练习 7.02: 成本与节省的描述

模型输出与客户将做出的业务决策之间的联系，归结为为预测概率选择一个阈值。因此，在本次练习中，我们将描述咨询程序的预期成本（以提供单独咨询会议的成本表示）以及预期节省（以预防违约的节省表示），这些将在一系列阈值下进行计算。每个阈值下会有不同的成本和节省，因为每个阈值预计会导致不同数量的正预测结果，并且在这些结果中会有不同数量的真正正例。第一步是创建一个潜在阈值的数组。我们将使用从 0 到 1，步长为 0.01 的值。执行以下步骤以完成本次练习：

注意

本次练习的 Jupyter notebook 可以在这里找到：packt.link/yiMEr。基于本章前面的结果，已在 notebook 中添加了准备数据的额外步骤。请确保在进行本次练习前，先执行 notebook 中展示的前提步骤。

1. 使用以下代码创建一个阈值范围，以计算咨询的预期成本和收益：

   thresholds = np.linspace(0, 1, 101)
   

   这会在 0 到 1 之间（包括 0 和 1）创建 101 个等间距的点。

   现在，我们需要了解预防违约的潜在节省。为了精确计算这一点，我们需要知道下个月的月度账单。然而，客户已经告知我们，在他们需要创建待联系账户持有者名单时，这些数据将无法提供。因此，为了估算潜在节省，我们将使用最近的月度账单。

   我们将使用测试数据来创建此分析，因为它提供了模型交付给客户后使用的模拟：在未用于模型训练的新账户上。

2. 确认测试数据特征数组中对应最近一个月账单的索引：

   features_response[5]
   

   输出应为：

   'BILL_AMT1'
   

   索引 5 对应的是最近几个月的账单，我们稍后会用到。

3. 将咨询成本存储在一个变量中，以便用于分析：

   cost_per_counseling = 7500
   

   我们还知道，客户告知我们，咨询程序的效果并不是 100%的。我们应该在分析中考虑这一点。

4. 存储客户提供的有效性率，以便用于分析：

   effectiveness = 0.70
   

   现在，我们将计算每个阈值的成本和节省。我们会逐步讲解每个计算过程，但目前我们需要创建空数组来存储每个阈值的结果。

5. 创建空数组来存储分析结果。我们将在后续步骤中解释每个数组的作用：

   n_pos_pred = np.empty_like(thresholds)
   total_cost = np.empty_like(thresholds)
   n_true_pos = np.empty_like(thresholds)
   total_savings = np.empty_like(thresholds)
   

   这些创建了与我们分析中的阈值数量相同的空数组。我们将遍历每个阈值值来填充这些数组。

6. 创建一个counter变量并打开一个for循环来遍历各个阈值：

   counter = 0
   for threshold in thresholds:
   

   对于每个阈值，根据预测概率高于该阈值的数量，会有不同数量的正向预测。这些预测对应的是被认为会违约的账户。每个被预测为违约的账户都会接到一次咨询电话，这会产生相应的费用。因此，这是成本计算的第一部分。

7. 确定在该阈值下哪些账户获得了正向预测：

       pos_pred = test_set_pred_proba > threshold
   

   pos_pred 是一个布尔数组。pos_pred 的总和表示在该阈值下预测的违约数量。

8. 计算给定阈值下的正向预测数量：

       n_pos_pred[counter] = sum(pos_pred)
   

9. 计算给定阈值下的咨询总成本：

       total_cost[counter] \
           = n_pos_pred[counter] * cost_per_counseling
   

   现在，我们已经了解了咨询项目在每个阈值下的可能成本，我们需要查看预期的节省。在提供咨询服务给正确的账户持有人时，会获得节省：即那些本来会违约的账户。从分类问题的角度来看，这些是正向预测，其响应变量的真实值也是正向的——换句话说，是真正的正向预测。

10. 基于正向预测数组和响应变量，确定哪些账户是真正的正向预测：

        true_pos = pos_pred & y_test_all.astype(bool)
    

11. 通过对真正的正向预测数组求和来计算真正的正向预测数量：

        n_true_pos[counter] = sum(true_pos)
    

    我们从成功为原本会违约的账户持有人提供咨询所能获得的节省，取决于每个预防违约的节省金额以及咨询的有效性率。我们无法预防每一个违约。

12. 使用真正的正向预测数量、预防违约的节省（通过上个月的账单估算）和咨询的有效性率来计算每个阈值下的预期节省：

        total_savings[counter] = np.sum(
            true_pos.astype(int)
            * X_test_all[:,5]
            * effectiveness
            ) 
    

13. 增加计数器：

    counter += 1
    

    步骤 5 到 13 应该作为 for 循环在 Jupyter Notebook 的一个单元格中运行。之后，可以通过将节省减去成本来计算每个阈值的净节省。

14. 通过将节省和成本数组相减，计算所有阈值的净节省：

    net_savings = total_savings - total_cost
    

    现在，我们可以可视化通过为合适的账户持有人提供咨询服务，我们可能帮助客户节省多少钱。让我们来可视化一下。

15. 按照如下方式绘制净节省与阈值的关系：

    mpl.rcParams['figure.dpi'] = 400
    plt.plot(thresholds, net_savings)
    plt.xlabel('Threshold')
    plt.ylabel('Net savings (NT$)')
    plt.xticks(np.linspace(0,1,11))
    plt.grid(True)
    

    结果图应如下所示：

    

    图 7.12：净节省与阈值的关系图

    图表显示，阈值的选择非常重要。虽然在许多不同的阈值下都能实现净节省，但看起来通过将阈值设置在大约 0.25 到 0.5 之间的某个范围内，可以获得最大的净节省。

    让我们确认产生最大节省的最佳阈值，并查看节省的具体数额。

16. 使用 NumPy 的argmax查找净节省数组中最大元素的索引：

    max_savings_ix = np.argmax(net_savings)
    

17. 显示产生最大净节省的阈值：

    thresholds[max_savings_ix]
    

    输出应如下所示：

    0.36
    

18. 显示最大的净节省：

    net_savings[max_savings_ix]
    

    输出应如下所示：

    13415710.0
    

我们发现，最大净节省发生在 0.36 的阈值下。在这个阈值下，所实现的净节省金额超过 1300 万新台币，针对这个测试数据集中的账户。这些节省需要根据客户所服务的账户数量进行规模化估算，以推算出总的可能节省金额，前提是我们所使用的数据能够代表所有这些账户。

然而需要注意的是，节省金额在约 0.5 的阈值之前大致相同，如图 7.12所示。

随着阈值的增加，我们在“提高标准”，即要求客户的风险更高，才能与他们联系并提供咨询服务。从 0.36 提高到 0.5 意味着我们只会联系那些风险较高、概率大于 0.5 的客户。这意味着联系的客户较少，从而减少了项目的前期成本。图 7.12 表明，即使我们联系较少的客户，可能仍然能够创造大致相同的净节省。虽然净效应相同，但咨询的初始支出较小。这对客户来说可能是更可取的。我们将在接下来的活动中进一步探讨这个概念。

活动 7.01：得出财务见解

财务分析的原始数据已完成。然而，在本活动中，你的目标是从这些结果中生成一些额外的见解，为客户提供更多的背景信息，帮助他们理解我们构建的预测模型如何为他们创造价值。特别地，我们已经查看了模型构建中保留的测试集结果。客户可能拥有比他们提供给我们的更多的账户，这些账户能代表他们的业务。你应该向他们报告可以轻松扩展到他们业务规模的结果，具体来说，是以账户数量为基础。

我们还可以帮助客户了解这个项目的成本；虽然净节省是需要考虑的一个重要数字，但客户必须在实现这些节省之前为咨询项目提供资金。最后，我们将把财务分析与标准的机器学习模型性能指标关联起来。

一旦完成活动，你应该能够向客户传达咨询项目的初步成本，并获得诸如下图所示的精确率和召回率的图表：

图 7.13：预期的精确率-召回率曲线

这个曲线在解释模型在不同阈值下创造的价值时非常有用。

执行以下步骤以完成活动：

注意

包含此活动代码的 Jupyter 笔记本可以在此处找到：packt.link/2kTVB。根据本章先前的结果，已向笔记本添加了为此活动准备数据的额外步骤。请按照笔记本中呈现的先决步骤执行这些步骤。

1. 使用测试集，计算如果没有辅导计划，所有违约的成本。

2. 计算辅导计划可以减少违约成本的百分比。

3. 计算在最优阈值下每个帐户的净储蓄，考虑可能能够辅导的所有帐户，换句话说，相对于整个测试集。

4. 绘制每个阈值下每个帐户的净储蓄与每个帐户的辅导成本。

5. 绘制在每个阈值下预测为正的帐户的比例（这称为“标志率”）。

6. 绘制测试数据的精确率-召回率曲线。

7. 在y-轴上分别绘制精确率和召回率，对应x-轴上的阈值。

   注意

   通过此链接可找到此活动的解决方案。

关于向客户交付预测模型的最终思考

我们现在已经完成了建模活动，并创建了财务分析，以指示客户如何使用模型。虽然我们已完成了数据科学家的基本智力贡献，但需要与客户达成一致意见，以确定所有这些贡献的交付形式。

关键贡献是训练模型中体现的预测能力。假设客户可以使用我们用 XGBoost 创建的训练模型对象，此模型可以像我们所做的那样保存到磁盘并发送给客户。然后，客户可以在其工作流程中使用它。这种模型交付路径可能需要数据科学家与客户组织的工程师合作，以在客户基础设施中部署模型。

或者，可能需要将模型表达为数学方程（例如，使用逻辑回归系数）或一组 if-then 语句（如决策树或随机森林），客户可以使用 SQL 实现预测能力。虽然由于可能存在许多树和许多级别，将随机森林表达为 SQL 代码是很繁琐的，但有软件包可以根据训练好的 scikit-learn 模型为您创建此表示（例如，pypi.org/project/SKompiler/）。

注意：云平台用于模型开发和部署

在本书中，我们使用 scikit-learn 和 XGBoost 包在本地计算机上构建预测模型。最近，云平台如亚马逊网络服务（AWS）通过 Amazon SageMaker 等服务提供了机器学习能力。SageMaker 包含了一个 XGBoost 的版本，您可以使用类似于我们在这里所做的语法来训练模型。在本书中展示的方法和亚马逊 SageMaker 的模型训练实现之间可能存在细微差异，建议您在每一步检查您的工作，以确保结果符合预期。例如，在 SageMaker 中使用早停止来拟合 XGBoost 模型可能需要额外的步骤，以确保训练的模型使用最佳迭代进行预测，而不是在训练停止时的最后迭代。

云平台如 AWS 很有吸引力，因为它们可能极大地简化将训练好的机器学习模型集成到客户的技术堆栈中的过程，而这在许多情况下可能已经建立在云平台上。

在使用模型进行预测之前，客户需要确保数据的准备方式与我们之前构建模型时相同。例如，删除所有特征值为0的样本以及清理EDUCATION和MARRIAGE特征的操作必须与我们在本章前面演示的方式相同。另外，还有其他可能的交付模型预测的方式，比如客户将特征提供给数据科学家，然后接收预测结果。

讨论交付内容的另一个重要考虑因素是：预测结果应该以什么格式交付？ 二元分类模型的典型交付格式，例如我们为案例研究创建的模型，是按照预测违约概率对账户进行排名。预测的概率应该与账户 ID 以及客户希望的其他列一起提供。这样，当呼叫中心逐个联系账户持有人以提供咨询时，他们可以首先联系那些违约风险最高的人，然后根据时间和资源的允许继续联系优先级较低的账户持有人。客户应该被告知使用哪个预测概率阈值，以获得最高的净节省。这个阈值将代表按照违约概率排名的账户持有人联系的终止点。

模型监控

根据客户与数据科学家的合作时间，持续监控模型性能随时间的变化始终是有益的，因为它在使用过程中是否保持稳定或逐渐下降？在评估此案例时，需要牢记，如果账户持有人正在接受辅导，他们的违约概率可能会低于预测的概率，因为新辅导计划的预期效果。因此，为了测试辅导计划的有效性，最好将一部分随机选择的账户持有人保留下来，这部分人将不会接受任何辅导，无论他们的违约风险如何。这部分人被称为对照组，它应该相对于接受辅导的账户群体较小，但足够大以便得出统计学上显著的推论。

尽管本书的范围并不包括如何设计和使用对照组的细节，但在这里可以简要说明，模型的预测能力可以通过对照组进行评估，因为他们没有接受任何辅导，类似于模型训练时使用的账户群体。对照组的另一个好处是，可以将违约率及因违约而产生的财务损失与那些接受了模型引导辅导计划的账户进行比较。如果该计划按预期运行，接受辅导的账户应具有较低的违约率和较小的违约财务损失。对照组可以提供证据，证明该计划实际上是有效的。

注意：选择性治疗的高级建模技术——提升建模

当一家企业考虑是否有选择地为其客户提供昂贵的治疗方案（如案例研究中的辅导计划）时，应考虑使用一种被称为提升建模的技术。提升建模旨在基于个体来确定治疗的有效性。我们假设电话辅导对客户的有效性平均为 70%。然而，效果可能因客户而异；一些客户更容易接受，另一些则不太容易。有关提升建模的更多信息，请参见www.steveklosterman.com/uplift-modeling/。

一种相对简单的监控模型实施方法是查看模型预测的分布是否随着时间变化，与用于模型训练的人群相比有所变化。我们在图 7.2中绘制了测试集的预测概率直方图。如果预测概率的直方图形状发生了显著变化，这可能表明特征发生了变化，或者特征与响应之间的关系发生了变化，模型可能需要重新训练或重建。为了量化分布的变化，感兴趣的读者可以查阅统计资源，学习卡方拟合优度检验或 Kolmogorov-Smirnov 检验。如果模型预测的分布发生变化，或者根据所选择的阈值，预测违约账户的比例发生明显变化，也可能会显现出来。

本章及全书中呈现的所有其他模型评估指标也可以是监控模型在生产环境中表现的好方法：十分位数图和等间距图、校准、ROC AUC 等。

预测建模中的伦理问题

随着机器学习的应用扩展到大多数现代企业，模型是否做出公平预测的问题受到了更多关注。公平性可以通过模型是否同样擅长为不同受保护类群体的成员做出预测来评估，例如，不同性别群体。

在本书中，我们采用了将性别从模型特征中移除的方法。然而，其他特征可能有效地充当性别的代理，因此即使性别没有作为特征使用，模型也可能对不同性别群体产生偏见的结果。筛查这种偏见可能性的一个简单方法是检查模型中使用的特征是否与受保护类群体有特别高的关联性，例如，可以通过使用 t 检验来进行检查。如果是这样，最好将这些特征从模型中移除。

如何判断模型是否公平，以及如果不公平，应该如何处理，这个问题是当前研究的热点。我们鼓励你熟悉像 AI Fairness 360（aif360.mybluemix.net/）这样的努力，这些努力正在提供工具，以提高机器学习中的公平性。在开始与公平性相关的工作之前，理解客户对于公平的定义至关重要，因为不同国家的法律和客户组织的具体政策可能会导致公平的定义因地区而异。

总结

在本章中，你学习了几种分析技术，以提供对模型性能的洞察，例如按模型预测区间划分的违约率的十分位和等距图表，以及如何调查模型校准的质量。通过使用模型测试集得出这些洞察并计算诸如 ROC AUC 等指标是很好的，因为这旨在代表模型在真实世界中新数据上的表现。

我们还看到如何进行模型性能的财务分析。虽然我们将此部分留到了书的最后，但对与模型决策相关的成本和节省的理解应该从典型项目开始时就有所了解。这些使数据科学家能够朝着通过增加利润或节省成本的具体目标努力。对于二分类模型，这一过程的关键步骤是选择一个预测概率阈值，以此来宣布一个正向预测，从而最大化由模型指导决策所带来的利润或节省。

最后，我们考虑了与交付和监控模型相关的任务，包括建立对照组以监控模型性能并测试任何由模型输出指导的程序有效性的想法。对照组的结构和模型监控策略会因项目而异，因此你需要在每个新案例中确定合适的行动方案。为了进一步了解在现实世界中使用模型，你应继续研究如实验设计、可以用于训练和部署模型的云平台（例如 AWS）以及预测建模中的公平性问题等主题。

你现在已经完成了项目，并准备向客户交付你的研究成果。除了将训练好的模型保存到磁盘或你可能提供给客户的其他数据产品或服务外，你还可能希望创建一个演示文稿，通常是一个幻灯片展示，详细说明你的进展。此类演示文稿的内容通常包括问题陈述、数据探索和清理结果、你构建的不同模型的性能比较、模型解释（如 SHAP 值），以及展示你的工作价值的财务分析。在制作你的工作演示文稿时，通常更好的是通过图片讲述故事，而不是大量的文字。在整本书中，我们展示了许多可用于此目的的可视化技术，你应继续探索描绘数据和建模结果的方式。

始终确保询问客户他们可能希望在演示文稿中看到的具体内容，并确保回答他们所有的问题。当客户看到你能够以易于理解的方式为他们创造价值时，你就成功了。

附录

1. 数据探索与清理

活动 1.01：探索数据集中剩余的财务特征

解决方案：

在开始之前，设置好你的环境并按如下方式加载已清理的数据集：

import pandas as pd
import matplotlib.pyplot as plt #import plotting package
#render plotting automatically
%matplotlib inline
import matplotlib as mpl #additional plotting functionality
mpl.rcParams['figure.dpi'] = 400 #high resolution figures
mpl.rcParams['font.size'] = 4 #font size for figures
from scipy import stats
import numpy as np
df = pd.read_csv('../../Data/Chapter_1_cleaned_data.csv')

1. 为剩余的财务特征创建特征名称列表。

   这些可以分为两组，因此我们将像之前一样列出特征名称，以便一起分析。你可以使用以下代码来实现：

   bill_feats = ['BILL_AMT1', 'BILL_AMT2', 'BILL_AMT3', \
                 'BILL_AMT4', 'BILL_AMT5', 'BILL_AMT6']
   pay_amt_feats = ['PAY_AMT1', 'PAY_AMT2', 'PAY_AMT3', \
                    'PAY_AMT4', 'PAY_AMT5', 'PAY_AMT6']
   

2. 使用.describe()方法查看账单金额特征的统计摘要。反思你所看到的内容。这合理吗？

   使用以下代码查看摘要：

   df[bill_feats].describe()
   

   输出应该如下所示：

   

   图 1.47：过去 6 个月账单金额的统计描述

   我们看到平均每月账单大约是 40,000 到 50,000 新台币。建议你检查一下本地货币的汇率。例如，1 美元约等于 30 新台币。做一下换算，问问自己，这个月度支付是否合理？我们也应该向客户确认这一点，但看起来是合理的。

   我们还注意到有些账单金额为负。这似乎是合理的，因为可能是前一个月的账单超额支付了，或许是预期当前账单中会有某项购买。类似的情况会导致该账户余额为负，意味着该账户持有人有了一个信用额度。

3. 使用以下代码，按 2x3 的网格方式可视化账单金额特征的直方图：

   df[bill_feats].hist(bins=20, layout=(2,3))
   

   图表应该是这样的：

   

   图 1.48：账单金额的直方图

   图 1.48中的直方图从多个方面来看是有意义的。大多数账户的账单金额较小。随着账单金额的增加，账户的数量逐渐减少。看起来账单金额的分布在每个月之间大致相似，因此我们没有像处理支付状态特征时那样发现数据不一致问题。该特征似乎通过了我们的数据质量检查。现在，我们将继续分析最后一组特征。

4. 使用.describe()方法，通过以下代码获取支付金额特征的摘要：

   df[pay_amt_feats].describe()
   

   输出应如下所示：

   

   图 1.49：过去 6 个月账单支付金额的统计描述

   平均支付金额大约比我们在前面的活动中总结的平均账单金额低一个数量级（10 的幂）。这意味着“平均情况”是一个每月未还清全部余额的账户。从我们对 PAY_1 特征的探索来看，这很有意义，因为该特征中最常见的值是 0（账户至少支付了最低付款额，但没有支付全部余额）。没有负支付，这也似乎是合理的。

5. 绘制与账单金额特征类似的支付金额特征的直方图，但还要使用 xrot 关键字参数对 x 轴 标签进行旋转，以避免重叠。使用 xrot=<角度> 关键字参数按给定的角度（以度为单位）旋转 x 轴 标签，使用以下代码：

   df[pay_amt_feats].hist(layout=(2,3), xrot=30)
   

   在我们的案例中，我们发现 30 度的旋转效果很好。绘图应如下所示：

   

   py

   图 1.50：原始支付金额数据的直方图

   这张图的快速浏览表明，这不是一个非常有用的图形；大多数直方图只有一个区间的高度较为显著。这不是可视化这些数据的有效方式。看起来，月度支付金额主要集中在包含 0 的区间中。那么，实际上有多少项是 0 呢？

6. 使用布尔掩码来查看支付金额数据中有多少项恰好等于 0，使用以下代码：使用以下代码执行此操作：

   pay_zero_mask = df[pay_amt_feats] == 0
   pay_zero_mask.sum()
   ```py
   
   输出应如下所示：
   
   ![图 1.51：支付金额等于 0 的账单计数    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/77ea5f1a4fb04f8d8a28f03933bcd2ea~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=eR5MjFPgkMyfcqJem0DrzEYnrgY%3D)
   
   图 1.51：支付金额等于 0 的账单计数
   
   `pay_zero_mask` 是一个包含 `True` 和 `False` 值的 DataFrame，表示支付金额是否等于 0。第二行对该 DataFrame 进行列求和，将 `True` 视为 1，将 `False` 视为 0，因此列的和表示每个特征中支付金额为 0 的账户数。
   
   我们可以看到，约 20-25% 的账户在任何给定的月份里账单支付额为 0。然而，大多数账单支付额大于 0。那么，为什么我们在直方图中看不到它们呢？这是由于账单支付额的**范围**相对于大多数账单支付额的值。
   
   在统计摘要中，我们可以看到一个月中的最大账单支付额通常比平均账单支付额大两个数量级（100 倍）。看起来这些非常大的账单支付可能只有少数几个。然而，由于直方图的创建方式，使用相同大小的区间，几乎所有数据都被聚集在最小的区间中，较大的区间几乎不可见，因为它们的账户数太少。我们需要一种有效的策略来可视化这些数据。
   
   

7. 忽略前一步创建的掩码中的 0 支付值，使用 pandas 的.apply()方法和 NumPy 的np.log10()方法绘制非零支付的对数转换直方图。你可以使用.apply()将任何函数（包括log10）应用到 DataFrame 的所有元素。使用以下代码来完成此操作：

   df[pay_amt_feats][~pay_zero_mask].apply(np.log10)\
                                    .hist(layout=(2,3))
   ```py
   
   这是 pandas 的一个相对高级的使用方法，所以如果你自己没弄明白也不必担心。然而，开始理解如何用相对少量的代码在 pandas 中做很多事情是很有帮助的。
   
   输出应如下所示：
   
   ![图 1.52：非零账单支付金额的 10 为基对数    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/35510d8a60304f7b877356875bb8f2d4~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=RNPAg9JszlTgPG97iK1Is2xjXQM%3D)
   

图 1.52：非零账单支付金额的 10 为基对数

虽然我们本可以尝试创建不同宽度的区间来更好地可视化支付金额，但另一种常用且便捷的方法是对数变换，或称为对数变换。我们使用了 10 为基的对数变换。大致来说，这种变换告诉我们一个数值中有多少个零。换句话说，一个余额至少为 100 万美元但不到 1000 万美元的账户，其对数变换结果会是 6 到 7 之间，因为 106 = 1,000,000（而log10(1,000,000) = 6），而 107 = 10,000,000。

为了将此变换应用到我们的数据，首先需要屏蔽掉零支付值，因为log10(0)是未定义的（另一种常见方法是对所有值加上一个非常小的数字，例如 0.01，这样就没有零值）。我们使用了 Python 逻辑运算符not（~）和我们已经创建的零值掩码。然后，我们使用了 pandas 的.apply()方法，它可以将我们喜欢的任何函数应用到我们选择的数据上。在这种情况下，我们希望应用的是一个基于 10 的对数，使用np.log10来计算。最后，我们对这些值绘制了直方图。

结果是更加有效的数据可视化：这些值在直方图的区间中分布得更具信息量。我们可以看到，最常出现的账单支付金额位于千元范围内（log10(1,000) = 3），这与我们在统计摘要中观察到的平均账单支付金额一致。也有一些非常小的账单支付金额，以及少数较大的支付金额。总体来看，账单支付金额的分布从每月来看似乎非常一致，因此我们没有发现该数据中存在任何潜在问题。

2. Scikit-Learn 简介与模型评估

活动 2.01：使用新特征执行逻辑回归并创建精准率-召回率曲线

解决方案：

1. 使用 scikit-learn 的train_test_split生成新的训练和测试数据集。这次，使用LIMIT_BAL，即账户的信用额度，作为特征，而不是EDUCATION。

   执行以下代码来实现：

   X_train_2, X_test_2, y_train_2, y_test_2 = train_test_split\
                                             (df['LIMIT_BAL']\
                                              .values\
                                              .reshape(-1,1),\
                                              df['default'\
                                                 'payment next'\
                                                 'month'].values,\
                                              test_size=0.2,\
                                              random_state=24))
   ```py
   
   请注意，这里我们创建了新的训练和测试数据集，并且变量名称也发生了变化。
   
   

2. 使用你拆分后的训练数据训练一个逻辑回归模型。

   以下代码实现了这个功能：

   example_lr.fit(X_train_2, y_train_2)
   ```py
   
   如果你在一个单一的笔记本中运行整个章节，可以重新使用之前使用的模型对象`example_lr`。你可以**重新训练**这个对象，以学习这个新特征与响应之间的关系。如果你愿意，也可以尝试不同的训练/测试拆分，而不必创建新的模型对象。在这些场景中，现有的模型对象已经被**原地更新**。
   
   

3. 创建测试数据的预测概率数组。

   这里是此步骤的代码：

   y_test_2_pred_proba = example_lr.predict_proba(X_test_2)
   ```py
   
   

4. 使用预测的概率和测试数据的真实标签计算 ROC AUC。将其与使用EDUCATION特征的 ROC AUC 进行比较。

   运行以下代码进行这一步操作：

   metrics.roc_auc_score(y_test_2, y_test_2_pred_proba[:,1])
   ```py
   
   输出结果如下：
   
   

   0.6201990844642832

   
   请注意，我们对预测的概率数组进行了索引，以便从第二列获取正类的预测概率。与`EDUCATION`的逻辑回归 ROC AUC 相比，结果如何？AUC 更高。这可能是因为现在我们使用的是与账户财务状况（信用额度）相关的特征，来预测与账户财务状况相关的另一项内容（是否违约），而不是使用与财务关系较弱的特征。
   
   

5. 绘制 ROC 曲线。

   这里是实现此功能的代码；它与我们在前一个练习中使用的代码类似：

   fpr_2, tpr_2, thresholds_2 = metrics.roc_curve\
                                (y_test_2, \
                                 y_test_2_pred_proba[:,1])
   plt.plot(fpr_2, tpr_2, '*-')
   plt.plot([0, 1], [0, 1], 'r--')
   plt.legend(['Logistic regression', 'Random chance'])
   plt.xlabel('FPR')
   plt.ylabel('TPR')
   plt.title('ROC curve for logistic regression with '\
             'LIMIT_BAL feature')
   ```py
   
   图形应如下所示：
   
   ![图 2.30：LIMIT_BAL 逻辑回归的 ROC 曲线    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/b017ce5967914d34938152483ffa3ada~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=%2F0lhBVQdPfhWDd2L7GaylWH8BLA%3D)
   
   图 2.30：LIMIT_BAL 逻辑回归的 ROC 曲线
   
   这看起来有点像我们希望看到的 ROC 曲线：它比仅使用`EDUCATION`特征的模型更远离随机机会线。还注意到，真实和假阳性率的变化在阈值范围内更加平滑，反映了`LIMIT_BAL`特征具有更多不同的值。
   
   

6. 使用 scikit-learn 的功能计算测试数据上精确度-召回率曲线的数据。

   精确度通常与召回率一起考虑。我们可以使用sklearn.metrics中的precision_recall_curve来自动调整阈值，并计算每个阈值下的精确度和召回率对。以下是提取这些值的代码，类似于roc_curve：

   precision, recall, thresh_3 = metrics.precision_recall_curve\
                                 (y_test_2,\
                                  y_test_2_pred_proba[:,1])
   ```py
   
   

7. 使用 matplotlib 绘制精确度-召回率曲线：我们可以通过以下代码实现这一点。

   注意，我们将召回率放在x轴，将精确度放在y轴，并将坐标轴的限制设置为[0, 1]范围：

   plt.plot(recall, precision, '-x')
   plt.xlabel('Recall')
   plt.ylabel('Precision')
   plt.title('Precision and recall for the logistic'\
             'regression 'with LIMIT_BAL')
   plt.xlim([0, 1])
   plt.ylim([0, 1])
   ```py
   
   ![图 2.31：精确度-召回率曲线图    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/683d99ba90a943db938e50d77ec018e8~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=%2F3vza2sIQBMJdLB7sB44JosTWak%3D)
   
   图 2.31：精确度-召回率曲线图
   
   

8. 使用 scikit-learn 计算精确度-召回率曲线下的面积。

   以下是实现此操作的代码：

   metrics.auc(recall, precision)
   ```py
   
   你将获得以下输出：
   
   

   0.31566964427378624

   
   我们看到，精确率-召回率曲线表明，该模型的精确率通常相对较低；在几乎所有阈值范围内，精确率（即正确的正类分类所占比例）都不到一半。我们可以通过计算精确率-召回率曲线下面积来比较这个分类器与我们可能考虑的其他模型或特征集。
   
   Scikit-learn 提供了一个计算任何`x-y`数据 AUC 的功能，使用的是梯形规则，你可能还记得这个方法来自微积分：`metrics.auc`。我们使用这个功能来获取精确率-召回率曲线下面积。
   
   

9. 现在重新计算 ROC AUC，不过这次要计算训练数据的 ROC AUC。这在概念上和定量上与之前的计算有何不同？

   首先，我们需要使用训练数据而不是测试数据来计算预测概率。然后，我们可以使用训练数据标签来计算 ROC AUC。以下是代码：

   y_train_2_pred_proba = example_lr.predict_proba(X_train_2)
   metrics.roc_auc_score(y_train_2, y_train_2_pred_proba[:,1])
   ```py
   
   你应该得到以下输出：
   
   

   0.6182918113358344

定量来看，我们可以看到这个 AUC 与我们之前计算的测试数据 ROC AUC 差别不大。两者都大约是 0.62。概念上，这有什么不同？当我们在训练数据上计算这个指标时，我们衡量的是模型在预测“教会”模型如何进行预测的相同数据上的能力。我们看到的是模型如何拟合数据。另一方面，测试数据的指标表示模型在未见过的外部样本数据上的表现。如果这些得分差异很大，通常表现为训练得分高于测试得分，那就意味着虽然模型很好地拟合了数据，但训练好的模型无法很好地泛化到新的、未见过的数据。

在这种情况下，训练得分和测试得分相似，这意味着模型在训练数据和未见过的数据（外部样本数据）上的表现差不多。我们将在第四章，偏差-方差权衡中学习更多关于通过比较训练和测试得分可以获得的见解。

3. 逻辑回归和特征探索的详细信息

活动 3.01：拟合逻辑回归模型并直接使用系数

解答：

前几个步骤与我们在之前的活动中做的类似：

1. 创建一个训练/测试数据集（80/20），以PAY_1和LIMIT_BAL作为特征：

   from sklearn.model_selection import train_test_split
   X_train, X_test, y_train, y_test = train_test_split(
       df[['PAY_1', 'LIMIT_BAL']].values,
       df['default payment next month'].values,
       test_size=0.2, random_state=24)
   ```py
   
   

2. 导入LogisticRegression，使用默认选项，但将求解器设置为'liblinear'：

   from sklearn.linear_model import LogisticRegression
   lr_model = LogisticRegression(solver='liblinear')
   ```py
   
   

3. 在训练数据上进行训练，并使用测试数据获取预测类别以及类别概率：

   lr_model.fit(X_train, y_train)
   y_pred = lr_model.predict(X_test)
   y_pred_proba = lr_model.predict_proba(X_test)
   ```py
   
   

4. 从训练好的模型中提取系数和截距，并手动计算预测概率。你需要在特征中添加一列值为 1 的列，以便与截距相乘。

   首先，让我们创建特征数组，添加一列 1 值，使用水平堆叠：

   ones_and_features = np.hstack\
                       ([np.ones((X_test.shape[0],1)), X_test])
   ```py
   
   现在我们需要截距和系数，我们将从 scikit-learn 输出中重新调整形状并连接它们：
   
   

   intercept_and_coefs = np.concatenate
   ([lr_model.intercept_.reshape(1,1),
   lr_model.coef_], axis=1)

   
   为了反复将截距和系数乘以 `ones_and_features` 的所有行，并求每行的和（也就是求线性组合），你可以使用乘法和加法将这些全部写出来。不过，使用点积会更快：
   
   

   X_lin_comb = np.dot(intercept_and_coefs,
   np.transpose(ones_and_features))

   
   现在 `X_lin_comb` 包含了我们需要传递给我们定义的 sigmoid 函数的参数，以计算预测概率：
   
   

   y_pred_proba_manual = sigmoid(X_lin_comb)

5. 使用 0.5 的阈值，手动计算预测的类别。与 scikit-learn 输出的类别预测进行比较。

   手动预测的概率 y_pred_proba_manual 应该与 y_pred_proba 相同，我们马上检查这一点。首先，使用阈值手动预测类别：

   y_pred_manual = y_pred_proba_manual >= 0.5
   ```py
   
   这个数组的形状将与 `y_pred` 不同，但它应该包含相同的值。我们可以像这样检查两个数组的所有元素是否相等：
   
   

   np.array_equal(y_pred.reshape(1,-1), y_pred_manual)

   
   如果数组相等，这应该返回一个逻辑值`True`。
   
   

6. 使用 scikit-learn 的预测概率和手动预测的概率计算 ROC AUC，并进行比较。

   首先，导入以下内容：

   from sklearn.metrics import roc_auc_score
   ```py
   
   然后，在两个版本上计算此指标，确保访问正确的列，或根据需要调整形状：
   
   ![图 3.37：从预测概率计算 ROC AUC    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/9ce6c9b7bbfb423897a054fba4e9aee5~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=Cv%2FkryaSjY1X8Txvj4Sl7NQdP7M%3D)
   

图 3.37：从预测概率计算 ROC AUC

实际上，AUC 是相同的。我们在这里做了什么？我们已经确认，从这个拟合的 scikit-learn 模型中，我们实际上只需要三个数字：截距和两个系数。一旦我们得到了这些，就可以使用几行代码，借助数学函数，来创建模型预测，这与直接从 scikit-learn 生成的预测是等效的。

这有助于确认你对知识的理解，但除此之外，你为什么要这么做呢？我们将在最后一章讨论模型部署。不过，根据你的具体情况，你可能会遇到一种情形，在其中你没有 Python 环境来为模型输入新的特征进行预测。例如，你可能需要完全在 SQL 中进行预测。虽然这在一般情况下是一个限制，但使用逻辑回归时，你可以利用 SQL 中可用的数学函数重新创建逻辑回归预测，只需要将截距和系数粘贴到 SQL 代码的某个地方即可。点积可能不可用，但你可以使用乘法和加法来实现相同的目的。

那么，结果如何呢？我们看到的是，通过稍微提高模型的表现，我们可以超过之前的尝试：在上一章的活动中，仅使用 LIMIT_BAL 作为特征时，ROC AUC 稍低为 0.62，而此处为 0.63。下一章我们将学习使用逻辑回归的高级技术，进一步提升性能。

4. 偏差-方差权衡

活动 4.01：使用案例研究数据进行交叉验证和特征工程

解答：

1. 从案例研究数据的 DataFrame 中选择特征。

   你可以使用我们在本章中已经创建的特征名称列表，但请确保不要包括响应变量，它本可以是一个非常好的（但完全不合适的）特征：

   features = features_response[:-1]
   X = df[features].values
   ```py
   
   

2. 使用随机种子 24 进行训练/测试集拆分：

   X_train, X_test, y_train, y_test = \
   train_test_split(X, df['default payment next month'].values,
                    test_size=0.2, random_state=24)
   ```py
   
   我们将在后续中使用这个数据，并将其作为未见测试集进行保留。通过指定随机种子，我们可以轻松创建包含其他建模方法的独立笔记本，使用相同的训练数据。
   
   

3. 实例化 MinMaxScaler 来缩放数据，如下所示的代码：

   from sklearn.preprocessing import MinMaxScaler
   min_max_sc = MinMaxScaler()
   ```py
   
   

4. 实例化一个逻辑回归模型，使用 saga 求解器，L1 惩罚，并将 max_iter 设置为 1000，因为我们希望允许求解器有足够的迭代次数来找到一个好的解：

   lr = LogisticRegression(solver='saga', penalty='l1',
                           max_iter=1000)
   ```py
   
   

5. 导入 Pipeline 类并创建一个包含缩放器和逻辑回归模型的管道，分别使用 'scaler' 和 'model' 作为步骤的名称：

   from sklearn.pipeline import Pipeline
   scale_lr_pipeline = Pipeline(
       steps=[('scaler', min_max_sc), ('model', lr)])
   ```py
   
   

6. 使用 get_params 和 set_params 方法查看每个阶段的参数，并更改它们（在你的笔记本中分别执行以下每行代码并观察输出）：

   scale_lr_pipeline.get_params()
   scale_lr_pipeline.get_params()['model__C']
   scale_lr_pipeline.set_params(model__C = 2)
   ```py
   
   

7. 创建一个较小范围的 C 值进行交叉验证测试，因为这些模型相比我们之前的练习，在训练和测试时需要更多的数据，因此训练时间会更长；我们推荐的 C = [102*, 10, 1, 10*-1*, 10*-2*, 10*-3*]*：

   C_val_exponents = np.linspace(2,-3,6)
   C_vals = np.float(10)**C_val_exponents
   ```py
   
   

8. 创建 cross_val_C_search 函数的新版本，命名为 cross_val_C_search_pipe。这个函数将接受一个 pipeline 参数，而不是 model 参数。函数内的更改是使用 set_params(model__C = <你想测试的值>) 设置 C 值，替换模型为管道，并在 fit 和 predict_proba 方法中使用管道，且通过 pipeline.get_params()['model__C'] 获取 C 值以打印状态更新。

   更改如下：

   def cross_val_C_search_pipe(k_folds, C_vals, pipeline, X, Y):
   ##[…]
   pipeline.set_params(model__C = C_vals[c_val_counter])
   ##[…]
   pipeline.fit(X_cv_train, y_cv_train)
   ##[…]
   y_cv_train_predict_proba = pipeline.predict_proba(X_cv_train)
   ##[…]
   y_cv_test_predict_proba = pipeline.predict_proba(X_cv_test)
   ##[…]
   print('Done with C = {}'.format(pipeline.get_params()\
                                   ['model__C']))
   ```py
   
   注意
   
   有关完整的代码，请参阅 [`packt.link/AsQmK`](https://packt.link/AsQmK)。
   
   

9. 按照上一个练习的方式运行这个函数，但使用新的C值范围、你创建的管道以及从案例研究数据训练集分割得到的特征和响应变量。你可能会在此或后续步骤中看到关于求解器未收敛的警告；你可以尝试调整tol或max_iter选项以尝试实现收敛，尽管使用max_iter = 1000时得到的结果应该是足够的。以下是执行此操作的代码：

   cv_train_roc_auc, cv_test_roc_auc, cv_test_roc = \
   cross_val_C_search_pipe(k_folds, C_vals, scale_lr_pipeline,
                           X_train, y_train)
   ```py
   
   你将获得以下输出：
   
   

   Done with C = 100.0 Done with C = 10.0 Done with C = 1.0 Done with C = 0.1 Done with C = 0.01 Done with C = 0.001

10. 使用以下代码绘制每个C值下，跨折叠的平均训练和测试 ROC AUC：

    plt.plot(C_val_exponents, np.mean(cv_train_roc_auc, axis=0),
             '-o', label='Average training score')
    plt.plot(C_val_exponents, np.mean(cv_test_roc_auc, axis=0),
             '-x', label='Average testing score')
    plt.ylabel('ROC AUC')
    plt.xlabel('log$_{10}$(C)')
    plt.legend()
    plt.title('Cross-validation on Case Study problem')
    ```py
    
    你将获得以下输出：
    
    ![图 4.25：交叉验证测试性能    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/af5cd91267fc43ec9de62316089ca87e~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=yEY5jhNr%2BTlAGPEhmw9CA7LIUlg%3D)
    
    图 4.25：交叉验证测试性能
    
    你应该注意到，正则化在这里并没有带来太多好处，正如预期的那样：对于较低的*C*值，即较强的正则化，模型的测试（以及训练）性能下降。虽然通过使用所有可用特征，我们能够提高模型的性能，但似乎并没有出现过拟合。相反，训练和测试分数差不多。与其说是过拟合，不如说我们可能存在欠拟合的情况。让我们尝试构建一些交互特征，看看它们是否能提升性能。
    
    

11. 使用以下代码为案例研究数据创建交互特征，并确认新特征的数量是合理的：

    from sklearn.preprocessing import PolynomialFeatures
    make_interactions = PolynomialFeatures(degree=2,
                                           interaction_only=True,
                                           include_bias=False)
    X_interact = make_interactions.fit_transform(X)
    X_train, X_test, y_train, y_test = train_test_split(
        X_interact, df['default payment next month'].values,
        test_size=0.2, random_state=24)
    print(X_train.shape)
    print(X_test.shape)
    ```py
    
    你将获得以下输出：
    
    

    (21331, 153) (5333, 153)

    
    从中你应该看到新特征的数量是 153，计算方式是*17 + "17 选 2" = 17 + 136 = 153*。*“17 选 2”*部分来自于从 17 个原始特征中选择所有可能的两个特征的组合进行交互。
    
    

12. 重复交叉验证过程，并观察使用交互特征时模型的表现；也就是说，重复步骤 9和步骤 10。请注意，由于特征数量增加，这将需要更多时间，但应该不会超过 10 分钟。

    你将获得以下输出：

    

图 4.26：通过添加交互特征改善的交叉验证测试性能

那么，交互特征是否改善了平均交叉验证测试性能？正则化是否有用？

构造交互特征使得最佳模型测试分数在各折叠之间平均约为ROC AUC = 0.74，而没有交互特征时大约是 0.72。这些分数出现在C = 100时，即几乎没有正则化。在包含交互特征的模型的训练与测试分数对比图上，你可以看到训练分数稍高于测试分数，因此可以说存在一定程度的过拟合。然而，在这里我们无法通过正则化来提高测试分数，因此这可能不是过拟合的一个问题实例。在大多数情况下，任何能够产生最高测试分数的策略就是最佳策略。

总结来说，添加交互特征提高了交叉验证的表现，而正则化在目前使用逻辑回归模型的案例中似乎并不有用。我们将在稍后进行全量训练数据拟合的步骤，当我们尝试其他模型并在交叉验证中找到最佳模型时再进行。

5. 决策树与随机森林

活动 5.01：使用随机森林进行交叉验证网格搜索

解决方案：

1. 创建一个字典，表示将要搜索的max_depth和n_estimators超参数的网格。包括深度为 3、6、9 和 12，以及树的数量为 10、50、100 和 200。保持其他超参数为默认值。使用以下代码创建字典：

   rf_params = {'max_depth':[3, 6, 9, 12],
                'n_estimators':[10, 50, 100, 200]}
   ```py
   
   注意
   
   还有许多其他可能的超参数需要搜索。特别是，scikit-learn 文档中关于随机森林指出：“使用这些方法时，主要调整的参数是`n_estimators`和`max_features`”，并且“经验上良好的默认值是……分类任务时，`max_features=sqrt(n_features)`。”
   
   来源：https://scikit-learn.org/stable/modules/ensemble.html#parameters
   
   出于本书的目的，我们将使用`max_features='auto'`（等同于`sqrt(n_features)`）并将探索限制在`max_depth`和`n_estimators`上，以缩短运行时间。在实际情况中，你应根据可承受的计算时间探索其他超参数。记住，为了在特别大的参数空间中搜索，你可以使用`RandomizedSearchCV`，以避免为网格中每个超参数组合计算所有指标。
   
   

2. 使用我们在本章之前使用的相同选项实例化一个GridSearchCV对象，但使用步骤 1 中创建的超参数字典。设置verbose=2以查看每次拟合的输出。你可以重用我们一直在使用的相同随机森林模型对象rf，或者创建一个新的。创建一个新的随机森林对象，并使用以下代码实例化GridSearchCV类：

   rf = RandomForestClassifier(n_estimators=10,\
                               criterion='gini',\
                               max_depth=3,\
                               min_samples_split=2,\
                               min_samples_leaf=1,\
                               min_weight_fraction_leaf=0.0,\
                               max_features='auto',\
                               max_leaf_nodes=None,\
                               min_impurity_decrease=0.0,\
                               min_impurity_split=None,\
                               bootstrap=True,\
                               oob_score=False,\
                               n_jobs=None,
                               random_state=4,\
                               verbose=0,\
                               warm_start=False,\
                               class_weight=None)
   cv_rf = GridSearchCV(rf, param_grid=rf_params,\
                        scoring='roc_auc',\
                        n_jobs=-1,\
                        refit=True,\
                        cv=4,\
                        verbose=2,\
                        error_score=np.nan,\
                        return_train_score=True)
   ```py
   
   

3. 在训练数据上拟合GridSearchCV对象。使用以下代码执行网格搜索：

   cv_rf.fit(X_train, y_train)
   ```py
   
   由于我们选择了`verbose=2`选项，你将看到笔记本中相对较多的输出。每个超参数组合都会有输出，并且对于每个折叠，都有拟合和测试的输出。以下是输出的前几行：
   
   ![图 5.22：交叉验证的冗长输出    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/60610ca0048c46e98a1a05e80eec73e5~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=25OHvC7JH6VXb6KdBIroE3rAsp4%3D)
   
   图 5.22：交叉验证的冗长输出
   
   虽然在较短的交叉验证过程中不必查看所有这些输出，但对于较长的交叉验证，看到这些输出可以让你确认交叉验证正在进行，并且能给你一些关于不同超参数组合下拟合所需时间的概念。如果某些操作花费的时间过长，你可能需要通过点击笔记本顶部的停止按钮（方形）中断内核，并选择一些运行时间更短的超参数，或者使用一个更有限的超参数集。
   
   当这一切完成时，你应该会看到以下输出：
   
   ![图 5.22：交叉验证完成后的输出    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/c7652790593c4702806096401a8a0787~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=3Cnv%2FIBv29kf1GjAEWYw9arXyXs%3D)
   
   图 5.23：交叉验证完成后的输出
   
   这个交叉验证任务大约运行了 2 分钟。随着任务规模的增大，你可能希望通过`n_jobs`参数探索并行处理，以查看是否能够加速搜索。使用`n_jobs=-1`进行并行处理时，运行时间应该比串行处理更短。然而，在并行处理中，你将无法看到每个单独模型拟合操作的输出，如*图 5.23*所示。
   
   

4. 将网格搜索结果放入 pandas DataFrame 中。使用以下代码将结果放入 DataFrame：

   cv_rf_results_df = pd.DataFrame(cv_rf.cv_results_)
   ```py
   
   

5. 创建一个pcolormesh可视化图，显示每个超参数组合的平均测试分数。以下是创建交叉验证结果网格图的代码。它与我们之前创建的示例图类似，但带有特定于此处执行的交叉验证的注释：

   ax_rf = plt.axes()
   pcolor_graph = ax_rf.pcolormesh\
                  (xx_rf, yy_rf,\
                   cv_rf_results_df['mean_test_score']\
                   .values.reshape((4,4)), cmap=cm_rf)
   plt.colorbar(pcolor_graph, label='Average testing ROC AUC')
   ax_rf.set_aspect('equal')
   ax_rf.set_xticks([0.5, 1.5, 2.5, 3.5])
   ax_rf.set_yticks([0.5, 1.5, 2.5, 3.5])
   ax_rf.set_xticklabels\
   ([str(tick_label) for tick_label in rf_params['n_estimators']])
   ax_rf.set_yticklabels\
   ([str(tick_label) for tick_label in rf_params['max_depth']])
   ax_rf.set_xlabel('Number of trees')
   ax_rf.set_ylabel('Maximum depth')
   ```py
   
   我们与之前示例的主要区别在于，我们不是绘制从 1 到 16 的整数，而是绘制我们通过`cv_rf_results_df['mean_test_score'].values.reshape((4,4))`检索并重塑的平均测试分数。这里的另一个新内容是，我们使用列表推导式基于网格中超参数的数值，创建字符串的列表作为刻度标签。我们从定义的字典中访问这些数值，然后在列表推导式中将它们逐一转换为`str`（字符串）数据类型，例如，`ax_rf.set_xticklabels([str(tick_label) for tick_label in rf_params['n_estimators']])`。我们已经使用`set_xticks`将刻度位置设置为我们希望显示刻度的位置。此外，我们使用`ax_rf.set_aspect('equal')`创建了一个正方形图形。图形应该如下所示：
   
   ![图 5.24：对随机森林进行交叉验证的结果    在包含两个超参数的网格上    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/81c6f0a92b9e4f68b372a090e3217b06~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=RucoGnH31s9KhPDvfLd8tgnrfDw%3D)
   
   图 5.24：在两个超参数的网格上进行随机森林交叉验证的结果
   
   

6. 确定使用哪一组超参数。

   从我们的网格搜索中可以得出什么结论？显然，使用深度大于 3 的树有一定的优势。在我们尝试的参数组合中，max_depth=9且树木数量为 200 的组合提供了最佳的平均测试得分，你可以在数据框中查找并确认其 ROC AUC = 0.776。

   这是我们迄今为止所有努力中找到的最佳模型。

   在实际场景中，我们可能会进行更彻底的搜索。接下来的好步骤是尝试更多的树木数量，并且不再花费时间在n_estimators < 200 上，因为我们知道至少需要 200 棵树才能获得最佳性能。你还可以更精细地搜索max_depth的空间，而不是像我们这里一样按 3 递增，并尝试其他超参数，比如max_features。不过为了我们的目的，我们将假设已经找到了最佳超参数并继续前进。

6. 梯度提升、XGBoost 和 SHAP 值

活动 6.01：使用 XGBoost 建模案例研究数据并通过 SHAP 解释模型

解决方案：

在本次活动中，我们将结合本章学习的内容，使用合成数据集并将其应用到案例研究数据中。我们将观察 XGBoost 模型在验证集上的表现，并使用 SHAP 值解释模型预测。我们已经通过替换之前忽略的PAY_1特征缺失值的样本来准备该数据集，同时保留了没有缺失值的样本的训练/测试分割。你可以在本活动的笔记本附录中查看数据是如何准备的。

1. 加载为本次练习准备的案例研究数据。文件路径为../../Data/Activity_6_01_data.pkl，变量包括：features_response, X_train_all, y_train_all, X_test_all, y_test_all：

   with open('../../Data/Activity_6_01_data.pkl', 'rb') as f:
       features_response, X_train_all, y_train_all, X_test_all,\
       y_test_all = pickle.load(f)
   ```py
   
   

2. 定义一个验证集，用于训练 XGBoost 并进行提前停止：

   from sklearn.model_selection import train_test_split
   X_train_2, X_val_2, y_train_2, y_val_2 = \
   train_test_split(X_train_all, y_train_all,\
                    test_size=0.2, random_state=24)
   ```py
   
   

3. 实例化一个 XGBoost 模型。我们将使用lossguide增长策略，并检查不同max_leaves值下的验证集表现：

   xgb_model_4 = xgb.XGBClassifier(
       n_estimators=1000,
       max_depth=0,
       learning_rate=0.1,
       verbosity=1,
       objective='binary:logistic',
       use_label_encoder=False,
       n_jobs=-1,
       tree_method='hist',
       grow_policy='lossguide')
   ```py
   
   

4. 搜索max_leaves的值，从 5 到 200，步长为 5：

   max_leaves_values = list(range(5,205,5))
   ```py
   
   

5. 创建用于提前停止的评估集：

   eval_set_2 = [(X_train_2, y_train_2), (X_val_2, y_val_2)]
   ```py
   
   

6. 遍历超参数值，并创建一个验证集的 ROC AUC 列表，采用与练习 6.01：随机化网格搜索调优 XGBoost 超参数相同的方法：

   %%time
   val_aucs = []
   for max_leaves in max_leaves_values:
       #Set parameter and fit model
       xgb_model_4.set_params(**{'max_leaves':max_leaves})
       xgb_model_4.fit(X_train_2, y_train_2,\
                       eval_set=eval_set_2,\
                       eval_metric='auc',\
                       verbose=False,\
                       early_stopping_rounds=30)
       #Get validation score
       val_set_pred_proba = xgb_model_4.predict_proba(X_val_2)[:,1]
       val_aucs.append(roc_auc_score(y_val_2, val_set_pred_proba))
   ```py
   
   

7. 创建一个数据框，记录超参数搜索结果，并绘制验证 AUC 与max_leaves的关系：

   max_leaves_df_2 = \
   pd.DataFrame({'Max leaves':max_leaves_values,\
                 'Validation AUC':val_aucs})
   mpl.rcParams['figure.dpi'] = 400
   max_leaves_df_2.set_index('Max leaves').plot()
   ```py
   
   图表应该看起来像这样：
   
   ![图 6.15：案例研究数据中验证 AUC 与 max_leaves 的关系    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/2e79319787834b98a14b1e10950cc632~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=cX%2FFIRW7zqYZuQCX92gAvYduNoE%3D)
   
   图 6.15：案例研究数据中验证 AUC 与`max_leaves`的关系
   
   尽管关系有些噪声，我们可以看到，通常较低的`max_leaves`值会导致较高的验证集 ROC AUC。这是因为通过允许较少的叶子来限制树的复杂度，从而减少了过拟合，增加了验证集得分。
   
   

8. 观察max_leaves对应于验证集上最高 ROC AUC 的数量：

   max_auc_2 = max_leaves_df_2['Validation AUC'].max()
   max_auc_2
   max_ix_2 = max_leaves_df_2['Validation AUC'] == max_auc_2
   max_leaves_df_2[max_ix_2]
   ```py
   
   结果应如下所示：
   
   ![图 6.16：案例研究数据的最佳 max_leaves 和验证集 AUC    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/a02b5ff25b7449f4b3e5014febdfd837~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=9Xml8C1AjziTKMlBLl5FATmVgd8%3D)
   
   图 6.16：案例研究数据的最佳`max_leaves`和验证集 AUC
   
   我们希望结合之前在建模案例数据方面的努力来解释这些结果。这不是一个完美的对比，因为这里我们在训练和验证数据中有缺失值，而之前我们忽略了这些缺失值，并且这里只有一个验证集，而不像之前使用的 k 折交叉验证（尽管有兴趣的读者可以尝试在 XGBoost 中使用 k 折交叉验证进行多个训练/验证分割，并加上早停）。
   
   然而，即使考虑到这些限制，下面的验证结果应提供一个类似于我们之前进行的 k 折交叉验证的样本外表现度量。我们注意到，这里的验证 ROC AUC 值为 0.779，比之前在*活动 5.01*中使用随机森林获得的 0.776 略高，*随机森林的交叉验证网格搜索*，出自*第五章，决策树与随机森林*。这些验证分数相当接近，实际上使用这两种模型都应该是可以的。接下来我们将继续使用 XGBoost 模型。
   
   

9. 使用最佳超参数重新拟合 XGBoost 模型：

   xgb_model_4.set_params(**{'max_leaves':40})
   xgb_model_4.fit(X_train_2, y_train_2, eval_set=eval_set_2,
                   eval_metric='auc',
                   verbose=False, early_stopping_rounds=30)
   ```py
   
   

10. 为了能够检查验证集的 SHAP 值，我们需要将此数据做成一个数据框：

    X_val_2_df = pd.DataFrame(data=X_val_2,
                              columns=features_response[:-1])
    ```py
    
    

11. 创建一个 SHAP 解释器，使用验证数据作为背景数据集，获取 SHAP 值，并制作一个摘要图：

    explainer_2 = shap.explainers.Tree(xgb_model_4, data=X_val_2_df)
    shap_values_2 = explainer_2(X_val_2_df)
    mpl.rcParams['figure.dpi'] = 75
    shap.summary_plot(shap_values_2.values, X_val_2_df)
    ```py
    
    图应如下所示：
    
    ![图 6.17：案例研究数据在验证集上的 XGBoost 模型 SHAP 值    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/a211539857c34074afaf2f8d2534937c~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=vPArQDc7LCMsUSmskwD6T4qMmXM%3D)
    
    图 6.17：案例研究数据在验证集上的 XGBoost 模型 SHAP 值
    
    从*图 6.17*中，我们可以看到，XGBoost 模型中最重要的特征与我们在*第五章，决策树与随机森林*中探索的随机森林模型中的特征略有不同（参见*图 5.15*）。`PAY_1`不再是最重要的特征，尽管它仍然在第三位，仍然很重要。现在，`LIMIT_BAL`（借款人的信用额度）是最重要的特征。这个结果是合理的，因为贷方可能根据借款人的风险来设定信用额度，因此它应该是一个很好的违约风险预测因子。
    
    让我们探索`LIMIT_BAL`是否与其他特征有任何有趣的 SHAP 交互。我们可以通过不为颜色参数指定特征，允许`shap`包自动选择与其他特征交互最强的特征，从而绘制散点图。
    
    

12. 绘制LIMIT_BAL的 SHAP 值散点图，按最强交互特征进行着色：

    shap.plots.scatter(shap_values_2[:,'LIMIT_BAL'],
                       color=shap_values_2)
    ```py
    
    该图应如下所示：
    
    ![图 6.18：LIMIT_BAL 和其他特征的 SHAP 值散点图    与特征最强交互的特征    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/59247609abc043cc9dfd58534872a790~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=BryVekpAwODTKXsbDp52QgQguGQ%3D)
    
    图 6.18：`LIMIT_BAL`和与之有最强交互特征的 SHAP 值散点图
    
    `BILL_AMT2`，即两个月前的账单金额，和`LIMIT_BAL`的交互最强。我们可以看到，对于大多数`LIMIT_BAL`的值，如果账单特别高，这会导致更高的 SHAP 值，意味着违约风险增加。通过观察*图 6.18*中最红的点分布，可以发现这些点集中在点带的顶部。这符合直觉：即使借款人获得了较大的信用额度，如果账单金额过大，也可能意味着违约风险增加。
    
    最后，我们将保存模型以及训练和测试数据以供分析，并交付给我们的商业伙伴。我们通过使用 Python 的`pickle`功能来完成此操作。
    
    

13. 将训练好的模型以及训练和测试数据保存到文件中：

    with open('../Data/xgb_model_w_data.pkl', 'wb') as f:
        pickle.dump([X_train_all, y_train_all,\
                     X_test_all, y_test_all,\
                     xgb_model_4], f)
    ```py
    

7. 测试集分析、财务洞察和交付给客户

活动 7.01：推导财务洞察

解决方案：

1. 使用测试集，计算如果没有辅导程序的所有违约成本。

   使用此代码进行计算：

   cost_of_defaults = np.sum(y_test_all * X_test_all[:,5])
   cost_of_defaults 
   ```py
   
   输出应为：
   
   

   60587763.0

2. 计算辅导项目可以降低违约成本的百分比。

   违约成本的潜在降低是辅导项目可能带来的最大净节省量，除以在没有辅导项目情况下所有违约的成本：

   net_savings[max_savings_ix]/cost_of_defaults
   ```py
   
   输出应为：
   
   

   0.2214260658542551

   
   结果表明，使用辅导程序可以将违约成本降低 22%，这一结果是通过预测建模得出的。
   
   

3. 计算在最优阈值下每个账户的净节省（考虑所有可能进行辅导的账户，也就是相对于整个测试集）。

   使用此代码进行计算：

   net_savings[max_savings_ix]/len(y_test_all)
   ```py
   
   输出应如下所示：
   
   

   2259.2977433479286

   
   这样的结果帮助客户根据辅导项目的潜在节省量，推算其可以为服务的所有账户节省的金额。
   
   

4. 绘制每个阈值下，每个账户的净节省与每个账户辅导成本的关系图。

   使用此代码创建图表：

   plt.plot(total_cost/len(y_test_all),
            net_savings/len(y_test_all))
   plt.xlabel\
   ('Upfront investment: cost of counselings per account (NT$)')
   plt.ylabel('Net savings per account (NT$)')
   ```py
   
   生成的图应如下所示：
   
   ![图 7.14：需要的辅导项目初始成本    达到给定节省量    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/25bc7b6b389840fa9c8940150dbaadb1~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=DyyRBUFiLDllUhHJHqhwfOL0fhQ%3D)
   
   图 7.14：实现一定节省金额所需的咨询项目初始成本
   
   这表明客户需要在某个月份为咨询项目预算多少资金，以实现一定的节省金额。看起来最大利益可以通过为每个账户预算大约 NT$1300 来实现（你可以使用 `np.argmax` 找到对应于最大净节省的精确预算金额）。然而，在 NT$1000 到 2000 之间，前期投资的净节省相对平稳，超出此范围则较低。客户实际上可能无法为该项目预算这么多资金，但这个图表为他们提供了证据，如果需要的话可以为更大预算进行辩护。
   
   这个结果与我们在上一个练习中的图形相对应。虽然我们显示了最优阈值是 0.36，但对客户而言，使用更高的阈值（高至 0.5）可能是可以接受的，这样就会做出更少的正预测，向更少的账户持有者提供咨询，并且前期项目成本较小。*图 7.14* 显示了这一点在成本和每个账户的净节省方面的体现。
   
   

5. 在每个阈值下，绘制预测为正的账户比例（这称为“标记率”）。

   使用此代码绘制标记率与阈值的关系：

   plt.plot(thresholds, n_pos_pred/len(y_test_all))
   plt.ylabel('Flag rate')
   plt.xlabel('Threshold')
   ```py
   
   图表应如下所示：
   
   ![图 7.15：信用咨询项目的标记率与阈值关系    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/58becb4d320c48f08c626c3738345ecd~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=Px%2Bp60fGRxpxjD4rVyW%2F%2F06BfgM%3D)
   
   图 7.15：信用咨询项目的标记率与阈值关系
   
   这个图表显示了在每个阈值下，将被预测为违约并因此会被推荐接触的账户比例。看起来在最优阈值 0.36 时，只有大约 20% 的账户会被标记为需要咨询。这表明，使用模型优先考虑需要咨询的账户可以帮助集中资源，减少资源浪费。更高的阈值，可能会导致几乎最优的节省，直到约 0.5 的阈值，如 *图 7.12* 所示（*第七章*，*测试集分析、财务洞察及交付给客户*）会导致更低的标记率。
   
   

6. 使用以下代码为测试数据绘制精度-召回率曲线：

   plt.plot(n_true_pos/sum(y_test_all),\
            np.divide(n_true_pos, n_pos_pred)) 
   plt.xlabel('Recall')
   plt.ylabel('Precision')
   ```py
   
   该图应如下所示：
   
   ![图 7.16：精度-召回率曲线    ](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/3f942e7a089643ff98953d41eb532ea3~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1771551349&x-signature=k%2BmGVuSRNDQGvOcKzE2w14C%2BuJg%3D)
   
   图 7.16：精度-召回率曲线
   
   *图 7.16* 显示，为了开始获得高于 0 的真实正率（即召回率），我们需要接受约 0.8 或更低的精准度。
   
   精度和召回率与项目的成本和节省直接相关：我们预测越精准，由于模型预测错误而浪费在咨询上的钱就越少。而且，召回率越高，我们就能通过成功识别即将违约的账户创造更多节省。将此步骤中的代码与上一个练习中用于计算成本和节省的代码进行比较，即可看到这一点。
   
   为了查看精度和召回率与定义正负预测的阈值之间的关系，分别绘制它们可能会很有启发性。
   
   

7. 将精度和召回率分别绘制在y轴上，阈值绘制在x轴上。

   使用这个代码来生成图表：

   plt.plot(thresholds, np.divide(n_true_pos, n_pos_pred),
            label='Precision')
   plt.plot(thresholds, n_true_pos/sum(y_test_all),
            label='Recall')
   plt.xlabel('Threshold')
   plt.legend()
   

   图表应如下所示：

   

图 7.17：精度和召回率分别绘制与阈值的关系

这个图表揭示了为何最佳阈值为 0.36 的原因。虽然最佳阈值还取决于成本和节省的财务分析，但我们可以看到，精度最初增加的陡峭部分，代表了正向预测的准确性，因此也反映了模型引导的咨询有多具成本效益，发生在大约 0.36 的阈值之前。

嘿！

我是 Stephen Klosterman，本书的作者。我非常希望你喜欢阅读我的书，并且觉得它对你有所帮助。

如果你能在亚马逊上留下关于《使用 Python 进行数据科学项目》第二版的评论，分享你的想法，将对我（和其他潜在读者！）非常有帮助。

请访问链接packt.link/r/1800564481。

或者

扫描二维码留下你的评论。

你的评论将帮助我了解这本书中哪些地方做得好，哪些地方可以改进，以便为未来的版本做出改进，因此我非常感激。

最好的祝愿，

Stephen Klosterman