Unlocking Black‑Box Models: A Practical Guide to PDP, ICE, and Post‑Hoc Interpretation

Post‑hoc Interpretation Methods for Complex Models

Although intrinsically interpretable machine‑learning models offer good explainability, they still cannot match the predictive accuracy of complex models. When both high accuracy and interpretability are required, post‑hoc interpretation methods can bridge the gap. Common techniques include sensitivity analysis, Partial Dependence Plot (PDP), Accumulated Local Effects (ALE), Shapley‑based fair contribution, and recent surrogate‑model approaches.

When a black‑box model must be used, scientists apply various post‑hoc methods to open the box and help users understand its internal mechanisms. LIME (Local Interpretable Model‑agnostic Explanation) was introduced by Ribeiro, Singh and Guestrin in 2016, using surrogate models to explain individual predictions. SHAP (Shapley‑Additive Explanations) was proposed by Lundberg and Lee in 2017, leveraging Shapley values to attribute feature contributions.

Partial Dependence Plot

Linear regression models struggle with complex data, while modern models such as SVM, random forest, and gradient‑boosting achieve higher accuracy but sacrifice interpretability. PDP, introduced by Friedman in 2001, visualizes the marginal effect of one or two target features on model predictions, providing an intuitive way to diagnose model behavior.

Partial Dependence Function

Consider a credit‑scoring task where a model predicts the probability of default. Business users need indicators that clearly separate high‑risk from low‑risk customers. The PDP shows how one or two target features marginally affect the prediction, indicating the importance of those features.

Formally, the partial dependence function is defined as the expectation of the model output over the distribution of non‑target features while fixing the target feature(s):

<code>f_{PD}(x_T) = \mathbb{E}_{X_{\neg T}}[f(x_T, X_{\neg T})]</code>

Here, f is the trained model, x_T denotes the target feature(s), and X_{\neg T} represents all other features.

Estimation Method

Directly computing the expectation is often infeasible because many combinations of x_T do not appear in the finite dataset. In practice, the marginal expectation is approximated by averaging the model predictions over the observed data:

<code>\hat{f}_{PD}(x_T) = \frac{1}{n}\sum_{i=1}^{n} f(x_T, x_{i,\neg T})</code>

where n is the number of samples and x_{i,\neg T} are the non‑target feature values of the i‑th sample.

For example, to study the effect of weight on food intake, we fix weight at a grid of values, replace all other features (e.g., height, age) with their average values, and compute the model prediction for each grid point. The resulting curve is the PDP.

Limitations of Partial Dependence Plots

PDPs can display at most two target features simultaneously; higher‑dimensional relationships become impractical to visualize. Moreover, the method assumes independence between target and non‑target features. When this assumption is violated, the estimated marginal effect may be biased because the integration includes unrealistic feature combinations.

Individual Conditional Expectation (ICE) Plots

While PDP provides a global, averaged view, ICE (Individual Conditional Expectation) visualizes the dependence of the prediction on a target feature for each individual sample. ICE curves are generated by applying the same estimation procedure to each observation, and the average of all ICE curves equals the PDP.

The following figure illustrates an ICE plot:

Example Demonstration

We use a small dataset (≈500 samples, 13 numeric features) to predict house prices with XGBoost. The target variable is continuous, so we employ XGBRegressor . The workflow includes data loading, train‑test split, model training, and then generating PDP and ICE visualizations for the feature NOX .

Code Implementation

<code>from sklearn.datasets import load_boston
import pandas as pd
from sklearn.model_selection import train_test_split
import xgboost as xgb

data = load_boston()
df  = pd.DataFrame(data['data'], columns=data['feature_names'])
df['target'] = data['target']
X, y = df.iloc[:, :-1], df['target']
X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.7)

xgb_reg = xgb.XGBRFRegressor(random_state=1)
xgb_reg.fit(X_train, y_train)

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn

XC_name = X.drop('NOX', axis=1).columns
XC_mean = X[XC_name].mean().values
XS_min, XS_max = X['NOX'].min(), X['NOX'].max()
XS_grid = np.arange(XS_min, XS_max, (XS_max - XS_min) / 1000)

pdp_data = pd.DataFrame(np.tile(XC_mean, (1000, 1)))
pdp_data.columns = XC_name
pdp_data.insert(loc=4, column='NOX', value=XS_grid)
</code>

Using the trained model, we predict on pdp_data and plot the PDP for NOX . The resulting curve shows a rapid increase in price when NOX is between 0.6 and 0.7, followed by a decrease.

To explore individual sample behavior, we generate ICE curves:

<code>fig = plt.figure()
for i in range(50):
    XC_i = X[XC_name].iloc[i, :].values
    ICE_data = pd.DataFrame(np.tile(XC_i, (1000, 1)))
    ICE_data.columns = XC_name
    ICE_data.insert(loc=4, column='NOX', value=XS_grid)
    preds = xgb_reg.predict(ICE_data)
    sn.lineplot(x=XS_grid, y=preds)
plt.xlabel('NOX')
plt.ylabel('Price')
plt.show()
</code>

References

Shaw Ping, Yang Jianying, Su Sida. Interpretable Machine Learning: Models, Methods, and Practices .