A next-generation Generalized Additive Model (GAM) library for Python.
Developers
Richard Iannone
Maintainer
Next-generation GAMs for Python
Whittaker is a modern Generalized Additive Model (GAM) library for Python. It aims to hold a candle to R’s mgcv in functionality, performance, and usability. It’s built on a modern data stack with Narwhals for backend-agnostic compute, NumPy and SciPy for numerical work, and Altair for publication-quality visualization.
Named after the Whittaker smoother (one of the earliest penalized smoothing methods and a direct mathematical ancestor of P-splines) the library honors that lineage while delivering a modern implementation.
Why Whittaker?
- Works with your dataframe library: Pandas, Polars, PyArrow, or any Narwhals-compatible frame.
- Full smooth catalog: thin plate regression splines, cubic splines, P-splines, tensor products, cyclic splines, random effects, factor smooths, and more.
- Principled smoothness selection: REML by default, with GCV, ML, and fREML as alternatives.
- Beautiful diagnostics:
wk.check(model) and model.plot() produce interactive Altair charts right in your notebook.
- Goes beyond the mean: distributional regression (GAMLSS), quantile regression, conformal prediction, causal inference, streaming GAMs, and functional regression are all built in.
- Principled Bayesian inference: variational inference for fast approximate posteriors, and NUTS MCMC for exact posterior sampling with R-hat, ESS, and divergence diagnostics.
- Type-safe: full annotations, strict Pyright-clean codebase.
The analytical pipeline
- Specify your model with a formula string or programmatic term list.
- Fit with
model.fit(data): automatic smoothness selection via REML.
- Inspect with
model.summary(): EDF, p-values, and smoothing parameters for every term.
- Check with
wk.check(model): basis dimension adequacy, residual plots, convergence info.
- Predict with
model.predict(new_data): point estimates, standard errors, or term contributions.
- Visualize with
model.plot(): partial-effect plots with confidence bands for every smooth.
A first look
Here’s a complete example: generate data with a smooth pattern, fit a GAM, inspect the results, and predict on new observations.
Step 1: Fit a GAM to noisy sinusoidal data
import numpy as np
import whittaker as wk
rng = np.random.default_rng(23)
x = np.linspace(0, 2 * np.pi, 200)
y = np.sin(x) + rng.normal(0, 0.3, 200)
model = wk.GAM("y ~ s(x)")
model.fit({"x": x, "y": y}, method="REML")
model.summary()
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Gaussian(link='identity')
Inference: REML
Observations: 200
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) -0.0160 0.0226 -0.709 0.4789
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 6.99 7 896.560 < 1e-16
Total EDF: 7.99
Scale est: 0.102191
Deviance: 19.6215
Null dev: 111.9322
Dev. expl: 82.5%
GCV score: 0.106445
AIC: 119.39
BIC: 145.74
Step 2: Check basis dimension adequacy
Step 3: Predict with standard errors and visualize the fit
import altair as alt
x_grid = np.linspace(0, 2 * np.pi, 200)
preds = model.predict({"x": x_grid}, se=True)
plot_data = [
{"x": float(x_grid[i]), "fit": float(preds.values[i]),
"lower": float(preds.values[i] - 1.96 * preds.se[i]),
"upper": float(preds.values[i] + 1.96 * preds.se[i])}
for i in range(len(x_grid))
]
obs_data = [{"x": float(x[i]), "y": float(y[i])} for i in range(len(x))]
band = alt.Chart({"values": plot_data}).mark_area(
opacity=0.2, color="darkorange"
).encode(x=alt.X("x:Q", title="x"), y=alt.Y("lower:Q", title="y"), y2="upper:Q")
line = alt.Chart({"values": plot_data}).mark_line(
color="darkorange", strokeWidth=2
).encode(x="x:Q", y="fit:Q")
pts = alt.Chart({"values": obs_data}).mark_circle(
size=15, opacity=0.3, color="steelblue"
).encode(x="x:Q", y="y:Q")
(band + pts + line).properties(width="container", height=300, title="Fitted GAM with 95% confidence band")
Step 4: Try a non-Gaussian family with multiple predictors
n = 500
x1 = rng.uniform(0, 2 * np.pi, n)
x2 = rng.uniform(0, 1, n)
mu = np.exp(0.5 + 0.8 * np.sin(x1) + 2 * x2)
counts = rng.poisson(mu).astype(float)
model_pois = wk.GAM("y ~ s(x1) + s(x2)", family=wk.Poisson())
model_pois.fit({"x1": x1, "x2": x2, "y": counts}, method="REML")
model_pois.summary()
GAM fit summary
============================================================
Formula: y ~ s(x1) + s(x2)
Family: Poisson(link='log')
Inference: REML
Observations: 500
Coefficients: 19
Parametric coefficients:
Term Estimate Std.Err z value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 1.4616 0.0238 61.286 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x1) 6.94 7 713.509 < 1e-16
s(x2) 1.00 2 884.135 < 1e-16
Total EDF: 8.94
Scale est: 1.000000
Deviance: 560.2627
Null dev: 2324.0584
Dev. expl: 75.9%
GCV score: 1.161715
AIC: 2149.21
BIC: 2186.91
Step 5: Cross-validate and save the model
cv = wk.cross_validate(
"y ~ s(x1) + s(x2)",
{"x1": x1, "x2": x2, "y": counts},
family=wk.Poisson(),
n_folds=5,
)
print(f"CV deviance: {cv.cv_score:.4f} (SE: {cv.cv_se:.4f})")
wk.save_gam(model_pois, "/tmp/poisson_model.npz")
loaded = wk.load_gam("/tmp/poisson_model.npz")
print(f"Loaded model EDF: {loaded.edf_total:.1f}")
CV deviance: 1.1691 (SE: 0.0356)
Loaded model EDF: 8.9
What’s inside
Core GAM fitting: fit penalized regression splines with automatic smoothness selection (REML, GCV, ML, fREML). R-style formula syntax with smooth terms, tensor products, linear terms, interactions, offsets, and by-variable smooths. Full summary with EDF and significance tests.
Smooth basis types: thin plate regression splines (TPRS), cubic regression splines, P-splines, cyclic variants, shrinkage smooths, Duchon splines, Gaussian processes, soap film smooths, Markov random fields, adaptive TPRS, random effects, factor smooths, and tensor products.
Response families: Gaussian, Poisson, Binomial, Gamma, Negative Binomial, Beta, Tweedie, Inverse Gaussian, Cox PH, Ordered Categorical, and Multinomial. Each with appropriate link functions and variance structure.
Shape constraints: monotone increasing/decreasing, convex, and concave smooths via constrained P-splines with PAVA projection.
Prediction and inference: point estimates, standard errors, confidence intervals (pointwise and simultaneous), prediction intervals, and term-level contributions. All on response or link scale.
Model diagnostics: model.summary() for EDF and significance tests, wk.check(model) for basis dimension adequacy (k-index test), concurvity analysis, and residual plots.
Bayesian inference: variational inference (method="VI") produces a fast Gaussian approximation to the posterior with richer uncertainty than the Laplace approximation. MCMC (method="MCMC") uses the No-U-Turn Sampler (NUTS) to draw exact posterior samples (no trajectory-length tuning required!). Both methods report R-hat convergence statistics, effective sample sizes, and (for MCMC) divergent transition counts.
Distributional regression (GAMLSS): model location, scale, and shape simultaneously with Gaussian, Gamma, Beta, zero-inflated Poisson, and zero-inflated Negative Binomial families.
Quantile regression: fit conditional quantiles with ELF loss, optional non-crossing constraints, and sigma calibration.
Conformal prediction: distribution-free prediction intervals via split, CV+, and jackknife+ methods.
Causal inference: double/debiased machine learning for ATE and CATE estimation, with mediation analysis.
Streaming and online GAMs: incremental fitting via sufficient statistics with exponential decay for tracking distribution shift.
Multi-response models: joint fitting of multiple responses with optional residual correlation modeling.
Functional regression: scalar-on-function regression with B-spline or Fourier bases for functional covariates.
Large datasets: BigGAM (discretized P-IRLS), PolarsGAM (streaming from Polars/files), and DuckDBGAM (SQL-native streaming) for datasets that exceed memory.
Cross-validation: k-fold CV with deviance, MSE, or MAE scoring.
Serialization: save and load fitted models as compact .npz archives, and convert to/from mgcv-compatible dictionaries for R interoperability.
scikit-learn integration: GAMRegressor and GAMClassifier for use in pipelines and grid search.