Partial dependence as data

The partial_effects() method produces Altair charts that are ready to display, but sometimes you need the underlying numbers. You might want to plot with matplotlib, export the curves to a CSV, run a custom analysis on the estimated effects, or combine smooth estimates across models. The partial_dependence() method returns the same data that partial_effects() plots, but as structured arrays rather than chart objects.

Why partial dependence as data?

partial_effects() is the fast path: one call gives you a polished visualization of every smooth term. But Altair charts are not always what you need:

  • Matplotlib or seaborn workflows: your project already uses matplotlib and you want a consistent style.
  • Export: you want to write the estimated curves to a file for a collaborator or a table in a paper.
  • Custom analysis: you want to compute the area under a smooth, find the x-value where the effect crosses zero, or compare effects across models numerically.
  • Fine-grained control: you need to overlay observed data, add reference lines, or combine panels in ways that go beyond what the built-in plots offer.

partial_dependence() gives you structured PartialDependenceResult objects with arrays you can manipulate freely.

Basic usage

Fit a model on the motorcycle crash dataset and call partial_dependence():

import numpy as np
import whittaker as wk

data = wk.load_dataset("mcycle")

model = wk.GAM("accel ~ s(times, k=15)")
model.fit(data, method="REML")

results = model.partial_dependence()

The return value is a list of PartialDependenceResult objects, one per smooth term in formula order:

print(f"Number of results: {len(results)}")
print(results[0])
Number of results: 1
PartialDependenceResult(term='s(times, k=15)', n_grid=200, edf=9.7, level=0.95)

Each result is a dataclass with the following fields:

r = results[0]

print(f"term:   {r.term}")
print(f"x keys: {list(r.x.keys())}")
print(f"effect: shape {r.effect.shape}, dtype {r.effect.dtype}")
print(f"se:     shape {r.se.shape}")
print(f"lower:  shape {r.lower.shape}")
print(f"upper:  shape {r.upper.shape}")
print(f"edf:    {r.edf:.2f}")
print(f"level:  {r.level}")
print(f"n_grid: {r.n_grid}")
term:   s(times, k=15)
x keys: ['times']
effect: shape (200,), dtype float64
se:     shape (200,)
lower:  shape (200,)
upper:  shape (200,)
edf:    9.71
level:  0.95
n_grid: 200

The x dictionary maps covariate names to their evaluation grids. For a 1-D smooth like s(times), there is a single key. The effect, se, lower, and upper arrays all have the same length as the grid.

Plotting with the raw arrays

With the arrays in hand, building a custom plot is straightforward. Here we plot the estimated effect as a line with a shaded confidence band using Altair:

import altair as alt

r = results[0]
x_vals = r.x["times"]

plot_data = [
    {"times": float(x_vals[i]), "effect": float(r.effect[i]),
     "lower": float(r.lower[i]), "upper": float(r.upper[i])}
    for i in range(len(x_vals))
]

band = alt.Chart({"values": plot_data}).mark_area(
    opacity=0.25, color="steelblue"
).encode(x="times:Q", y="lower:Q", y2="upper:Q")

line = alt.Chart({"values": plot_data}).mark_line(
    color="steelblue", strokeWidth=2
).encode(x=alt.X("times:Q", title="times"), y=alt.Y("effect:Q", title=f"Effect ({r.term})"))

ref = alt.Chart({"values": [{"y": 0}]}).mark_rule(
    color="gray", strokeDash=[4, 4]
).encode(y="y:Q")

(band + ref + line).properties(
    width="container", height=300,
    title=f"Partial dependence: {r.term} (EDF = {r.edf:.1f})"
)

You can overlay the raw data by adding a scatter layer before the effect line (this is something that is easy with arrays but would require more work if you only had a chart object).

Multi-term models

When a model has multiple smooth terms, partial_dependence() returns one result per term. You can iterate over them to build a panel of plots:

data_wages = wk.load_dataset("wages")

model_wages = wk.GAM("wage ~ s(age) + s(experience)")
model_wages.fit(data_wages, method="REML")

results_wages = model_wages.partial_dependence()

panels = []
for r in results_wages:
    var_name = list(r.x.keys())[0]
    x_vals = r.x[var_name]

    pd_data = [
        {"x": float(x_vals[i]), "effect": float(r.effect[i]),
         "lower": float(r.lower[i]), "upper": float(r.upper[i])}
        for i in range(len(x_vals))
    ]

    band = alt.Chart({"values": pd_data}).mark_area(
        opacity=0.25, color="steelblue"
    ).encode(x="x:Q", y="lower:Q", y2="upper:Q")

    line = alt.Chart({"values": pd_data}).mark_line(
        color="steelblue", strokeWidth=2
    ).encode(
        x=alt.X("x:Q", title=var_name),
        y=alt.Y("effect:Q", title="Partial effect"),
    )

    ref = alt.Chart({"values": [{"y": 0}]}).mark_rule(
        color="gray", strokeDash=[4, 4]
    ).encode(y="y:Q")

    panels.append(
        (band + ref + line).properties(
            width=300, height=250, title=f"{r.term} (EDF = {r.edf:.1f})"
        )
    )

alt.hconcat(*panels)

Controlling the grid and confidence level

Two keyword-only parameters let you adjust the evaluation:

  • n_points (default 200): the number of evenly spaced points along each covariate’s range. More points give smoother curves whereas fewer points are faster.
  • level (default 0.95): the confidence level for the lower and upper bounds.
# Coarse grid with 90% confidence bands
results_coarse = model.partial_dependence(n_points=50, level=0.90)
r_coarse = results_coarse[0]

print(f"Grid points: {r_coarse.n_grid}")
print(f"Confidence level: {r_coarse.level}")
Grid points: 50
Confidence level: 0.9

Compare two confidence levels side by side:

r_95 = model.partial_dependence(n_points=200, level=0.95)[0]
r_80 = model.partial_dependence(n_points=200, level=0.80)[0]

x_vals = r_95.x["times"]

ci_data = []
for i in range(len(x_vals)):
    ci_data.append({
        "times": float(x_vals[i]), "effect": float(r_95.effect[i]),
        "lower_95": float(r_95.lower[i]), "upper_95": float(r_95.upper[i]),
        "lower_80": float(r_80.lower[i]), "upper_80": float(r_80.upper[i]),
    })

band_95 = alt.Chart({"values": ci_data}).mark_area(
    opacity=0.15, color="steelblue"
).encode(x="times:Q", y="lower_95:Q", y2="upper_95:Q")

band_80 = alt.Chart({"values": ci_data}).mark_area(
    opacity=0.30, color="steelblue"
).encode(x="times:Q", y="lower_80:Q", y2="upper_80:Q")

line = alt.Chart({"values": ci_data}).mark_line(
    color="steelblue", strokeWidth=2
).encode(
    x=alt.X("times:Q", title="times"),
    y=alt.Y("effect:Q", title="Partial effect"),
)

ref = alt.Chart({"values": [{"y": 0}]}).mark_rule(
    color="gray", strokeDash=[4, 4]
).encode(y="y:Q")

(band_95 + band_80 + ref + line).properties(
    width="container", height=300,
    title="Effect of confidence level on band width"
)

2-D smooths

For a 2-D smooth like s(x, y), the x dictionary has two keys corresponding to the two covariates. Each value is a 1-D marginal grid of length approximately sqrt(n_points). The effect, se, lower, and upper arrays are flattened over the full grid (length = n_side^2).

To reconstruct the 2-D surface, use np.meshgrid on the marginal grids and reshape the effect:

data_meuse = wk.load_dataset("meuse")

model_2d = wk.GAM("zinc ~ s(x, y, k=25)")
model_2d.fit(data_meuse, method="REML")

results_2d = model_2d.partial_dependence(n_points=225)
r2 = results_2d[0]

print(f"Term: {r2.term}")
print(f"x keys: {list(r2.x.keys())}")
print(f"Marginal grid lengths: {[len(v) for v in r2.x.values()]}")
print(f"Effect length: {len(r2.effect)}")
Term: s(x, y, k=25)
x keys: ['x', 'y']
Marginal grid lengths: [15, 15]
Effect length: 225

Plot the 2-D partial effect as a heatmap:

x_grid = r2.x["x"]
y_grid = r2.x["y"]

# Build a flat list of grid cells with their effect values
heat_data = []
idx = 0
for yi in range(len(y_grid)):
    for xi in range(len(x_grid)):
        heat_data.append({
            "x": float(x_grid[xi]),
            "y": float(y_grid[yi]),
            "effect": float(r2.effect[idx]),
        })
        idx += 1

alt.Chart({"values": heat_data}).mark_rect().encode(
    x=alt.X("x:O", title="x", axis=alt.Axis(labelAngle=0, values=x_grid[::3].tolist())),
    y=alt.Y("y:O", title="y", sort="descending", axis=alt.Axis(values=y_grid[::3].tolist())),
    color=alt.Color("effect:Q", scale=alt.Scale(scheme="redblue", domainMid=0),
                     title="Partial effect"),
).properties(
    width=400, height=350,
    title=f"2-D partial dependence: {r2.term}"
)

When to use partial_dependence() vs partial_effects()

Scenario Method
Quick visualization during exploration partial_effects()
Publication-quality matplotlib figure partial_dependence()
Export smooth curves to CSV or DataFrame partial_dependence()
Overlay observed data on effect plots partial_dependence()
Numerical analysis (zero crossings, AUC) partial_dependence()
Interactive notebook with Altair tooltips partial_effects()
Compare effects across multiple models partial_dependence()

Both methods compute the same underlying quantities. The difference is purely in what they return: chart objects vs. arrays.

Where to go next