Diagnostic data for custom plots

Why diagnostic data?

The check() function from whittaker.plotting produces a set of four Altair diagnostic charts (QQ plot, residuals vs. fitted, histogram of deviance residuals, and response vs. fitted) with a single call. This is the fastest way to visually assess a model. But sometimes you need the raw arrays behind those plots:

  • Matplotlib or seaborn workflows: your project already uses matplotlib and you want a consistent style across all figures.
  • Automated checks: you want to run programmatic tests on the residuals (e.g., normality tests, heteroscedasticity detection) rather than eyeballing charts.
  • Export: you need to write the diagnostic values to a CSV or pass them to another tool.
  • Fine-grained control: you want to customize axis limits, colors, annotations, or panel layout beyond what the built-in plots offer.

The check_data() method on a fitted GAM returns a CheckDataResult dataclass containing all the arrays that check() would plot, without producing any charts.

Basic usage

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

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")

cd = model.check_data()
print(cd)
CheckDataResult(n_obs=133)

The return value is a CheckDataResult dataclass with the following fields:

print(f"n_obs:              {cd.n_obs}")
print(f"deviance_residuals: shape {cd.deviance_residuals.shape}")
print(f"pearson_residuals:  shape {cd.pearson_residuals.shape}")
print(f"fitted_values:      shape {cd.fitted_values.shape}")
print(f"response:           shape {cd.response.shape}")
print(f"qq_theoretical:     shape {cd.qq_theoretical.shape}")
print(f"qq_observed:        shape {cd.qq_observed.shape}")
n_obs:              133
deviance_residuals: shape (133,)
pearson_residuals:  shape (133,)
fitted_values:      shape (133,)
response:           shape (133,)
qq_theoretical:     shape (133,)
qq_observed:        shape (133,)

Every field is a NumPy array of length n_obs (except qq_theoretical and qq_observed, which are sorted for direct plotting). The arrays correspond exactly to the data rendered by the four panels in check().

Recreating the four diagnostic plots

With the arrays in hand, you can reproduce all four diagnostic plots. Here they are arranged in a 2x2 grid using Altair, mirroring the layout of check():

import altair as alt

# 1. QQ plot of deviance residuals
qq_data = [
    {"theoretical": float(cd.qq_theoretical[i]), "observed": float(cd.qq_observed[i])}
    for i in range(len(cd.qq_theoretical))
]
qq_min = min(cd.qq_theoretical.min(), cd.qq_observed.min())
qq_max = max(cd.qq_theoretical.max(), cd.qq_observed.max())
qq_ref = [{"x": float(qq_min), "y": float(qq_min)}, {"x": float(qq_max), "y": float(qq_max)}]

qq_plot = (
    alt.Chart({"values": qq_data}).mark_circle(size=15, opacity=0.6, color="steelblue").encode(
        x=alt.X("theoretical:Q", title="Theoretical quantiles"),
        y=alt.Y("observed:Q", title="Observed quantiles"),
    )
    + alt.Chart({"values": qq_ref}).mark_line(color="firebrick", strokeDash=[4, 4]).encode(
        x="x:Q", y="y:Q"
    )
).properties(width=300, height=250, title="QQ plot")

# 2. Residuals vs. fitted values
resid_data = [
    {"fitted": float(cd.fitted_values[i]), "residual": float(cd.deviance_residuals[i])}
    for i in range(cd.n_obs)
]

resid_plot = (
    alt.Chart({"values": resid_data}).mark_circle(size=15, opacity=0.6, color="steelblue").encode(
        x=alt.X("fitted:Q", title="Fitted values"),
        y=alt.Y("residual:Q", title="Deviance residuals"),
    )
    + alt.Chart({"values": [{"y": 0}]}).mark_rule(color="firebrick", strokeDash=[4, 4]).encode(
        y="y:Q"
    )
).properties(width=300, height=250, title="Residuals vs. fitted")

# 3. Histogram of deviance residuals
hist_data = [{"residual": float(r)} for r in cd.deviance_residuals]

hist_plot = alt.Chart({"values": hist_data}).mark_bar(
    opacity=0.8, color="steelblue"
).encode(
    x=alt.X("residual:Q", bin=alt.Bin(maxbins=30), title="Deviance residuals"),
    y=alt.Y("count():Q", title="Frequency"),
).properties(width=300, height=250, title="Histogram of residuals")

# 4. Response vs. fitted values
resp_data = [
    {"fitted": float(cd.fitted_values[i]), "response": float(cd.response[i])}
    for i in range(cd.n_obs)
]
mn, mx = float(cd.fitted_values.min()), float(cd.fitted_values.max())
resp_ref = [{"x": mn, "y": mn}, {"x": mx, "y": mx}]

resp_plot = (
    alt.Chart({"values": resp_data}).mark_circle(size=15, opacity=0.6, color="steelblue").encode(
        x=alt.X("fitted:Q", title="Fitted values"),
        y=alt.Y("response:Q", title="Response"),
    )
    + alt.Chart({"values": resp_ref}).mark_line(color="firebrick", strokeDash=[4, 4]).encode(
        x="x:Q", y="y:Q"
    )
).properties(width=300, height=250, title="Response vs. fitted")

# 2x2 grid
(qq_plot | resid_plot) & (hist_plot | resp_plot)

Because you have full control over the chart specifications, you can adjust colors, add annotations, change bin widths, or rearrange the panel layout to suit your needs.

Automated residual checks

Having the residuals as arrays makes it straightforward to run programmatic diagnostics. For example, you can test whether the deviance residuals are approximately normally distributed using the Shapiro-Wilk test, and look for heteroscedasticity by checking whether the variance of residuals changes across the range of fitted values:

from scipy import stats

# Shapiro-Wilk test for normality of deviance residuals
stat, p_value = stats.shapiro(cd.deviance_residuals)
print(f"Shapiro-Wilk statistic: {stat:.4f}")
print(f"p-value:                {p_value:.4f}")

if p_value < 0.05:
    print("Evidence against normality of residuals (p < 0.05)")
else:
    print("No strong evidence against normality (p >= 0.05)")
Shapiro-Wilk statistic: 0.9698
p-value:                0.0047
Evidence against normality of residuals (p < 0.05)

To check for heteroscedasticity, split the residuals into groups by fitted value and compare their variances:

# Split residuals into lower and upper halves by fitted value
median_fitted = np.median(cd.fitted_values)
lower = cd.deviance_residuals[cd.fitted_values <= median_fitted]
upper = cd.deviance_residuals[cd.fitted_values > median_fitted]

# Levene's test for equal variances
stat_lev, p_lev = stats.levene(lower, upper)
print(f"Levene's test statistic: {stat_lev:.4f}")
print(f"p-value:                 {p_lev:.4f}")

if p_lev < 0.05:
    print("Evidence of heteroscedasticity (p < 0.05)")
else:
    print("No strong evidence of heteroscedasticity (p >= 0.05)")
Levene's test statistic: 11.3780
p-value:                 0.0010
Evidence of heteroscedasticity (p < 0.05)

You can combine these tests into a reusable function that flags potential problems automatically, rather than relying on visual inspection every time you fit a model.

check_data() vs check()

Scenario Method
Quick visual diagnostic in a notebook check()
Publication-quality matplotlib figure check_data()
Automated test suite or CI pipeline check_data()
Export residuals to CSV or DataFrame check_data()
Statistical tests on residuals check_data()
Interactive exploration with Altair tooltips check()
Custom panel layout or overlays check_data()

Both methods compute the same underlying quantities. check() returns Altair chart objects; check_data() returns arrays in a CheckDataResult dataclass.

Where to go next