# Diagnostic data for custom plots


# Why diagnostic data?

The [check()](../reference/check.md#whittaker.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()](../reference/GAM.md#whittaker.GAM.check_data) method on a fitted [GAM](../reference/GAM.md#whittaker.GAM) returns a [CheckDataResult](../reference/CheckDataResult.md#whittaker.CheckDataResult) dataclass containing all the arrays that [check()](../reference/check.md#whittaker.check) would plot, without producing any charts.


# Basic usage

Fit a model on the motorcycle crash dataset and call [check_data()](../reference/GAM.md#whittaker.GAM.check_data):


``` python
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](../reference/CheckDataResult.md#whittaker.CheckDataResult) dataclass with the following fields:


``` python
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()](../reference/check.md#whittaker.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()](../reference/check.md#whittaker.check):


``` python
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)
```


<style>
  #altair-viz-8d49eb3522944d0585286c9a0542b8ed.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-8d49eb3522944d0585286c9a0542b8ed.vega-embed details,
  #altair-viz-8d49eb3522944d0585286c9a0542b8ed.vega-embed details summary {
    position: relative;
  }
</style>


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:


``` python
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:


``` python
# 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()](../reference/GAM.md#whittaker.GAM.check_data) vs [check()](../reference/check.md#whittaker.check)

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

Both methods compute the same underlying quantities. [check()](../reference/check.md#whittaker.check) returns Altair chart objects; [check_data()](../reference/GAM.md#whittaker.GAM.check_data) returns arrays in a [CheckDataResult](../reference/CheckDataResult.md#whittaker.CheckDataResult) dataclass.


# Where to go next

- [Model diagnostics](diagnostics.md) covers the full diagnostic workflow, including [gam_check()](../reference/GAM.md#whittaker.GAM.gam_check) for basis dimension adequacy checks.
- [Advanced diagnostics](advanced-diagnostics.md) dives deeper into residual types and influence measures.
- [Partial dependence as data](partial-dependence.md) follows the same pattern: raw arrays for the data behind [partial_effects()](../reference/partial_effects.md#whittaker.partial_effects) plots.