Partial Dependence Plots and ICE Curves: See How Your Model Really Uses Each Feature

1. The Problem: Your Model Works, But Why?

Last time, the team compared LIME and SHAP to explain Maria’s individual denial — weighing speed against consistency. Mr. Owens saw both tools explain why the model said no to Maria specifically, but then he tapped the Income column and asked a broader question: “I can see why the model said no to Maria. But how does income, in general, affect approval odds across all our applicants? Does the relationship go up? Does it flatten out? I want to see the shape of it, not just one person’s number.”

You could show him a feature importance bar chart. Income ranks highest — great. But that still doesn’t answer his question. Does approval probability climb linearly with income? Level off after some threshold? Could the model actually treat an extremely high income as a risk signal that hurts someone’s odds?

Feature importance tells you which features matter, not how they matter. To see inside the black box, we need a way to visualize the relationship between a single feature and the model’s prediction. Partial Dependence Plots (PDP) give us that. They let us ask: if I change this one feature, holding everything else constant, how does the prediction change?

2. Partial Dependence Plots: The Core Idea

Think of a Partial Dependence Plot (PDP) as a way to surface the average behavior of a single feature. Say you have 1,000 rows of data. To see how ‘Income’ affects ‘Home Price,’ you run a small experiment:

1. Pick a specific income value — say $50,000.
2. Walk through every row in the dataset and set income to $50,000, leaving location, age, and the rest untouched.
3. Ask the model to predict the home price for each modified row.
4. Average those predictions. That single number is one point on the plot.
5. Repeat for $60,000, $70,000, and so on.

Connect the points and you get a curve. This is the marginal effect of income on the predicted outcome — the model’s logic collapsed to one dimension.

The partial dependence of the model’s prediction $f$ on a subset of features $x_S$ is defined (Friedman, 2001) as:

$$\bar{f}_S(x_S) = \mathbb{E}_{x_C}\!\left[f(x_S, x_C)\right] = \int f(x_S, x_C)\, p_C(x_C)\, dx_C$$

In plain English: hold the features of interest ($x_S$) at fixed values, average the model’s prediction over the observed distribution of all other features ($x_C$), and that average is one point on the PDP curve. In practice we can’t integrate analytically, so we estimate it by Monte Carlo — just average over the data we already have:

$$\hat{\bar{f}}_S(x_S) = \frac{1}{n} \sum_{i=1}^{n} f\!\left(x_S,\; x_C^{(i)}\right)$$

An ICE curve for a single row $i$ is simply $f\!\left(x_S,\; x_C^{(i)}\right)$ evaluated at each grid value of $x_S$ — without the averaging step.

Plain English	Statistical symbol	Python equivalent
The feature(s) we’re sweeping (held at grid values)	$x_S$ (the subset $S$)	features=['MedInc']
All other features (kept at their original values)	$x_C$ (the complement set)	Every column in X except MedInc
The trained model’s prediction function	$f$	model.predict
Partial dependence (average prediction at a given $x_S$)	$\bar{f}_S(x_S)$	The y-value on the PDP curve
Distribution of the “other” features	$p_C(x_C)$	Approximated by the empirical distribution of X
Monte Carlo estimate (average over $n$ data rows)	$\frac{1}{n}\sum_{i=1}^{n} f(x_S, x_C^{(i)})$	Computed internally by PartialDependenceDisplay.from_estimator
One row’s prediction as $x_S$ varies (an ICE curve)	$f(x_S, x_C^{(i)})$ for fixed $i$	One line in the kind='both' plot

Key distinction: The PDP averages over $i$; the ICE curve does not. That’s why the PDP is one smooth line while ICE is a tangle of individual trajectories.

3. Code: Computing and Plotting Partial Dependence

Let’s apply this to the California Housing dataset. We’ll train a Random Forest—a classic “black box” model—then use Scikit-Learn to visualize what it learned about income.

import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import RandomForestRegressor
from sklearn.inspection import PartialDependenceDisplay

# 1. Load data
data = fetch_california_housing()
X, y = data.data, data.target
feature_names = data.feature_names

# 2. Train a black-box model
model = RandomForestRegressor(n_estimators=50, random_state=42)
model.fit(X, y)

# 3. Compute and plot PDP for 'MedInc' (Median Income)
fig, ax = plt.subplots(figsize=(8, 6))
display = PartialDependenceDisplay.from_estimator(
    model, X, features=['MedInc'], 
    feature_names=feature_names, 
    ax=ax
)
plt.title("Partial Dependence of House Value on Median Income")
plt.show()

• fetch_california_housing() — Loads the California Housing dataset bundled with scikit-learn. Returns a Bunch object with .data (the feature matrix, ~20,600 rows × 8 features), .target (median house values), and .feature_names (column labels like 'MedInc', 'HouseAge', 'AveRooms', etc.).
• RandomForestRegressor(n_estimators=50, random_state=42) — Creates (but does not yet train) a random forest regressor with 50 decision trees. random_state=42 seeds the bootstrapping so the forest is identical on every run.
• model.fit(X, y) — Trains the forest on the full dataset (no train/test split here because we’re interested in interpreting the model’s learned logic, not evaluating generalisation).
• PartialDependenceDisplay.from_estimator(model, X, features=['MedInc'], ...) — The workhorse call. Internally it builds a grid of MedInc values (e.g., 0.5, 1.0, 1.5, … up to ~15), replaces the MedInc column in every row of X with each grid value, asks the model to predict, and averages the predictions across all rows. The resulting curve is rendered on the supplied axis.
• feature_names=feature_names — Passes the column names so the plot’s x-axis is labelled 'MedInc' instead of 'Feature 0'.
• ax=ax — Directs the plot to the Matplotlib axis we created with plt.subplots, so it renders inside our figure rather than creating a new one.
• display = ... — The return value (a PartialDependenceDisplay object) is stored in display but never used afterwards; the plot is already rendered on ax, so the variable isn’t needed. It could be used to call display.plot() again later if desired.

4. What the Curve Actually Means: Interpreting Partial Dependence

Look at the plot we just generated. Median Income sits on the x-axis; the model’s predicted house value sits on the y-axis.

The curve probably climbs steeply at first, then flattens. So as income rises from 2 to 6, predicted price jumps significantly. Past 8, extra income barely moves the prediction.

• A steep upward curve means the feature has a strong positive effect.
• A flat line means the feature doesn’t change the prediction at all (on average).
• A downward curve means the feature has a negative effect.

5. The Catch: Partial Dependence Assumes Features Are Independent

Here’s the tricky part about PDP: it assumes your features don’t move together. In practice, they often do. In our housing data, ‘Average Rooms’ and ‘Average Bedrooms’ are tightly linked.

When PDP pushes ‘Average Rooms’ high while holding ‘Average Bedrooms’ low, it can create a “Frankenstein” house—a 10-room mansion with 0.5 bedrooms. The model has never seen such a house. Its prediction could be nonsense. Before trusting a PDP, always check whether your features are heavily correlated.

import pandas as pd
# Check correlation between rooms and bedrooms
df = pd.DataFrame(X, columns=feature_names)
print(f"Correlation: {df['AveRooms'].corr(df['AveBedrms']):.2f}")
# A correlation of 0.85 means these features are tightly linked!

• pd.DataFrame(X, columns=feature_names) — Wraps the NumPy array X in a pandas DataFrame with human-readable column names from feature_names. This lets us select columns by name (e.g., df['AveRooms']) instead of by integer index.
• df['AveRooms'].corr(df['AveBedrms']) — Computes the Pearson correlation coefficient between the AveRooms and AveBedrms columns. Returns a float between −1 and 1: 1 means perfectly positively correlated, 0 means no linear relationship, −1 means perfectly negatively correlated.
• :.2f — Formats the correlation to two decimal places (e.g., 0.85).
• The output (≈ 0.85) confirms that these two features move together tightly — a red flag for PDP, because the algorithm will create unrealistic combinations of the two when it sweeps one while holding the other constant.

6. Individual Conditional Expectation (ICE) Curves: One Row at a Time

PDP shows the average effect. But averages can lie. Picture a feature that helps half your rows and hurts the other half. The average flattens to a straight line, and the feature looks inert.

Individual Conditional Expectation (ICE) curves address this. Instead of one averaged line, an ICE plot draws a line for every single row in your data. You see how the prediction for one specific house moves as you sweep the feature value.

7. Code: Computing and Plotting ICE Curves

Same Scikit-Learn function, just a different parameter. Switch kind to 'both' and you get the individual lines and the average PDP line together in one view.

fig, ax = plt.subplots(figsize=(8, 6))
# 'kind="both"' plots both individual (ICE) and average (PDP) lines
PartialDependenceDisplay.from_estimator(
    model, X, features=['MedInc'], 
    feature_names=feature_names, 
    kind='both', 
    subsample=50, # Plot 50 individual lines to keep it clean
    ax=ax
)
plt.title("ICE Curves and PDP for Median Income")
plt.show()

• kind='both' — Tells scikit-learn to overlay both the average PDP line (usually shown in a distinct colour or thicker stroke) and the individual ICE lines (one per row) on the same axes. Other valid values: 'individual' (ICE lines only) or 'average' (PDP only, the default).
• subsample=50 — Randomly selects 50 rows from X to draw ICE lines for, rather than all ~20,600. Without subsampling, the plot would be an unreadable tangle of thousands of overlapping lines; 50 keeps it visually tractable while still showing the spread.
• features=['MedInc'] — Same as before: the feature whose effect we’re visualising. The ICE lines show how each individual row’s prediction changes as MedInc is swept from low to high, while the row’s other features stay at their original values.
• ax=ax — Renders the plot on our pre-created Matplotlib axis.

8. ICE Curves Reveal Interactions: When Features Work Together

Here’s the payoff. Parallel ICE lines mean the feature affects every row the same way. When lines cross or fan out, you’ve found an interaction.

For example, ‘Income’ might have a huge effect on house price in San Francisco but barely matter in a rural area. The city ICE curves would be steep; the rural ones flat. The PDP just averages across both and hides the difference. ICE curves show you this heterogeneity (a fancy word for “differences”) that the PDP misses.

9. Comparing PDP and ICE: When to Use Each

So what does that mean in practice?

• Use PDP when you need to explain the general trend to a stakeholder. It’s clean and easy to read.
• Use ICE when you’re debugging the model. It shows if the model behaves inconsistently or if there are hidden interactions you missed.

I almost always plot them together. The PDP gives you the main story. The spread of the ICE curves tells you how much to trust that story for any individual case.

PDP / ICE vs. SHAP Dependence Plots

Both PDP/ICE and SHAP dependence plots visualise how a single feature relates to the model’s output — but they answer subtly different questions and come with different trade-offs.

Aspect	PDP / ICE Curves	SHAP Dependence Plots
What it shows	Model’s predicted value (y-axis) vs. feature value (x-axis)	SHAP value for that feature (y-axis) vs. feature value (x-axis)
Question it answers	“How does the prediction change as this feature changes?”	“How much did this feature push each individual prediction?”
Granularity	PDP = average; ICE = per-row lines	Every dot is one prediction — no averaging
Interactions	ICE fanning/crossing hints at interactions, but doesn’t identify which feature	Automatically colours dots by the most interacting second feature
Correlated features	Creates unrealistic “Frankenstein” rows (sweep one, hold the other) — can mislead	Each dot is a real prediction from real data — no synthetic combinations
Speed	Fast; built into scikit-learn, no extra library	Requires computing SHAP values first (fast with TreeExplainer, slower with KernelExplainer)
Stakeholder readability	High — a simple line graph is intuitive	Medium — scatter plots with colour gradients need explanation

When to reach for PDP / ICE

• ✅ You want to see the overall shape of a feature’s effect (linear, saturating, non-monotonic) — the “how” question.
• ✅ You’re explaining to a non-technical stakeholder who wants a clean, single-line graph.
• ✅ You want to check for heterogeneity (ICE) — do different rows respond differently to the same feature?
• ✅ You don’t have SHAP values computed and want a quick visual without extra dependencies.

When to reach for SHAP dependence plots

• ✅ You’ve already computed SHAP values and want to visualise them.
• ✅ You want to automatically surface interactions — the colour encoding pinpoints the second feature that matters most.
• ✅ You’re worried about correlated features and want to avoid the “Frankenstein” problem entirely.
• ✅ You want to connect the global trend to individual predictions — each dot ties back to a specific row’s SHAP decomposition.

The short version

Use PDP/ICE when you want to see the shape of a feature’s effect (the “how does it curve?” question). Use SHAP dependence plots when you want to see how a feature drives individual predictions and spot interactions (the “who and why” question). They’re complementary — PDP shows the forest, SHAP dependence shows both the forest and the trees, with each tree labelled by what it is.

10. Common Pitfalls and How to Avoid Them

1. Extrapolation: Don’t trust the ends of the plot. If your data only has incomes between 2 and 10, the model’s prediction for an income of 50 is a guess — nothing more.
2. The Flat Curve Trap: A flat PDP doesn’t always mean a feature is useless. It might matter a great deal, but only when combined with another feature. That’s an interaction, and ICE curves help you spot it.
3. Correlation: As mentioned, when features are highly correlated, the “hold everything else constant” logic breaks down.

11. Putting It Together: A Real-World Example

A final example: ‘House Age’ and ‘Average Occupancy’.

fig, ax = plt.subplots(1, 2, figsize=(12, 5))
PartialDependenceDisplay.from_estimator(model, X, ['HouseAge'], feature_names=feature_names, kind='both', subsample=50, ax=ax[0])
PartialDependenceDisplay.from_estimator(model, X, ['AveOccup'], feature_names=feature_names, kind='both', subsample=50, ax=ax[1])
plt.tight_layout()
plt.show()

• plt.subplots(1, 2, figsize=(12, 5)) — Creates a 1×2 grid of subplots (two panels side by side). fig is the overall figure; ax is now an array of two axis objects (ax[0] and ax[1]).
• ['HouseAge'] and ['AveOccup'] — Each call specifies a different feature to compute the PDP/ICE for, passed positionally (the third positional argument is features) rather than as features=.
• feature_names=feature_names — Still required here, and easy to forget once the feature list is passed positionally. X is a plain NumPy array, not a DataFrame, so scikit-learn has no column names to fall back on — drop this argument and the call raises ValueError: Feature 'HouseAge' not in feature_names, because there’s nothing for the string 'HouseAge' to resolve against.
• ax=ax[0] / ax=ax[1] — Directs the first PDP/ICE plot to the left panel and the second to the right panel.
• plt.tight_layout() — Adjusts spacing between the two subplots to prevent axis labels and titles from overlapping. Always call this before plt.show() when rendering multiple panels.

The ‘HouseAge’ lines stay mostly flat. On average, the model doesn’t lean much on house age. The ‘AveOccup’ plot shows a sharp drop-off instead. As occupancy rises, the predicted value falls—high-occupancy homes tend to be smaller, crowded rentals, which likely explains the drop.

12. Next Steps: From Interpretation to Explanation

You can now see the “shape” of your model’s logic. Not just that a feature matters, but how it shifts the outcome.

• PDP gives you the average forest view.
• ICE gives you the individual tree view.

So the team has three tools, each answering a different question. SHAP values explain a single prediction. LIME gives a fast approximation. PDP and ICE curves show the overall shape of a feature’s effect. Mr. Owens looks at the growing collection of plots and says: “This is thorough — but I need one explanation I can actually give Maria and the regulator. Not three different charts. One story.” That’s what the final article tackles: turning these tools into one human-readable explanation.

Check Your Understanding

Questions below move from simple recall up to open-ended design, roughly following Bloom’s Taxonomy.

Remember What are the five steps the article describes for manually computing one point on a Partial Dependence Plot?

Understand In your own words, explain the “half the people helped, half hurt” scenario the article uses to show why a flat PDP can be misleading, and how ICE curves reveal what PDP hides in that case.

Apply Using the article’s interpretation rules (steep upward = strong positive effect, flat = no effect, downward = negative effect), what would you conclude about a feature whose PDP curve rises steeply from x=0 to x=5 and then stays completely flat from x=5 to x=10?

Analyze The article explains that PDP creates “Frankenstein” rows when features are correlated (a 10-room house with 0.5 bedrooms). Walk through why this specific problem doesn’t occur when you plot the ICE curve for a single row varying only one feature — is the correlation problem eliminated, or does it still apply to the “modified” values within that one row’s sweep?

Evaluate The article recommends checking feature correlation before trusting a PDP (using .corr()) but doesn’t specify a threshold for when correlation is “too high to trust.” Critique the article’s example (AveRooms/AveBedrms at 0.85) — what would you actually do differently with the PDP interpretation once you know two features are that correlated, beyond just noting the number?

Create Design a PDP/ICE investigation for a new scenario: a hospital readmission model where you suspect medication_count matters differently for patients under 65 versus over 65. Describe how you’d use ICE curves (not just the average PDP) to test this specific hypothesis, and what pattern in the ICE lines would confirm it.

---

Related articles

• SHAP Values: Why Your Model Predicted That
• LIME vs SHAP: Choosing the Right Translator for You)
• Why Did the Model Say No? Explaining Black-Box Decisions to the Business)

---

References & Further reading

• Friedman, J. H. (2001). “Greedy Function Approximation: A Gradient Boosting Machine.” Annals of Statistics, 29(5), 1189–1232. https://projecteuclid.org/euclid.aos/1013203451 — The foundational paper that introduced Partial Dependence Plots as part of the Gradient Boosting Machine framework. The PDP definition and Monte Carlo estimation procedure used by sklearn.inspection originate here.
• scikit-learn inspection module documentation. https://scikit-learn.org/stable/modules/partial_dependence.html — Official API reference for PartialDependenceDisplay, partial_dependence, and the kind parameter that switches between PDP, ICE, and combined plots.
• Kaggle: “House Prices: Advanced Regression Techniques” competition. https://www.kaggle.com/c/house-prices-advanced-regression-techniques — A real-world housing dataset with 79 features (square footage, neighbourhood, year built, etc.) where PDP and ICE curves are widely used in top-scoring notebooks to visualise how each feature drives predicted sale price — an excellent sandbox for practising the techniques in this article.

---