Partial dependence plots and what they quietly assume

A churn model arrives at the review meeting with one plot: partial dependence of predicted churn on tenure_days. The curve drops steeply across the first year, then flattens out. Someone says the obvious thing — “so get people past twelve months and churn falls; discount the first year.” Nobody objects, because the plot looks exactly like a picture of what would happen. It is not one.

What the curve is made of

A partial dependence plot — Friedman’s device, from the 2001 gradient boosting paper — answers a mechanical question. Take one feature. Take a grid of values across its range. For each grid value, overwrite that column in every row of the dataset, score all the rows, and average the predictions. Plot those averages against the grid.

That is the entire definition — the same one sklearn.inspection.partial_dependence implements — and it fits in six lines without a library:

The model below is a churn model on the public dataset, and the feature is the subscription’s monthly price.

import numpy as np
import pandas as pd
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import train_test_split

subs = pd.read_csv("https://dataacademy.ai/data/subscriptions.csv",
                   parse_dates=["signed_up_at", "cancelled_at"])
events = pd.read_csv("https://dataacademy.ai/data/events.csv",
                     parse_dates=["occurred_at"])

ev = events.merge(subs[["subscription_id", "customer_id", "signed_up_at"]],
                  on="customer_id")
month1 = ev[(ev.occurred_at - ev.signed_up_at).dt.days.between(0, 29)]
counts = month1.pivot_table(index="subscription_id", columns="event_type",
                            values="event_id", aggfunc="count", fill_value=0)

d = subs.set_index("subscription_id").join(counts)
d[counts.columns] = d[counts.columns].fillna(0)
y = d.cancelled_at.notna().astype(int)
X = d[["monthly_price", "view", "login", "export", "purchase"]].astype(float)
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.3,
                                                  random_state=0, stratify=y)
model = GradientBoostingClassifier(random_state=0).fit(X_train, y_train)

grid = np.linspace(X_val["monthly_price"].min(), X_val["monthly_price"].max(), 20)
curve = []
for v in grid:
    Xg = X_val.copy()
    Xg["monthly_price"] = v
    curve.append(model.predict_proba(Xg)[:, 1].mean())

print(pd.DataFrame({"price": grid.round(1), "pdp": np.round(curve, 3),
                    "rows_at_this_price":
                        [(X_val.monthly_price == v).sum() for v in grid]}))

    price    pdp  rows_at_this_price
0    12.0  0.585                 145
1    24.5  0.534                   0
2    36.9  0.534                   0
3    49.4  0.534                   0
4    61.9  0.381                   0
5    74.4  0.381                   0
...
17  224.1  0.363                   0
18  236.5  0.363                   0
19  249.0  0.363                  38

The .astype(float) is not cosmetic: scikit-learn refuses to compute partial dependence on integer columns, because the grid values would be silently rounded back onto the integers.

PartialDependenceDisplay does the same thing with better grids and a plot attached:

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

fig, ax = plt.subplots(figsize=(9, 4))
PartialDependenceDisplay.from_estimator(
    model,
    X_val,
    features=["monthly_price", "view"],
    kind="average",
    grid_resolution=20,
    ax=ax,
)

Two things follow straight away. The curve describes the model, not the world — a badly fitted model produces a clean, confident, meaningless plot. And every point on it is an average over rows that were forced to a value, whether or not those rows could ever hold it.

The average runs through rows that cannot exist

Look at the third column of that table. This product has four prices — 12, 29, 79 and 249 — and a twenty-point linear grid lands on two of them. The other eighteen rows of the curve are averages over a validation set in which every subscription has been repriced to a number nobody is charged, and 29 and 79, which between them cover most of the customer base, never appear on the plot at all. The curve is real output, computed correctly, describing a price list that does not exist.

The same failure is less visible when the feature is continuous. Fix tenure_days at 30 and the copied dataset now contains customers with a month of history, forty lifetime orders and a loyalty tier that takes three years to reach. Fix age at 25 and it contains a 25-year-old with 30 years of service. These rows get scored like any other, and their predictions go into the average that becomes the height of the curve at that grid point.

The model has no training data anywhere near them, so what it returns there is whatever its shape happens to extend to: gradient boosting holds the value of the last split it saw, a linear model keeps going in a straight line, a neural network does something unadvertised. None of that is evidence about anything.

How much of the curve is fabricated depends entirely on how correlated the feature is with the rest of the row. For a feature that is close to independent, almost every synthetic row is plausible and the PDP is fine. For two features correlated at 0.9, most of the joint grid is empty space, and the plot spends most of its width in territory the model was never fitted on. This is the same objection that makes permutation importance unreliable under correlation — both methods break the joint distribution and then trust what comes back.

Check the region before trusting the curve. A rug plot of the feature, or a scatter of the feature against its strongest correlate, shows where the data actually lives. Read the PDP only across that stretch. Print the support beside the curve and the stretch names itself:

for v in np.linspace(X_val["view"].min(), X_val["view"].max(), 20):
    Xg = X_val.copy()
    Xg["view"] = v
    print(f"{v:5.1f} {model.predict_proba(Xg)[:, 1].mean():.3f} "
          f"{(X_val['view'] >= v).sum():4d}")

  0.0 0.614  360
  3.9 0.562  237
  7.8 0.599  139
 11.7 0.601   87
...
 38.9 0.449    7
 42.8 0.449    5
 46.7 0.662    4
 50.6 0.665    2

The jump from 0.449 to 0.662 at the right-hand end is not a finding about heavy users. It is four rows, then two, and the model has nothing to go on out there. The left half of that curve is worth reading; the right half is the model’s extrapolation shape drawn as if it were evidence.

ICE curves show the spread the average hides

An individual conditional expectation curve is the same computation without the final averaging step: one line per row, showing how that row’s prediction moves as the feature sweeps the grid. The PDP is the mean of the ICE curves, so everything the mean throws away is visible in the bundle.

What it throws away is often the finding. Suppose a price feature pushes predicted churn up for month-to-month customers and down for annual contracts. Both effects are real, both are large, and the PDP is a flat line. Read alone, that flat line says the model ignores price. The ICE plot says the model has learnt two opposite behaviours and the segments cancel.

fig, ax = plt.subplots(figsize=(9, 4))
PartialDependenceDisplay.from_estimator(
    model,
    X_val,
    features=["monthly_price"],
    kind="both",          # ICE lines plus the average
    subsample=200,        # 200 rows, otherwise it is a smear
    centered=True,        # start every curve at zero
    random_state=0,
    ax=ax,
)

centered=True matters more than it sounds. Uncentred, the lines are stacked by their base prediction and you compare vertical positions instead of shapes. Centred, every curve starts at zero and the plot shows only the change, which is the thing under discussion. Fanning lines mean the effect depends on the rest of the row; roughly parallel lines mean it does not, and only then does the average stand in for the individuals.

The honest alternative when features are correlated

Accumulated local effects fix the impossible-row problem by never leaving the data. Split the feature into intervals. Inside each interval, take only the rows that genuinely fall there, score each of them at the interval’s lower and upper edge, and average the difference. Then accumulate those local differences across intervals to build the curve. Every prediction used comes from a row that is already in the right neighbourhood, so nothing is extrapolated. The cost is that the result is a relative effect — the level of the curve is arbitrary, only the shape is meaningful — and narrow intervals with few rows get noisy. It is not in scikit-learn; alibi and PyALE implement it. When ALE and PDP tell the same story, the correlation was not doing damage. When they disagree, believe ALE.

The sentence that gets people in trouble

Everything above is a technical caveat. This one is the failure that reaches production, and it survives all the caveats being fixed.

A PDP is not a policy simulator. The curve says how the model’s output varies across a grid, given the data as it stands. It does not say what happens if you change the feature, because changing a feature in the world changes the others along with it — the same trap as reading a SHAP value as a lever. Discounting the first year does not convert three-month customers into fifteen-month customers. It creates discounted three-month customers, a population the model has never scored, whose other columns now move too.

And tenure carries survivorship in it. Customers with long tenure are, by construction, the ones who did not leave — people who liked the product, who picked the right plan, whose use case fitted. The curve slopes down partly because staying causes long tenure, not the other way round. No amount of gridding will separate those.

So read the plot with the model named out loud. “Across this range, the model’s predicted churn falls with tenure” is defensible and often useful — it tells you what the model learnt, and whether that matches what the business believes. Drop those two words, and the same picture becomes a promise about next quarter that nothing in the computation ever supported.