用 PDP 与 ICE 读懂模型:平均趋势、个体差异与特征交互

模型认为气温变高时租车量会怎样变化?仅看一个特征的重要性排名,回答不了这个问题。部分依赖图(Partial Dependence Plot,PDP)展示所关注特征与拟合预测函数之间的关系,同时对其他特征取平均;个体条件期望图(Individual Conditional Expectation,ICE)则为每个样本画出一条曲线,展示平均值可能掩盖的差异。

本文译编自 scikit-learn 官方完整实例,作者为 The scikit-learn developers,许可为 BSD-3-Clause。2026 年 10 月 5 日核验的页面标为 scikit-learn 1.9.1;已连读该页与可下载 Python 源文件。下文图像和数值均来自官方实例,编辑过程没有运行训练或绘图。API 以该版本页面为依据,旧版本未必接受所有参数。

PDP 平均了什么,ICE 保留了什么

以单个特征 x 为例,选定一个网格值,把参考数据集中每条样本的这个特征替换为该值,其余特征保持原样,再计算模型预测。将所有预测平均,就得到这个网格点的 PDP;沿多个网格点重复,即得到一条曲线。ICE 保留每条样本自己的预测轨迹,不在样本维度上平均。PDP 常展示一到两个特征,ICE 在此接口中只支持单个关注特征。

这些图解释的是拟合函数,不是对现实世界进行干预的因果效应。如果特征相关,替换单个特征可能制造现实中几乎不存在的组合,例如温度与体感温度明显不匹配。图上得到的是模型在这些合成输入上的行为,不自动等于真实租车需求会怎样改变。分类任务的目标响应还涉及决策分数或概率等选择;本文实例使用回归。

准备共享单车数据,并按年份留出测试集

实例使用 OpenML 的 Bike_Sharing_Demand 第 2 版,根据天气、季节和日期时间预测租车数量。为缩短示例耗时,每隔五条取一条。显式复制特征表是为了避免 pandas 的切片赋值歧义。源实例中 heavy_rain 在抽样后只出现一次,因此合并为 rain。

from sklearn.datasets import fetch_openml

bikes = fetch_openml("Bike_Sharing_Demand", version=2, as_frame=True)
X, y = bikes.data.copy(), bikes.target
X = X.iloc[::5, :]
y = y[::5]
X["weather"] = (
    X["weather"].astype(object)
    .replace(to_replace="heavy_rain", value="rain")
    .astype("category")
)

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]

numerical_features = ["temp", "feel_temp", "humidity", "windspeed"]
categorical_features = X_train.columns.drop(numerical_features)

用第一个年份代码对应的数据训练,第二个年份代码对应的数据测试,比把相邻小时随机混在训练集与测试集中更接近跨时间预测。这里不根据脚本图题推断实际公历年份;源脚本画图的年份文字是硬编码的,开展自己的复现时应与数据说明核实。月份、小时、星期、假期和工作日等变量在这个示例里被视为类别;温度、体感温度、湿度和风速为数值特征。

在建模前,原文还按 year、season、weekday、hour 汇总平均租车量:

average_bike_rentals = bikes.frame.groupby(
    ["year", "season", "weekday", "hour"], observed=True
).mean(numeric_only=True)["count"]

这一描述性图使用完整原始表,而模型用的是每五条取一条的子集。原文观察到第二年的租车量更高,因此用第一年训练的模型可能低估测试期;春季租车少,工作日早上约 6–7 点和傍晚约 5–6 点出现峰值。这些先验观察有助于理解模型图,但使用测试期探索后,不应把全部后续选择称作完全盲测。

给两类模型准备不同的输入变换

神经网络对不同特征的数值尺度敏感。官方实例给数值列使用 QuantileTransformer(n_quantiles=100),给类别列使用 one-hot 编码,并让未知类别在转换时被忽略。预处理器放进 Pipeline,确保变换只在训练数据上拟合。

from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder, QuantileTransformer, OrdinalEncoder

mlp_preprocessor = ColumnTransformer([
    ("num", QuantileTransformer(n_quantiles=100), numerical_features),
    ("cat", OneHotEncoder(handle_unknown="ignore"), categorical_features),
])

hgbdt_preprocessor = ColumnTransformer(
    [
        ("cat", OrdinalEncoder(), categorical_features),
        ("num", "passthrough", numerical_features),
    ],
    sparse_threshold=1,
    verbose_feature_names_out=False,
).set_output(transform="pandas")

直方图梯度提升模型保持数值列不变,只给类别列做序数编码。输出设为 pandas,并保留原列名,让模型可以按名称识别类别特征;这些整数是类别标签,不应被误解为连续大小关系。原文说树模型“不做预处理”,应理解为数值列不做尺度变换:代码实际上仍然编码了类别列。这里的 OrdinalEncoder() 默认遇到未见类别会报错,真实应用需事先设计未知类别处理方式并验证,不能从示例数据没有报错推导出生产数据也不会报错。

先检查预测能力,再解释函数

from sklearn.neural_network import MLPRegressor
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.pipeline import make_pipeline

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)

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(mlp_model.score(X_test, y_test))
print(hgbdt_model.score(X_test, y_test))

所核验官方页面报告 MLP 的测试 R² 为 0.61,梯度提升为 0.62;这些是原文输出,不是本次测试结论。原文也报告梯度提升在这份表格数据上训练更快,但速度结论依赖硬件、依赖版本和具体模型设置,不能外推为所有任务的规则。

解释一个预测性能很差的模型,仍能发现它学到了什么,但不一定能说明问题本身的合理规律。原文明确承认神经网络大小和学习率曾按测试表现折中选择。这适合作为模型解释演示;若要正式报告泛化性能,应另设验证集调参,再用新的未参与选择的测试集评估。early_stopping=True 的内部验证也不自动变成时间顺序验证。

画数值与类别的一维 PDP

import matplotlib.pyplot as plt
from sklearn.inspection import PartialDependenceDisplay

common_params = {
    "subsample": 50,
    "n_jobs": 2,
    "grid_resolution": 20,
    "random_state": 0,
}
features_info = {
    "features": ["temp", "humidity", "windspeed", "season", "weather", "hour"],
    "kind": "average",
    "categorical_features": categorical_features,
}

for model, name in [(mlp_model, "MLP"), (hgbdt_model, "Gradient boosting")]:
    _, ax = plt.subplots(nrows=2, ncols=3, figsize=(9, 8), constrained_layout=True)
    display = PartialDependenceDisplay.from_estimator(
        model, X_train, **features_info, ax=ax, **common_params
    )
    display.figure_.suptitle(name)

这段把原文分别绘图的重复代码合并成循环,模型和关注特征不变。传入的是原始特征表与完整 Pipeline,而非已经编码过的矩阵,因此横轴仍可以用可理解的原特征名称。grid_resolution=20 控制连续特征网格;类别特征按类别显示。subsample=50 主要控制绘出的 ICE 曲线数量,不表示平均 PDP 只使用 50 个样本。

原文图中,两个模型总体都表现为气温升高时预测租车量增加,湿度或风速升高时减少。MLP 的曲线比梯度提升平滑;类别 PDP 中春季最低、雨天天气类别最低,小时特征在约 7 点与 18 点出现峰值。这些是经其他特征平均后的模型响应,不是对每个个体都成立的关系。

用 ICE 看见平均曲线隐藏的差异

_, ax = plt.subplots(ncols=2, figsize=(6, 4), sharey=True, constrained_layout=True)
display = PartialDependenceDisplay.from_estimator(
    hgbdt_model,
    X_train,
    features=["temp", "humidity"],
    kind="both",
    centered=True,
    ax=ax,
    **common_params,
)
官方实例中的居中ICE和平均PDP:温度与湿度各有多条个体曲线,部分曲线平坦,部分在高温或高湿区域明显下降。
原图:The scikit-learn developers,BSD-3-Clause;来自官方实例的 ICE and PDP representations,未重新绘制或伪造运行。

kind="both" 同时显示平均曲线和抽取的 50 条 ICE;centered=True 把曲线按起始网格点居中,让不同样本的变化形状更容易比较。官方图里,有的温度 ICE 近乎平坦,有的在较高温端下降;湿度超过约 80% 时,部分曲线出现明显下落。原文文字把温度下降描述为超过约 35°C,但所核验当前图的温度网格右端约为 33°C,因此不能仅凭这张图确认 35°C 这个阈值。PDP 把不同曲线平均后,也容易隐藏少数群体的特殊响应。

若各 ICE 不平行,说明所学预测函数里存在相应特征与其他特征的交互:改变气温的效果取决于样本的其他属性。这仍是模型交互,不是自动发现了现实中的因果机制。

禁止交互,做一组结构对照

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(hgbdt_model_without_interactions.score(X_test, y_test))

_, ax = plt.subplots(ncols=2, figsize=(6, 4), sharey=True, constrained_layout=True)
PartialDependenceDisplay.from_estimator(
    hgbdt_model_without_interactions,
    X_train,
    features=["temp", "humidity"],
    kind="both",
    centered=False,
    ax=ax,
    **common_params,
)

将每个特征单独放进一个交互约束集合,可以限制树模型学习跨特征交互。本例预处理没有改变列数,因此用转换前列数生成这些单元素集合恰好适用;若改成 one-hot 展开,就必须根据转换后的特征重新建立约束,不能复用这个索引假设。

官方页面报告受限模型的测试 R² 降到 0.38。它的一维 PDP 出现局部尖峰,尤其在湿度方向上。原文提出一种解释:模型可能通过拟合某些训练点来补偿被禁止的交互;这是解释性判断,并非已经证明的机制。尖峰的可见数量也受 PDP 网格分辨率影响,不能单凭一张图确认“模型没有噪声”或“有多少转折点”。

二维 PDP 展示成对特征的联合响应

_, ax = plt.subplots(ncols=3, figsize=(10, 4), constrained_layout=True)
PartialDependenceDisplay.from_estimator(
    hgbdt_model,
    X_train,
    features=["temp", "humidity", ("temp", "humidity")],
    kind="average",
    ax=ax,
    **common_params,
)
官方梯度提升模型的温度、湿度一维PDP与二者的二维联合PDP,二维图显示约20摄氏度附近的变化形态随湿度改变。
原图:The scikit-learn developers,BSD-3-Clause;1-way vs 2-way of numerical PDP using gradient boosting。图中数值为原文模型输出。

二维图把温度和湿度同时设为网格值,再平均其他特征。原文观察到,在约 20°C 以上,湿度的影响看起来对温度较不敏感;较低温度时,两者都持续改变模型响应。20°C 附近的“脊”在干燥条件下较陡,在高于约 70% 的湿度下较平缓。

把相同绘图调用中的模型换成 hgbdt_model_without_interactions,可以得到禁止交互的对照。其二维图可能因单特征的高频尖峰形成网格状噪声,肉眼判断交互更困难,但前述随温度阈值变化的简单交互形状不再明显。不要把二维图的每处纹理都认作交互。

类别变量也能做二维 PDP。将关注特征改为 ["season", "weather", ("season", "weather")],同时传入 categorical_features=categorical_features,就会显示两个类别的一维柱状图及成对类别的离散热图。类别没有连续数值间距,不应把其编码整数画成连续梯度并进行插值解释。

需要三维视角时,使用同一份部分依赖结果

import numpy as np
from sklearn.inspection import partial_dependence

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
fig = plt.figure(figsize=(5.5, 5))
ax = fig.add_subplot(projection="3d")
surface = ax.plot_surface(XX, YY, Z, cmap=plt.cm.BuPu, edgecolor="k")
ax.set_xlabel(features[0])
ax.set_ylabel(features[1])
ax.view_init(elev=22, azim=122)
fig.colorbar(surface, ax=ax, pad=0.08, shrink=0.6, aspect=10)

这是原文三维绘图段的整理版:保留网格与转置关系,去掉针对 Matplotlib 3.2 以前版本的兼容性导入,以及重复添加已存在轴的语句;为 colorbar 显式指定轴。三维曲面只是同一类数值的另一种呈现,不增加新的统计证据。看不清的遮挡区域应回到二维等高或热图检查。

自定义网格,比较同一组输入点

默认网格根据输入数据的分位数生成。若要比较训练数据略有不同的模型,或者有意观察分布外行为,可通过 custom_values 指定检查点。它会覆盖对应特征的 grid_resolution 和 percentiles,其他特征仍依照默认规则。

_, ax = plt.subplots(ncols=2, figsize=(6, 4), sharey=True, constrained_layout=True)
PartialDependenceDisplay.from_estimator(
    hgbdt_model,
    X_train,
    features=["temp", "humidity"],
    kind="both",
    ax=ax,
    custom_values={"temp": np.linspace(0, 40, 10)},
    **common_params,
)

这个参数把温度网格设为 0–40 之间的十个值。网格包含一个值,并不代表训练数据在该区域提供了充分支持。记录参考样本、网格范围、抽样种子和模型版本,才能让两张图之间的差异可解释、可核查。

代码、数据与解释的边界

此次只做静态审查:fetch_openml 会联网下载并写入缓存,训练与并行部分依赖会消耗 CPU、内存和时间;本文未执行。示例没有发现硬编码秘密或拼接命令注入,但也未进行安全测试。实际任务需固定依赖与数据版本,验证缺失值和未知类别,控制对外下载源与本地输出路径。这里保留的计算代码并不把结果自动写入生产系统。

统计解释还需同时检查预测质量、训练/测试时间分布差异、相关特征的支持区域、ICE 异质性和网格敏感性。单个平均曲线不能取代个体曲线,漂亮的曲面也不能取代独立测试或因果设计。

原文参考:Christoph Molnar,Interpretable Machine Learning(2019);Goldstein、Kapelner、Bleich 与 Pitkin,Peeking Inside the Black Box: Visualizing Statistical Learning With Plots of Individual Conditional Expectation,Journal of Computational and Graphical Statistics 24(1), 44–65(2015)。

Copyright (c) 2007–2026 The scikit-learn developers。代码与原图保留 BSD-3-Clause 归属和许可,完整许可证见 sources/LICENSE.txt。中文译编与审查注:未完纪编辑;重复绘图代码适度合并、补充统计与版本边界,并标明三维段改动。委托方于 2026-10-05 确认全文翻译、转载及配图授权。本文未做运行测试;未发现问题不代表不存在漏洞。

保留的完整许可声明

以下为原项目适用许可文本,原作者与文档归属及本文改动说明见正文。

BSD 3-Clause License

Copyright (c) 2007-2026 The scikit-learn developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

* Neither the name of the copyright holder nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发

请登录后发表评论

    暂无评论内容