注意

转到末尾以下载完整示例代码。或者通过JupyterLite或Binder在浏览器中运行此示例

局部依赖图和个体条件期望图#

局部依赖图显示目标函数[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.core.frame.DataFrame'>
Index: 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: 115.4 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)
Bike rental for 2010 (train set), Bike rental for 2011 (test set)

训练集和测试集之间第一个显著差异是测试集中的单车租赁数量更高。因此，机器学习模型低估单车租赁数量也就不足为奇了。我们还观察到春季的单车租赁数量较低。此外，我们发现工作日期间，上午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

ColumnTransformer(transformers=[('num', QuantileTransformer(n_quantiles=100),
                                 ['temp', 'feel_temp', 'humidity',
                                  'windspeed']),
                                ('cat', OneHotEncoder(handle_unknown='ignore'),
                                 Index(['season', 'month', 'hour', 'holiday', 'weekday', 'workingday',
       'weather'],
      dtype='object'))])
在Jupyter环境中，请重新运行此单元格以显示HTML表示或信任此笔记本。
在GitHub上，HTML表示无法渲染，请尝试使用nbviewer.org加载此页面。


梯度提升模型预处理器#

对于梯度提升模型，我们保持数值特征不变，只使用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

ColumnTransformer(sparse_threshold=1,
                  transformers=[('cat', OrdinalEncoder(),
                                 Index(['season', 'month', 'hour', 'holiday', 'weekday', 'workingday',
       'weather'],
      dtype='object')),
                                ('num', 'passthrough',
                                 ['temp', 'feel_temp', 'humidity',
                                  'windspeed'])],
                  verbose_feature_names_out=False)
在Jupyter环境中，请重新运行此单元格以显示HTML表示或信任此笔记本。
在GitHub上，HTML表示无法渲染，请尝试使用nbviewer.org加载此页面。


使用不同模型的单向局部依赖#

在本节中，我们将使用两种不同的机器学习模型计算单向局部依赖：(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.591s
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,
)
Partial dependence of the number of bike rentals for the bike rental dataset with an MLPRegressor
Computing partial dependence plots...
done in 0.549s

梯度提升#

现在让我们拟合一个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.123s
Test R2 score: 0.62

在这里，我们使用了梯度提升模型的默认超参数，没有进行任何预处理，因为基于树的模型对数值特征的单调变换具有天然的鲁棒性。

请注意，在此表格数据集上，梯度提升机在训练速度和准确性方面都显著优于神经网络。调整其超参数的成本也显著降低（默认值通常效果良好，而神经网络则不然）。

我们将绘制一些数值和类别特征的局部依赖图。

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,
)
Partial dependence of the number of bike rentals for the bike rental dataset with a gradient boosting
Computing partial dependence plots...
done in 1.015s

图分析#

我们将首先查看数值特征的PDPs。对于这两个模型，温度PDP的总体趋势是单车租赁数量随温度升高而增加。我们可以对湿度特征进行类似但趋势相反的分析。湿度增加时，单车租赁数量减少。最后，我们看到风速特征也呈现相同趋势。当风速增加时，两个模型的单车租赁数量都减少。我们还观察到MLPRegressor的预测比HistGradientBoostingRegressor平滑得多。

现在，我们将查看类别特征的局部依赖图。

我们观察到，对于季节特征，春季是最低的。对于天气特征，降雨类别是最低的。关于小时特征，我们看到上午7点和下午6点左右有两个高峰。这些发现与我们之前在数据集上进行的观察结果一致。

然而，值得注意的是，如果特征相关，我们可能会创建潜在的无意义合成样本。

ICE vs. 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)
ICE and PDP representations
Computing partial dependence plots and individual conditional expectation...
done in 0.471s

我们看到温度特征的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)
ICE and PDP representations

2D交互图#

具有两个感兴趣特征的PDPs使我们能够可视化它们之间的交互作用。然而，ICEs无法以简单的方式绘制和解释。我们将展示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
)
1-way vs 2-way of numerical PDP using gradient boosting
Computing partial dependence plots...
done in 7.227s

双向局部依赖图显示了单车租赁数量对温度和湿度的联合值的依赖性。我们清楚地看到了这两个特征之间的交互作用。当温度高于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
)
1-way vs 2-way of numerical PDP using gradient boosting
Computing partial dependence plots...
done in 6.518s

限制不建模特征交互的模型的一维局部依赖图显示了每个特征的局部尖峰，特别是对于“湿度”特征。这些尖峰可能反映了模型的一种退化行为，它试图通过过拟合特定训练点来某种程度上补偿被禁止的交互。请注意，该模型在测试集上测量的预测性能显著低于原始的、不受约束的模型。

另请注意，这些图上可见的局部尖峰数量取决于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
)
1-way vs 2-way PDP of categorical features using gradient boosting
Computing partial dependence plots...
done in 0.317s

3D表示#

让我们为这两个特征交互绘制相同的局部依赖图，这次以3D形式呈现。

# 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()
PD of number of bike rentals on the temperature and humidity GBDT model, Partial dependence

自定义检查点#

到目前为止，所有示例都没有指定评估哪些点来创建局部依赖图。默认情况下，我们使用输入数据集定义的分位数。在某些情况下，指定您希望模型评估的精确点可能会有所帮助。例如，如果用户想要测试模型在分布外数据上的行为，或者比较在略微不同数据上拟合的两个模型。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,
)
Partial dependence of the number of bike rentals for the bike rental dataset with a gradient boosting
Computing partial dependence plots with custom evaluation values...
done in 0.425s

脚本总运行时间： (0分钟 22.278秒)

相关示例

局部依赖的高级绘图

局部依赖的高级绘图

时间相关特征工程

时间相关特征工程

scikit-learn 0.24 的发布亮点

scikit-learn 0.24 的发布亮点

scikit-learn 1.2 的发布亮点

scikit-learn 1.2 的发布亮点

由Sphinx-Gallery生成