注意
跳转至页面底部 下载完整的示例代码，或通过 JupyterLite 或 Binder 在浏览器中运行此示例。
偏相关与个体条件期望图#
偏相关图（Partial dependence plots）展示了目标函数 [2] 与一组关注特征之间的依赖关系，同时对所有其他特征（补充特征）的值进行边缘化处理。由于人类感知的局限性，关注特征集的规模必须很小（通常为一到两个），因此通常在最重要的特征中进行选择。
类似地，个体条件期望（ICE）图 [3] 展示了目标函数与某个关注特征之间的依赖关系。然而，与展示关注特征平均效应的偏相关图不同，ICE 图通过为每个 样本 单独绘制一条线，来可视化预测结果对该特征的依赖关系。ICE 图仅支持单个关注特征。
本示例展示了如何从在共享单车数据集上训练的 MLPRegressor 和 HistGradientBoostingRegressor 中获取偏相关图和 ICE 图。本示例受到 [1] 的启发。
# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause
共享单车数据集预处理#
我们将使用共享单车数据集。其目标是利用天气、季节数据以及日期时间信息来预测单车租赁数量。
from sklearn.datasets import fetch_openml
bikes = fetch_openml("Bike_Sharing_Demand", version=2, as_frame=True)
# Make an explicit copy to avoid "SettingWithCopyWarning" from pandas
X, y = bikes.data.copy(), bikes.target
# We use only a subset of the data to speed up the example.
X = X.iloc[::5, :]
y = y[::5]
"weather"（天气）特征有一个特性：类别 "heavy_rain"（大雨）是一个稀有类别。
X["weather"].value_counts()
weather
clear 2284
misty 904
rain 287
heavy_rain 1
Name: count, dtype: int64
由于这个稀有类别的存在，我们将其合并入 "rain"（雨天）。
X["weather"] = (
X["weather"]
.astype(object)
.replace(to_replace="heavy_rain", value="rain")
.astype("category")
)
现在我们更仔细地查看 "year"（年份）特征。
X["year"].value_counts()
year
1 1747
0 1729
Name: count, dtype: int64
我们可以看到数据涵盖了两年。我们使用第一年的数据来训练模型，使用第二年的数据来测试模型。
mask_training = X["year"] == 0.0
X = X.drop(columns=["year"])
X_train, y_train = X[mask_training], y[mask_training]
X_test, y_test = X[~mask_training], y[~mask_training]
通过查看数据集信息，我们可以确认存在异构数据类型。我们需要相应地预处理不同的列。
X_train.info()
<class 'pandas.DataFrame'>
RangeIndex: 1729 entries, 0 to 8640
Data columns (total 11 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 season 1729 non-null category
1 month 1729 non-null int64
2 hour 1729 non-null int64
3 holiday 1729 non-null category
4 weekday 1729 non-null int64
5 workingday 1729 non-null category
6 weather 1729 non-null category
7 temp 1729 non-null float64
8 feel_temp 1729 non-null float64
9 humidity 1729 non-null float64
10 windspeed 1729 non-null float64
dtypes: category(4), float64(4), int64(3)
memory usage: 101.5 KB
根据前面的信息，我们将 category（类别）列视为名义分类特征。此外，我们也将日期和时间信息视为分类特征。
我们手动定义包含数值特征和分类特征的列。
numerical_features = [
"temp",
"feel_temp",
"humidity",
"windspeed",
]
categorical_features = X_train.columns.drop(numerical_features)
在深入研究不同机器学习流水线的预处理细节之前，我们将尝试对数据集获得一些额外的直觉，这将有助于理解模型的统计表现以及偏相关分析的结果。
我们将通过按季节和年份对数据进行分组，绘制单车租赁的平均数量。
from itertools import product
import matplotlib.pyplot as plt
import numpy as np
days = ("Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat")
hours = tuple(range(24))
xticklabels = [f"{day}\n{hour}:00" for day, hour in product(days, hours)]
xtick_start, xtick_period = 6, 12
fig, axs = plt.subplots(nrows=2, figsize=(8, 6), sharey=True, sharex=True)
average_bike_rentals = bikes.frame.groupby(
["year", "season", "weekday", "hour"], observed=True
).mean(numeric_only=True)["count"]
for ax, (idx, df) in zip(axs, average_bike_rentals.groupby("year")):
df.groupby("season", observed=True).plot(ax=ax, legend=True)
# decorate the plot
ax.set_xticks(
np.linspace(
start=xtick_start,
stop=len(xticklabels),
num=len(xticklabels) // xtick_period,
)
)
ax.set_xticklabels(xticklabels[xtick_start::xtick_period])
ax.set_xlabel("")
ax.set_ylabel("Average number of bike rentals")
ax.set_title(
f"Bike rental for {'2010 (train set)' if idx == 0.0 else '2011 (test set)'}"
)
ax.set_ylim(0, 1_000)
ax.set_xlim(0, len(xticklabels))
ax.legend(loc=2)
训练集和测试集之间的第一个显著差异是测试集中的单车租赁数量更高。因此，得到一个低估单车租赁数量的机器学习模型并不令人意外。我们还观察到，春季的单车租赁数量较低。此外，我们看到在工作日期间，上午 6-7 点和下午 5-6 点左右存在特定的租赁高峰模式。我们可以记住这些不同的洞察，并利用它们来理解偏相关图。
机器学习模型的预处理器#
由于我们稍后将使用两种不同的模型，即 我们将使用 对于梯度提升模型，我们保留数值特征不变，仅使用 MLPRegressor 和 HistGradientBoostingRegressor，我们分别为每个模型创建了两个不同的预处理器。
神经网络模型的预处理器#
QuantileTransformer 来缩放数值特征，并使用 OneHotEncoder 对分类特征进行编码。from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder, QuantileTransformer
mlp_preprocessor = ColumnTransformer(
transformers=[
("num", QuantileTransformer(n_quantiles=100), numerical_features),
("cat", OneHotEncoder(handle_unknown="ignore"), categorical_features),
]
)
mlp_preprocessor
梯度提升模型的预处理器#
OrdinalEncoder 对分类特征进行编码。from sklearn.preprocessing import OrdinalEncoder
hgbdt_preprocessor = ColumnTransformer(
transformers=[
("cat", OrdinalEncoder(), categorical_features),
("num", "passthrough", numerical_features),
],
sparse_threshold=1,
verbose_feature_names_out=False,
).set_output(transform="pandas")
hgbdt_preprocessor
不同模型的 1 向偏相关分析#
在本节中，我们将使用两种不同的机器学习模型计算 1 向偏相关：(i) 多层感知机，(ii) 梯度提升模型。通过这两个模型，我们说明如何计算和解释数值特征和分类特征的偏相关图（PDP）以及个体条件期望（ICE）。
多层感知机#
让我们拟合一个 MLPRegressor 并计算单变量偏相关图。
from time import time
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
print("Training MLPRegressor...")
tic = time()
mlp_model = make_pipeline(
mlp_preprocessor,
MLPRegressor(
hidden_layer_sizes=(30, 15),
learning_rate_init=0.01,
early_stopping=True,
random_state=0,
),
)
mlp_model.fit(X_train, y_train)
print(f"done in {time() - tic:.3f}s")
print(f"Test R2 score: {mlp_model.score(X_test, y_test):.2f}")
Training MLPRegressor...
done in 0.745s
Test R2 score: 0.61
我们使用专门为神经网络创建的预处理器配置了一个流水线，并调整了神经网络的大小和学习率，以在训练时间和测试集预测性能之间取得合理的折中。
重要的是，此表格数据集的特征具有非常不同的动态范围。神经网络往往对具有不同尺度的特征非常敏感，如果忘记预处理数值特征，将导致模型性能非常差。
如果使用更大的神经网络，可能会获得更高的预测性能，但训练成本也会显著增加。
请注意，在绘制偏相关图之前，检查模型在测试集上的准确性是否足够高非常重要，因为解释预测性能较差的模型中给定特征对预测函数的影响意义不大。在这方面，我们的 MLP 模型工作得相当好。
我们将绘制平均偏相关。
import matplotlib.pyplot as plt
from sklearn.inspection import PartialDependenceDisplay
common_params = {
"subsample": 50,
"n_jobs": 2,
"grid_resolution": 20,
"random_state": 0,
}
print("Computing partial dependence plots...")
features_info = {
# features of interest
"features": ["temp", "humidity", "windspeed", "season", "weather", "hour"],
# type of partial dependence plot
"kind": "average",
# information regarding categorical features
"categorical_features": categorical_features,
}
tic = time()
_, ax = plt.subplots(ncols=3, nrows=2, figsize=(9, 8), constrained_layout=True)
display = PartialDependenceDisplay.from_estimator(
mlp_model,
X_train,
**features_info,
ax=ax,
**common_params,
)
print(f"done in {time() - tic:.3f}s")
_ = display.figure_.suptitle(
(
"Partial dependence of the number of bike rentals\n"
"for the bike rental dataset with an MLPRegressor"
),
fontsize=16,
)
Computing partial dependence plots...
done in 0.956s
梯度提升#
现在让我们拟合一个 HistGradientBoostingRegressor 并计算相同特征上的偏相关。我们同样使用为该模型创建的特定预处理器。
from sklearn.ensemble import HistGradientBoostingRegressor
print("Training HistGradientBoostingRegressor...")
tic = time()
hgbdt_model = make_pipeline(
hgbdt_preprocessor,
HistGradientBoostingRegressor(
categorical_features=categorical_features,
random_state=0,
max_iter=50,
),
)
hgbdt_model.fit(X_train, y_train)
print(f"done in {time() - tic:.3f}s")
print(f"Test R2 score: {hgbdt_model.score(X_test, y_test):.2f}")
Training HistGradientBoostingRegressor...
done in 0.118s
Test R2 score: 0.62
在这里，我们对梯度提升模型使用了默认超参数，且没有进行任何预处理，因为树模型天然对数值特征的单调变换具有鲁棒性。
请注意，在此表格数据集上，梯度提升机（Gradient Boosting Machines）的训练速度比神经网络快得多，且精度更高。调整其超参数也便宜得多（默认设置通常工作良好，而神经网络则不然）。
我们将绘制部分数值特征和分类特征的偏相关图。
print("Computing partial dependence plots...")
tic = time()
_, ax = plt.subplots(ncols=3, nrows=2, figsize=(9, 8), constrained_layout=True)
display = PartialDependenceDisplay.from_estimator(
hgbdt_model,
X_train,
**features_info,
ax=ax,
**common_params,
)
print(f"done in {time() - tic:.3f}s")
_ = display.figure_.suptitle(
(
"Partial dependence of the number of bike rentals\n"
"for the bike rental dataset with a gradient boosting"
),
fontsize=16,
)
Computing partial dependence plots...
done in 1.421s
图表分析#
我们将首先观察数值特征的 PDP。对于两个模型，温度 PDP 的总体趋势是单车租赁数量随温度升高而增加。我们可以对湿度特征进行类似的分析，但趋势相反：随着湿度增加，单车租赁数量会减少。最后，我们看到风速特征也呈现相同的趋势：对于两个模型，随着风速增加，单车租赁数量都会减少。我们还观察到 MLPRegressor 比 HistGradientBoostingRegressor 的预测平滑得多。
现在，我们将查看分类特征的偏相关图。
我们观察到春季是季节特征中的最低柱。在天气特征中，雨天类别是最低柱。关于小时特征，我们看到上午 7 点和下午 6 点左右有两个高峰。这些发现与我们之前在数据集上所做的观察一致。
然而，值得注意的是，如果特征之间存在相关性，我们可能会创建意义不明的合成样本。
ICE 与 PDP 对比#
PDP 是特征边缘效应的平均值。我们正在对提供集合的所有样本的响应进行平均。因此，一些效应可能会被隐藏。在这方面，可以绘制每个单独的响应。这种表示被称为个体效应图（ICE）。在下图中，我们绘制了针对温度和湿度特征随机选择的 50 个 ICE。
print("Computing partial dependence plots and individual conditional expectation...")
tic = time()
_, ax = plt.subplots(ncols=2, figsize=(6, 4), sharey=True, constrained_layout=True)
features_info = {
"features": ["temp", "humidity"],
"kind": "both",
"centered": True,
}
display = PartialDependenceDisplay.from_estimator(
hgbdt_model,
X_train,
**features_info,
ax=ax,
**common_params,
)
print(f"done in {time() - tic:.3f}s")
_ = display.figure_.suptitle("ICE and PDP representations", fontsize=16)
Computing partial dependence plots and individual conditional expectation...
done in 0.607s
我们看到温度特征的 ICE 为我们提供了一些额外信息：一些 ICE 线是平坦的，而另一些则显示在摄氏 35 度以上对温度的依赖性下降。我们观察到湿度特征有类似的模式：一些 ICE 线显示当湿度高于 80% 时有明显的下降。
并非所有的 ICE 线都是平行的，这表明模型发现了特征之间的相互作用。我们可以通过使用 interaction_cst 参数限制梯度提升模型不使用特征间的任何交互来重复该实验。
from sklearn.base import clone
interaction_cst = [[i] for i in range(X_train.shape[1])]
hgbdt_model_without_interactions = (
clone(hgbdt_model)
.set_params(histgradientboostingregressor__interaction_cst=interaction_cst)
.fit(X_train, y_train)
)
print(f"Test R2 score: {hgbdt_model_without_interactions.score(X_test, y_test):.2f}")
Test R2 score: 0.38
_, ax = plt.subplots(ncols=2, figsize=(6, 4), sharey=True, constrained_layout=True)
features_info["centered"] = False
display = PartialDependenceDisplay.from_estimator(
hgbdt_model_without_interactions,
X_train,
**features_info,
ax=ax,
**common_params,
)
_ = display.figure_.suptitle("ICE and PDP representations", fontsize=16)
2D 交互作用图#
具有两个关注特征的 PDP 使我们能够可视化它们之间的交互作用。然而，ICE 不能以简单的方式绘制和解释。我们将展示 from_estimator 中可用的 2D 热图表示。
print("Computing partial dependence plots...")
features_info = {
"features": ["temp", "humidity", ("temp", "humidity")],
"kind": "average",
}
_, ax = plt.subplots(ncols=3, figsize=(10, 4), constrained_layout=True)
tic = time()
display = PartialDependenceDisplay.from_estimator(
hgbdt_model,
X_train,
**features_info,
ax=ax,
**common_params,
)
print(f"done in {time() - tic:.3f}s")
_ = display.figure_.suptitle(
"1-way vs 2-way of numerical PDP using gradient boosting", fontsize=16
)
Computing partial dependence plots...
done in 9.457s
双向偏相关图显示了单车租赁数量对温度和湿度联合值的依赖关系。我们清楚地看到了这两个特征之间的相互作用。对于高于摄氏 20 度的温度，湿度对单车租赁数量的影响看起来与温度无关。
另一方面，对于低于摄氏 20 度的温度，温度和湿度都会持续影响单车租赁数量。
此外，摄氏 20 度阈值的影响脊线的斜率与湿度水平高度相关：在干燥条件下脊线陡峭，但在高于 70% 湿度的潮湿条件下则平缓得多。
我们现在将这些结果与为约束模型计算的相同图表进行对比，该模型被限制为学习不依赖此类非线性特征交互的预测函数。
print("Computing partial dependence plots...")
features_info = {
"features": ["temp", "humidity", ("temp", "humidity")],
"kind": "average",
}
_, ax = plt.subplots(ncols=3, figsize=(10, 4), constrained_layout=True)
tic = time()
display = PartialDependenceDisplay.from_estimator(
hgbdt_model_without_interactions,
X_train,
**features_info,
ax=ax,
**common_params,
)
print(f"done in {time() - tic:.3f}s")
_ = display.figure_.suptitle(
"1-way vs 2-way of numerical PDP using gradient boosting", fontsize=16
)
Computing partial dependence plots...
done in 9.690s
对于约束为不建模特征交互的模型，1D 偏相关图显示每个特征单独存在局部尖峰，特别是在“湿度”特征方面。这些尖峰可能反映了模型的降级行为，即模型试图通过过度拟合特定的训练点来以某种方式补偿被禁止的交互。请注意，该模型在测试集上测得的预测性能明显低于原始的无约束模型。
还要注意，这些图表上可见的局部尖峰数量取决于 PD 图本身的网格分辨率参数。
这些局部尖峰导致了 2D PD 图的网格化噪声。由于湿度特征中的高频振荡，很难断定这些特征之间是否没有交互作用。然而，可以清楚地看到，当温度跨越 20 度边界时观察到的简单交互效应在该模型中不再可见。
分类特征之间的偏相关将提供可以显示为热图的离散表示。例如，季节、天气和目标之间的交互如下：
print("Computing partial dependence plots...")
features_info = {
"features": ["season", "weather", ("season", "weather")],
"kind": "average",
"categorical_features": categorical_features,
}
_, ax = plt.subplots(ncols=3, figsize=(14, 6), constrained_layout=True)
tic = time()
display = PartialDependenceDisplay.from_estimator(
hgbdt_model,
X_train,
**features_info,
ax=ax,
**common_params,
)
print(f"done in {time() - tic:.3f}s")
_ = display.figure_.suptitle(
"1-way vs 2-way PDP of categorical features using gradient boosting", fontsize=16
)
Computing partial dependence plots...
done in 0.534s
3D 表示#
让我们为 2 特征交互制作相同的偏相关图，这次使用 3 维展示。
# unused but required import for doing 3d projections with matplotlib < 3.2
import mpl_toolkits.mplot3d # noqa: F401
import numpy as np
from sklearn.inspection import partial_dependence
fig = plt.figure(figsize=(5.5, 5))
features = ("temp", "humidity")
pdp = partial_dependence(
hgbdt_model, X_train, features=features, kind="average", grid_resolution=10
)
XX, YY = np.meshgrid(pdp["grid_values"][0], pdp["grid_values"][1])
Z = pdp.average[0].T
ax = fig.add_subplot(projection="3d")
fig.add_axes(ax)
surf = ax.plot_surface(XX, YY, Z, rstride=1, cstride=1, cmap=plt.cm.BuPu, edgecolor="k")
ax.set_xlabel(features[0])
ax.set_ylabel(features[1])
fig.suptitle(
"PD of number of bike rentals on\nthe temperature and humidity GBDT model",
fontsize=16,
)
# pretty init view
ax.view_init(elev=22, azim=122)
clb = plt.colorbar(surf, pad=0.08, shrink=0.6, aspect=10)
clb.ax.set_title("Partial\ndependence")
plt.show()
自定义检查点#
到目前为止的示例都没有指定为了创建偏相关图而评估了_哪些_点。默认情况下，我们使用输入数据集定义的百分位数。在某些情况下，指定您希望评估模型的精确点会很有帮助。例如，如果用户想要测试模型在分布外数据上的表现，或比较在稍微不同的数据上拟合的两个模型。custom_values 参数允许用户传入他们希望模型评估的值。这会覆盖 grid_resolution 和 percentiles 参数。让我们回到上面的梯度提升示例，但使用自定义值。
print("Computing partial dependence plots with custom evaluation values...")
tic = time()
_, ax = plt.subplots(ncols=2, figsize=(6, 4), sharey=True, constrained_layout=True)
features_info = {
"features": ["temp", "humidity"],
"kind": "both",
}
display = PartialDependenceDisplay.from_estimator(
hgbdt_model,
X_train,
**features_info,
ax=ax,
**common_params,
# we set custom values for temp feature -
# all other features are evaluated based on the data
custom_values={"temp": np.linspace(0, 40, 10)},
)
print(f"done in {time() - tic:.3f}s")
_ = display.figure_.suptitle(
(
"Partial dependence of the number of bike rentals\n"
"for the bike rental dataset with a gradient boosting"
),
fontsize=16,
)
Computing partial dependence plots with custom evaluation values...
done in 0.728s
脚本总运行时间： (0 分 30.470 秒)
相关示例
