GAM

Generalized Additive Model with automatic smoothness selection.

Usage

Source

GAM(
    formula,
    family=None,
)

A GAM extends the generalized linear model by replacing some or all linear predictor terms with smooth, data-driven functions of the covariates:

g(\mathbb{E}[y_i]) = \eta_i = \beta_0 + \sum_j \beta_j x_{ij} + \sum_k f_k(z_{ik})

where g is a link function, \beta_j x_{ij} are ordinary parametric (linear) terms, and each f_k is an unspecified smooth function represented by a spline basis (s()) or a tensor product of bases for multivariate smooths (te(), ti(), t2()). This lets the model capture nonlinear relationships without having to guess a parametric form ahead of time, while still supporting the full range of exponential-family response distributions (Gaussian, Binomial, Poisson, Gamma, Tweedie, and more) via a Family object.

Use GAM when you suspect a covariate’s effect on the response is nonlinear, when you want interaction surfaces between two or more continuous covariates, or when you want automatic, data-driven control of model complexity rather than manually choosing a polynomial degree or a fixed set of basis functions.

A GAM is specified with a formula string in an R/mgcv-like syntax, e.g. "y ~ s(x1) + s(x2, bs='cr', k=15) + te(x3, x4) + group", where s() denotes a univariate (or by=-varying) smooth, te()/ti()/t2() denote tensor-product smooths of two or more variables, and bare names denote ordinary parametric terms. See Formula, SmoothTerm, LinearTerm, InteractionTerm, and OffsetTerm for the term types this formula parses into.

Fitting (fit()) proceeds by Penalized Iteratively Reweighted Least Squares (P-IRLS): each smooth’s wiggliness is controlled by a quadratic penalty \lambda_k \boldsymbol{\beta}_k^T \mathbf{S}_k \boldsymbol{\beta}_k on its coefficients, and the smoothing parameters \lambda_k are themselves estimated from the data. By default this is via Generalized Cross-Validation (GCV), or via Restricted Maximum Likelihood (REML) or Marginal Likelihood (ML) when smooths are treated as correlated random effects. Larger \lambda_k shrinks a smooth toward a simpler (e.g. linear or constant) shape; smaller \lambda_k allows more flexibility. This automatic selection is what distinguishes a GAM from simply choosing a fixed spline basis: the effective complexity of each term (its effective degrees of freedom, or EDF) is learned rather than fixed in advance.

Once fitted, a GAM supports prediction with standard errors and intervals (predict()), partial-effect plotting (plot()), residual and basis-dimension diagnostics (check(), gam_check(), k_check()), hypothesis tests for parametric and smooth terms (parametric_tests(), smooth_tests()), and a text summary (summary()) analogous to summary.gam() in R’s mgcv.

Parameters

formula: str or Formula

Model formula, either as a string (e.g. "y ~ s(x1) + s(x2) + x3") or an already-parsed Formula object. The left-hand side names the response column; the right-hand side lists smooth terms (s(), te(), ti(), t2()), parametric terms (bare column names), interactions (x1 * x2), and optionally an offset(...) term. Use 0 + or - 1 on the right-hand side to suppress the intercept.

family: Family or None = None

Response distribution and link function. Defaults to Gaussian() (identity link) if not given. Other options include Binomial, Poisson, Gamma, and Tweedie-family classes, each defining the variance function, deviance, and link used during P-IRLS.

Examples

import numpy as np
from whittaker import GAM

rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 10, 200))
y = np.sin(x) + rng.normal(scale=0.2, size=200)

gam = GAM("y ~ s(x)").fit({"x": x, "y": y})
print(gam.summary())

GAM fit summary
============================================================
Formula:    y ~ s(x)
Family:     Gaussian(link='identity')
Inference:  GCV
Observations: 200
Coefficients: 10

Parametric coefficients:
  Term                       Estimate    Std.Err    t value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                  0.1325     0.0146      9.073  1.441e-16

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  s(x)                       8.35      9   1996.690    < 1e-16

Total EDF:  9.35
Scale est:  0.042647
Deviance:   8.1308
Null dev:   93.9889
Dev. expl:  91.3%
GCV score:  0.044739
AIC:        -54.03
BIC:        -23.20

Attributes

Name	Description
aic	Akaike Information Criterion of the fitted model.
bic	Bayesian Information Criterion of the fitted model.
coefficients	Estimated coefficients \boldsymbol{\beta}.
deviance	Model deviance at convergence.
deviance_explained	Proportion of null deviance explained by the model (analogous to R²).
edf	Effective degrees of freedom (EDF) for each smooth term.
edf_total	Total effective degrees of freedom across all model terms.
family	The response distribution family used to fit this model.
fitted_values	Fitted values \mu on the response scale.
formula	The parsed model formula.
gcv_score	Generalized Cross-Validation score at the fitted smoothing parameters.
is_fitted	Whether the model has been fitted.
mcmc_result	The MCMCResult when the model was fitted with method="MCMC", else None.
null_deviance	Deviance of the intercept-only (null) model.
residuals	Response residuals (y - \mu) on the training data.
scale	Estimated scale (dispersion) parameter \phi.
smoothing_params	Selected or fixed smoothing parameters \lambda_j, one per penalty.
vi_result	The VIResult when the model was fitted with method="VI", else None.

---

aic

Akaike Information Criterion of the fitted model.

aic: float

Balances goodness of fit against model complexity (using the total effective degrees of freedom in place of the raw parameter count). Lower values indicate a preferable trade-off. Use it to compare non-nested models fitted to the same data and family.

For Bayesian fits the log-likelihood is evaluated at the posterior mean.

---

bic

Bayesian Information Criterion of the fitted model.

bic: float

Like aic, but penalizes model complexity more heavily as sample size grows (using log(n) in place of 2 as the per-degree-of-freedom penalty), which tends to favor simpler models than AIC for larger datasets.

For Bayesian fits the log-likelihood is evaluated at the posterior mean.

---

coefficients

Estimated coefficients \boldsymbol{\beta}.

coefficients: NDArray

A single flat vector holding the intercept, parametric term coefficients, and every smooth term’s basis coefficients concatenated in formula order. Use self._model_matrix.smooths (or predict(type="terms")) to map sub-ranges of this vector back to individual terms.

---

deviance

Model deviance at convergence.

deviance: float

Twice the difference between the saturated log-likelihood and the fitted model’s log-likelihood, evaluated at the final coefficients. Lower values indicate a better fit to the training data; compare against null_deviance via deviance_explained.

For Bayesian fits (VI, MCMC) the deviance is evaluated at the posterior mean.

---

deviance_explained

Proportion of null deviance explained by the model (analogous to R²).

deviance_explained: float

Computed as 1 - deviance / null_deviance. Ranges from 0 (no improvement over an intercept-only model) up to 1 (a perfect fit), and provides a family-agnostic measure of goodness of fit that generalizes R² beyond the Gaussian case.

For Bayesian fits the deviance is evaluated at the posterior mean.

---

edf

Effective degrees of freedom (EDF) for each smooth term.

edf: list[float]

Each value is the trace of the portion of the hat matrix attributable to that term, reflecting how much shrinkage its smoothing parameter applied: values near the term’s basis dimension indicate little penalization, values near 1 indicate near-linear shrinkage.

---

edf_total

Total effective degrees of freedom across all model terms.

edf_total: float

The sum of the per-term EDF values (plus the intercept and parametric terms), i.e., the trace of the full hat (influence) matrix. Used in summary(), gam_check(), and in computing residual degrees of freedom for interval and test calculations.

---

family

The response distribution family used to fit this model.

family: Family

Determines the variance function, deviance, and link function used during P-IRLS. Defaults to Gaussian() when no family argument was given to GAM.__init__().

---

fitted_values

Fitted values \mu on the response scale.

fitted_values: NDArray

Equal to g^{-1}(\eta) where \eta = \mathbf{X}\boldsymbol{\beta} is the linear predictor evaluated on the training data used in fit().

---

formula

The parsed model formula.

formula: Formula

This is the Formula object produced by parsing the formula string passed to GAM.__init__() (or the Formula object passed directly). It lists the response column and the parsed smooth (s(), te(), ti(), t2()), parametric, interaction, and offset terms that make up the right-hand side.

---

gcv_score

Generalized Cross-Validation score at the fitted smoothing parameters.

gcv_score: float

Computed as \text{GCV} = n \cdot D / (n - \text{tr}(\mathbf{H}))^2, where D is the deviance and \mathbf{H} is the hat matrix. This is the criterion minimized when fit(method="GCV") selects smoothing parameters, and is reported even when a different method was used.

---

is_fitted

Whether the model has been fitted.

is_fitted: bool

True once fit() has completed successfully; False beforehand. Most other properties and methods (coefficients, predict(), summary(), etc.) require this to be True and raise a RuntimeError via _check_fitted() otherwise.

---

mcmc_result

The MCMCResult when the model was fitted with method="MCMC", else None.

mcmc_result: MCMCResult | None

---

null_deviance

Deviance of the intercept-only (null) model.

null_deviance: float

Fit on the same data and with the same family and weights, but with every covariate effect (smooth and parametric) removed. Serves as the baseline against which deviance_explained measures the reduction in deviance achieved by the fitted model.

---

residuals

Response residuals (y - \mu) on the training data.

residuals: NDArray

These are the raw (unstandardized) residuals. For Pearson, deviance, or working residuals (or residuals on new data) use get_residuals() instead.

---

scale

Estimated scale (dispersion) parameter \phi.

scale: float

For families with a known scale (Binomial, Poisson) this is fixed at 1.0. For families with unknown scale (Gaussian, Gamma, Tweedie) it is estimated from the Pearson residuals and is used to scale coefficient standard errors and prediction intervals.

---

smoothing_params

Selected or fixed smoothing parameters \lambda_j, one per penalty.

smoothing_params: list[float]

If fit() was called with smoothing_params=None (the default), these are the values chosen automatically via GCV, REML, or ML; otherwise they are the fixed values that were passed in. A te()/t2() term contributes more than one entry (one per marginal penalty), so this list is generally longer than the number of smooth terms.

---

vi_result

The VIResult when the model was fitted with method="VI", else None.

vi_result: VIResult | None

Methods

Name	Description
anova()	Compare this model with one or more other fitted GAMs via deviance-difference tests.
check_data()	Return diagnostic data for custom plotting.
concurvity()	Compute concurvity diagnostics for all smooth terms.
derivatives()	Estimate derivatives of smooth terms with respect to a variable.
dispersion_test()	Test for overdispersion in Poisson or Binomial models.
fit()	Fit the GAM to data via Penalized Iteratively Reweighted Least Squares (P-IRLS).
gam_check()	Run all-in-one GAM diagnostics.
get_residuals()	Compute residuals of the specified type.
goodness_of_fit()	Return all goodness-of-fit statistics in a single object.
influence()	Compute hat values and Cook’s distance for each observation.
k_check()	Check basis dimension adequacy for each smooth term.
loo()	Compute PSIS-LOO cross-validation for a Bayesian fit.
marginal_effects()	Compute marginal (partial) effects of a variable.
pairwise_comparisons()	Compute pairwise contrasts between conditions.
parametric_tests()	Compute Wald tests for parametric (non-smooth) coefficients.
partial_dependence()	Compute partial dependence data for each smooth term.
plot()	Plot the estimated partial effect of each smooth term, with a confidence band.
posterior_predict()	Draw from the posterior predictive distribution at new data points.
posterior_samples()	Draw coefficient vectors from the posterior.
ppc()	Run a posterior predictive check against the training data.
predict()	Predict from the fitted model on new data.
quantile_residuals()	Compute randomized quantile residuals (Dunn & Smyth 1996).
simulate()	Draw from the posterior distribution of the fitted model.
simultaneous_ci()	Compute simultaneous confidence bands for smooth terms.
smooth_tests()	Compute approximate p-values for all smooth terms.
smoothing_sensitivity()	Sweep smoothing parameters and record how predictions and fit statistics change.
summary()	Return a text summary of the fitted model, analogous to summary.gam() in R’s mgcv.
vif()	Compute variance inflation factors for parametric (linear) terms.
waic()	Compute the Widely Applicable Information Criterion (WAIC).

---

anova()

Compare this model with one or more other fitted GAMs via deviance-difference tests.

Usage

Source

anova(*others)

All models must use the same family and be fitted to the same data. Models are automatically sorted by complexity (edf). For known-scale families (Poisson, Binomial) a chi-squared test is used. For unknown-scale families (Gaussian, Gamma) an F-test is used.

Parameters

*others: GAM

One or more fitted GAM objects to compare against this model.

Returns

AnovaResult

Sequential deviance-comparison table.

---

check_data()

Return diagnostic data for custom plotting.

Usage

Source

check_data()

Provides the same data that check() renders as Altair charts (deviance and pearson residuals, fitted values, response, and QQ coordinates) as structured arrays.

Returns

CheckDataResult

---

concurvity()

Compute concurvity diagnostics for all smooth terms.

Usage

Source

concurvity(
    *,
    full=True,
)

Concurvity is the GAM analogue of collinearity. High values (> 0.8) indicate that a smooth’s effect may be confounded with other model terms, making its estimate unstable.

Parameters

full: bool = True

If True (default), measure each smooth against all other model terms combined. If False, compute pairwise concurvity between each pair of smooths.

Returns

ConcurvityResult

Object with worst, observed, and estimate arrays, plus labels.

---

derivatives()

Estimate derivatives of smooth terms with respect to a variable.

Usage

Source

derivatives(
    variable,
    *,
    order=1,
    n_points=200,
    level=0.95,
    eps=None,
    unconditional=False,
)

Uses central finite differences on the basis matrix with delta-method standard errors.

Parameters

variable: str

The covariate to differentiate with respect to.

order: int = 1

Derivative order: 1 for first derivative (rate of change), 2 for second derivative (curvature).

n_points: int = 200

Number of evaluation points along the variable’s range.

level: float = 0.95

Confidence level for the bands.

eps: float | None = None

Finite difference step size. If None, chosen automatically.

unconditional: bool = False

If True, use unconditional covariance (Marra & Wood 2012).

Returns

list[DerivativeResult]

One result per smooth term involving the variable.

---

dispersion_test()

Test for overdispersion in Poisson or Binomial models.

Usage

Source

dispersion_test()

Returns

DispersionTestResult

Object with dispersion, chi2_stat, and p_value.

---

fit()

Fit the GAM to data via Penalized Iteratively Reweighted Least Squares (P-IRLS).

Usage

Source

fit(
    data,
    *,
    smoothing_params=None,
    method="GCV",
    weights=None,
    select=False,
    vi_options=None,
    mcmc_options=None,
)

Builds the model matrix from the formula, then alternates between (a) an IRLS step that linearizes the exponential-family log-likelihood around the current fit and (b) a penalized weighted-least-squares solve that shrinks each smooth toward simplicity according to its smoothing parameter. When smoothing_params is not fixed, this inner P-IRLS loop is itself wrapped in an outer loop that re-estimates the smoothing parameters (by GCV, REML, or ML) at each iteration, until both the coefficients and the smoothing parameters converge.

Parameters

data: dict[str, numpy.ndarray]

Column-oriented data as {name: 1-D array}. All columns referenced by the formula (response, smooth covariates, parametric terms, and any by= factor) must be present.

smoothing_params: list[float] or None = None

Fixed smoothing parameters \lambda_k, one per penalty (a te()/t2() term contributes more than one). If None (default), smoothing parameters are selected automatically via method.

method: str = "GCV"

Criterion used to select smoothing parameters when smoothing_params is not fixed. One of:

• "GCV" (default): minimizes the Generalized Cross-Validation score \text{GCV} = n \cdot D / (n - \text{tr}(\mathbf{H}))^2, where D is the deviance and \mathbf{H} is the influence (hat) matrix. Fast and does not require treating smooths as random effects, but can occasionally undersmooth.
• "REML": maximizes the Restricted Maximum Likelihood, treating each smooth’s penalized coefficients as correlated Gaussian random effects and integrating out the fixed (unpenalized) effects. Generally the most reliable choice and the one recommended when using select=True.
• "ML": maximizes the Marginal Likelihood, similar to REML but without correcting for uncertainty in the fixed effects; tends to undersmooth slightly relative to REML.

weights: numpy.ndarray or None = None

Observation (prior) weights, shape (n,). Must be positive. When provided, the model minimizes the weighted deviance \sum_i w_i d_i and uses weighted IRLS throughout.

select: bool = False

If True, augment each smooth’s wiggliness penalty with a second penalty on its null space (the component, such as a pure linear trend, that the ordinary penalty never shrinks). With both penalties free, GCV/REML/ML can drive a term’s smoothing parameters high enough to remove it from the model entirely, giving automatic term selection analogous to the lasso (Marra & Wood, 2011). Recommended together with method="REML". Has no additional effect on bases whose null space is already zero (e.g. bs="re", bs="fs", or the shrinkage bases "ts"/"cs").

vi_options: dict or None = None

Extra keyword arguments forwarded to ~whittaker.fitting.vi.vi_fit when method="VI". Accepted keys: n_quad, lr, max_iter, tol, patience, seed, cov_structure, phi_inference. Ignored for all other methods.

mcmc_options: dict or None = None

Extra keyword arguments forwarded to ~whittaker.fitting.mcmc.mcmc_fit when method="MCMC". Accepted keys: n_samples, n_warmup, n_chains, leapfrog_steps, target_accept, seed. Ignored for all other methods.

Returns

GAM

Returns self for method chaining, e.g. GAM(formula).fit(data).predict(new_data).

---

gam_check()

Run all-in-one GAM diagnostics.

Usage

Source

gam_check(
    *,
    n_sim=100,
)

Returns a GamCheckResult containing deviance residuals, fitted values, response values, basis dimension checks, and summary statistics. Print the result for a quick diagnostic summary.

Parameters

n_sim: int = 100

Number of permutations for the k-check p-values (the default is 100).

Returns

GamCheckResult

Diagnostic results with a readable __repr__.

---

get_residuals()

Compute residuals of the specified type.

Usage

Source

get_residuals(type="deviance")

Parameters

type: str = "deviance"

One of "response", "pearson", "deviance", or "working".

Returns

NDArray

Residual vector of shape (n,).

---

goodness_of_fit()

Return all goodness-of-fit statistics in a single object.

Usage

Source

goodness_of_fit()

Collects deviance, null deviance, deviance explained, adjusted R-squared, AIC, BIC, GCV (when available), scale, EDF, and the number of observations into a GoodnessOfFit dataclass. This is a convenience method that avoids calling each property individually.

The adjusted R-squared is computed as 1 - (1 - deviance_explained) * (n - 1) / (n - edf_total - 1), generalizing the classical formula by using the effective degrees of freedom in place of the raw parameter count.

Returns

GoodnessOfFit

---

influence()

Compute hat values and Cook’s distance for each observation.

Usage

Source

influence()

Returns

InfluenceResult

Object with hat_values and cooks_distance arrays.

---

k_check()

Check basis dimension adequacy for each smooth term.

Usage

Source

k_check(
    *,
    n_sim=400,
)

For each smooth, computes a k-index based on the ratio of a neighbor-differencing variance estimate of the residuals (ordered by covariate) to the overall residual variance. A k-index well below 1 suggests that the basis dimension k may be too small. A simulation-based p-value is computed: low p-values indicate potential under-smoothing.

Parameters

n_sim: int = 400

Number of random permutations for the p-value simulation (the default is 400).

Returns

list[KCheckResult]

One result per smooth term.

---

loo()

Compute PSIS-LOO cross-validation for a Bayesian fit.

Usage

Source

loo(
    n_draws=2000,
    *,
    seed=None,
)

Uses Pareto-Smoothed Importance Sampling to approximate leave-one-out cross-validation from the existing posterior draws without refitting. Requires a model fitted with method="VI" or method="MCMC".

Parameters

n_draws: int = 2000

Number of posterior draws to use. For MCMC fits, all stored draws are used when n_draws= exceeds the available count. For VI fits, samples are drawn from the Gaussian posterior approximation. The default is 2000.

seed: int or None = None

Random seed for reproducibility when sampling from the posterior (VI fits only). MCMC fits use all stored draws directly without sampling.

Returns

LOOResult

Object with elpd_loo, se_elpd_loo, p_loo, pointwise, pareto_k, and n_bad_k. Inspect pareto_k to identify observations where the PSIS approximation may be unreliable (values > 0.7).

Raises

ValueError

If the model was not fitted with method="VI" or method="MCMC".

Examples

import numpy as np
from whittaker import GAM
from whittaker.families import Poisson

rng = np.random.default_rng(0)
x = np.linspace(0, 5, 100)
y = rng.poisson(np.exp(0.4 * x))

m1 = GAM("y ~ s(x)", family=Poisson()).fit({"x": x, "y": y}, method="VI")
m2 = GAM("y ~ x", family=Poisson()).fit({"x": x, "y": y}, method="VI")

loo1 = m1.loo()
loo2 = m2.loo()
print(loo1)

from whittaker import loo_compare
print(loo_compare(loo1, loo2))

LOOResult
  ELPD_LOO:    -195.11  (SE 7.68)
  p_LOO:       5.61
  Bad k > 0.7: 0 / 100 observations
LOOComparison
  ELPD diff: -3.34  (SE 0.54)
  model 2 preferred

---

marginal_effects()

Compute marginal (partial) effects of a variable.

Usage

Source

marginal_effects(
    variable,
    *,
    at=None,
    n_points=200,
    level=0.95,
    unconditional=False,
)

Evaluates the smooth term(s) involving variable over a grid while holding other variables at their means or at values specified via at.

Parameters

variable: str

The focal covariate.

at: dict | None = None

Dict mapping other variable names to fixed values (or lists of values for a grid). Variables not listed are held at their mean.

n_points: int = 200

Number of evaluation points along the variable’s range.

level: float = 0.95

Confidence level for the bands.

unconditional: bool = False

If True, use unconditional covariance.

Returns

list[MarginalEffectResult]

One result per smooth term per at combination.

---

pairwise_comparisons()

Compute pairwise contrasts between conditions.

Usage

Source

pairwise_comparisons(
    variable,
    pairs,
    *,
    n_points=200,
    level=0.95,
    unconditional=False,
)

Each pair is (condition1, condition2) where each condition is a dict of covariate values. The contrast f(x|cond1) - f(x|cond2) is evaluated over a grid of the focal variable.

Parameters

variable: str

The focal covariate (the x-axis for the contrast).

pairs: list[tuple[dict, dict]]

List of (cond1, cond2) dicts specifying the two conditions.

n_points: int = 200

Number of evaluation points.

level: float = 0.95

Confidence level for the bands.

unconditional: bool = False

If True, use unconditional covariance.

Returns

list[ContrastResult]

One result per smooth term per pair.

---

parametric_tests()

Compute Wald tests for parametric (non-smooth) coefficients.

Usage

Source

parametric_tests()

Uses the t-distribution for families with unknown scale (Gaussian, Gamma) and the z-distribution for known-scale families (Binomial, Poisson).

Returns

list[ParametricTestResult]

One result per parametric coefficient (intercept + linear terms).

---

partial_dependence()

Compute partial dependence data for each smooth term.

Usage

Source

partial_dependence(
    *,
    n_points=200,
    level=0.95,
)

Returns the same underlying data that partial_effects() plots, but as structured arrays rather than Altair chart objects. This is useful for custom plotting with matplotlib or other libraries, or for downstream numerical analysis of the smooth effects.

Parameters

n_points: int = 200

Number of evenly spaced evaluation points per smooth. For 2-D smooths, each marginal gets approximately sqrt(n_points) points.

level: float = 0.95

Confidence level for the bands (default 0.95).

Returns

list[PartialDependenceResult]

One result per smooth term, in formula order.

---

plot()

Plot the estimated partial effect of each smooth term, with a confidence band.

Usage

Source

plot(
    *,
    n_points=200,
    level=0.95,
)

For every smooth in the formula, evaluates the term’s contribution to the linear predictor over an evenly spaced grid spanning its covariate’s observed range (holding other terms out, i.e., this shows the additive component f_k(x) itself, not the full fitted response), together with a pointwise confidence band derived from the model’s coefficient covariance. This is the standard way to visually inspect the shape of each estimated smooth (e.g., whether it is roughly linear, monotonic, or has a distinct peak) without needing to call predict(type="terms") and plot manually.

Parameters

n_points: int = 200

Number of evenly spaced evaluation points per smooth (the default is 200).

level: float = 0.95

Confidence level for the bands, e.g. 0.95 for a 95% band (the default is 0.95).

Returns

altair.VConcatChart or altair.Chart

A vertically concatenated chart with one panel per smooth term (or a single Chart if the model has exactly one smooth).

Examples

import numpy as np
from whittaker import GAM

rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 10, 200))
y = np.sin(x) + rng.normal(scale=0.2, size=200)

gam = GAM("y ~ s(x)").fit({"x": x, "y": y})
gam.plot()

---

posterior_predict()

Draw from the posterior predictive distribution at new data points.

Usage

Source

posterior_predict(
    new_data=None,
    *,
    n_draws=1000,
    seed=None,
)

Generates the full (n, n_draws) posterior predictive sample by:

1. drawing coefficient vectors from the posterior (or Laplace approximation).
2. computing the linear predictor η = X @ β* (+ offset if present).
3. transforming to the response scale: μ = g⁻¹(η).
4. sampling observation noise from the family distribution at each μ.

The result includes both coefficient uncertainty and observation-level noise, giving the distribution of new observations rather than of the conditional mean.

Parameters

new_data: InputData | None = None

Column-oriented data for prediction. If None, uses the training data.

n_draws: int = 1000

Number of posterior draws (the default is 1000).

seed: int | None = None

Random seed for reproducibility.

Returns

PosteriorPredictResult

Object with the (n, n_draws) sample matrix and convenience methods for mean(), std(), quantile(), and interval().

Examples

Draw 2000 posterior predictive samples at five new points and compute the posterior predictive mean:

import numpy as np
from whittaker import GAM

rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 10, 200))
y = np.sin(x) + rng.normal(scale=0.2, size=200)

model = GAM("y ~ s(x)").fit({"x": x, "y": y}, method="VI")
x_new = np.linspace(0, 10, 5)
pp = model.posterior_predict({"x": x_new}, n_draws=2000, seed=42)
pp.mean()

array([-0.0643466 ,  0.52989434, -1.06311014,  0.91570952, -0.56499046])

The 95% equal-tailed posterior predictive interval:

lower, upper = pp.interval()
lower, upper

(array([-0.52793903,  0.09913929, -1.46438834,  0.50322525, -0.98361705]),
 array([ 0.3494651 ,  0.97353113, -0.64824351,  1.33621325, -0.13261815]))

The posterior predictive median:

pp.quantile(0.5)

array([-0.06112926,  0.53410627, -1.06429845,  0.91237061, -0.56671984])

---

posterior_samples()

Draw coefficient vectors from the posterior.

Usage

Source

posterior_samples(
    n=1000,
    *,
    seed=None,
)

For VI fits, samples from the variational posterior N(m, C). For Laplace fits (REML/GCV/ML), samples from N(β̂, V_β).

Parameters

n: int = 1000

Number of draws.

seed: int | None = None

Random seed.

Returns

NDArray

Shape (p, n), where p is the number of model coefficients.

---

ppc()

Run a posterior predictive check against the training data.

Usage

Source

ppc(
    n_sim=1000,
    *,
    seed=None,
)

Draws n_sim= datasets from the posterior predictive distribution and compares five test statistics (mean, sd, min, max, and proportion of zeros) between the observed data and the replicated datasets. The result is a PPCResult whose Bayesian p-values indicate whether the model generates data consistent with the observations.

Works with any fitting method. For Bayesian fits (method="VI" or method="MCMC") draws come from the full posterior. For frequentist fits the Laplace approximation to the posterior is used.

Parameters

n_sim: int = 1000

Number of replicated datasets to draw. The default is 1000.

seed: int or None = None

Random seed for reproducibility.

Returns

PPCResult

Contains y_rep (shape (n, n_sim)), observed (shape (n,)), and per-statistic Bayesian p-values.

Examples

result = model.ppc(n_sim=1000, seed=0)
print(result)
# PPCResult (n=300, n_sim=1000)
#
# Statistic        Observed  Mean(rep)    p-value
# -----------------------------------------------
# mean                 2.34       2.31      0.520
# sd                   1.45       1.43      0.480
# ...

---

predict()

Predict from the fitted model on new data.

Usage

Source

predict(
    new_data,
    *,
    se=False,
    type="response",
    interval=None,
    level=0.95,
    unconditional=False,
    n_sim=1000,
    seed=None,
)

Applies the estimated coefficients to a model matrix built from new_data, using the same basis transformations (knots, factor levels, centering constraints) fitted on the training data. Supports predictions on the response or linear-predictor scale, per-term decompositions, and pointwise or simultaneous uncertainty intervals.

Parameters

new_data: dict[str, numpy.ndarray]

Column-oriented new data. Must contain all covariate columns referenced by the formula (the response column is not needed).

se: bool = False

If True, compute standard errors on the linear predictor scale (for type="response" and type="link") or per-term standard errors (for type="terms"). Default False.

type: str = "response"

Prediction type. One of:

• "response" (default): predictions on the response scale, \mu = g^{-1}(\eta).
• "link": predictions on the linear predictor scale, \eta = \mathbf{X} \boldsymbol{\beta}.
• "terms": individual smooth term contributions to the linear predictor, returned separately rather than summed (see TermsPredictionResult).

interval: str or None = None

Interval type. None (default) returns no intervals. Ignored (and must be None) when type="terms". Otherwise one of:

• "confidence": interval for the mean response, reflecting uncertainty in \eta only.
• "prediction": interval for a new individual observation, adding the response-distribution variance on top of the uncertainty in \eta.
• "simultaneous": a band with level coverage for the entire curve simultaneously (via posterior simulation), rather than pointwise coverage.
• "credible": posterior credible interval for Bayesian fits (method="VI" or method="MCMC"). Draws n_sim coefficient vectors from the posterior, computes the predicted mean at each draw, and returns the (1-level)/2 and (1+level)/2 quantiles as lower and upper. More accurate than "confidence" for non-Gaussian families because it does not rely on normal approximation after the link-inverse transform. Raises ValueError for non-Bayesian fits.

All interval types are computed on the linear predictor scale and transformed to the response scale for type="response".

level: float = 0.95

Nominal coverage probability for the interval, e.g. 0.95 for a 95% interval (default).

unconditional: bool = False

If True, include smoothing-parameter uncertainty in standard errors and intervals (Marra & Wood, 2012), using the unconditional covariance matrix V_c in place of the conditional V_p. This produces wider, more honest intervals that account for the fact that \lambda was itself estimated from the data. Requires that the model was fitted with method="REML" or method="ML".

n_sim: int = 1000

Number of posterior draws used when interval="credible". Default 1000. Ignored for all other interval types.

seed: int or None = None

Random seed passed to the posterior sampler when interval="credible". Default None. Ignored for all other interval types.

Returns

PredictionResult or TermsPredictionResult

For type="response" or type="link", a PredictionResult with values, se, linear_predictor, lower, and upper. For type="terms", a TermsPredictionResult with per-smooth contributions and standard errors.

Examples

import numpy as np
from whittaker import GAM

rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 10, 200))
y = np.sin(x) + rng.normal(scale=0.2, size=200)

gam = GAM("y ~ s(x)").fit({"x": x, "y": y})
new_x = np.linspace(0, 10, 5)
result = gam.predict({"x": new_x}, se=True, interval="confidence")
result.values, result.lower, result.upper

(array([-0.03617302,  0.53232071, -1.04777748,  0.91697446, -0.57118698]),
 array([-0.18767157,  0.44050127, -1.13552273,  0.83216494, -0.71733199]),
 array([ 0.11532552,  0.62414015, -0.96003223,  1.00178398, -0.42504197]))

---

quantile_residuals()

Compute randomized quantile residuals (Dunn & Smyth 1996).

Usage

Source

quantile_residuals(
    *,
    seed=None,
)

For a correctly specified model, these should be approximately standard normal.

Parameters

seed: int | None = None

Random seed for the jittering step (discrete families).

Returns

NDArray

Quantile residuals, shape (n,).

---

simulate()

Draw from the posterior distribution of the fitted model.

Usage

Source

simulate(
    new_data=None,
    *,
    n_sim=1000,
    seed=None,
    unconditional=False,
)

Generates posterior simulations of the response by:

1. Drawing coefficient vectors β* ~ MVN(β̂, V_β) from the Bayesian posterior.
2. Computing η* = X @ β* (+ offset if present) for each draw.
3. Transforming to the response scale: μ* = g⁻¹(η*).

When unconditional=True, response noise is added by sampling from the family distribution at each μ*.

Parameters

new_data: InputData | None = None

Column-oriented data for prediction. If None, uses the training data.

n_sim: int = 1000

Number of posterior draws (the default is 1000).

seed: int | None = None

Random seed for reproducibility.

unconditional: bool = False

If True, add response-distribution noise on top of posterior uncertainty in the mean. This produces simulations of new observations rather than of the conditional mean.

Returns

NDArray

Simulated values on the response scale, shape (n, n_sim) where n is the number of observations in the prediction data.

---

simultaneous_ci()

Compute simultaneous confidence bands for smooth terms.

Usage

Source

simultaneous_ci(
    new_data,
    *,
    term=None,
    level=0.95,
    n_sim=10000,
    unconditional=False,
    seed=0,
)

Unlike pointwise intervals, these bands have (approximate) level= coverage probability for the entire function simultaneously, not just at individual points.

Parameters

new_data: InputData

Prediction data.

term: int | str | None = None

Which smooth term to compute bands for. An integer index (0-based) or the term label string. If None and the model has exactly one smooth, that term is used.

level: float = 0.95

Nominal simultaneous coverage probability (the default is 0.95).

n_sim: int = 10000

Number of posterior simulations for the critical value (the default is 10_000).

unconditional: bool = False

If True, include smoothing parameter uncertainty (the default is False).

seed: int = 0

Random seed for reproducibility (the default is 0).

Returns

SimultaneousCIResult

Contains estimate, se, lower, upper, term_label, crit_value.

---

smooth_tests()

Compute approximate p-values for all smooth terms.

Usage

Source

smooth_tests()

Uses the Wood (2013) approach: eigendecomposition of the Bayesian covariance block for each smooth, with a chi-squared reference distribution.

Returns

list[SmoothTestResult]

One result per smooth term.

---

smoothing_sensitivity()

Sweep smoothing parameters and record how predictions and fit statistics change.

Usage

Source

smoothing_sensitivity(
    new_data=None,
    *,
    multipliers=None,
    n_steps=11,
    log_range=(-2.0, 2.0),
)

Re-fits the model at each multiplier value, scaling all smoothing parameters uniformly. This reveals how sensitive the predictions are to the specific smoothing-parameter values chosen by the fitting criterion (GCV, REML, or ML).

Parameters

new_data: dict[str, numpy.ndarray] or None = None

Data at which to evaluate predictions. If None, uses the training data.

multipliers: sequence of float or None = None

Explicit multiplier values. If None, n_steps values are generated log-uniformly over log_range.

n_steps: int = 11

Number of log-spaced multiplier values when multipliers is None.

log_range: tuple[float, float] = (-2.0, 2.0)

(lo, hi) on the log10 scale for auto-generated multipliers. The default (-2, 2) sweeps from 0.01x to 100x.

Returns

SensitivityResult

Examples

import numpy as np
import whittaker as wk

rng = np.random.default_rng(0)
x = np.linspace(0, 2 * np.pi, 200)
y = np.sin(x) + rng.normal(0, 0.3, 200)
model = wk.GAM("y ~ s(x)").fit({"x": x, "y": y})

sens = model.smoothing_sensitivity()
print(sens)

SensitivityResult(11 steps, 200 observations)
  Multiplier range:    [0.01, 100] (baseline=1)
  EDF range:           [3.0, 9.9]
  Max |prediction change|: 0.7337

---

summary()

Return a text summary of the fitted model, analogous to summary.gam() in R’s mgcv.

Usage

Source

summary()

The summary reports, in order: the formula and family; a table of parametric (non-smooth) coefficients with estimates, standard errors, test statistics (t for unknown-scale families such as Gaussian and Gamma, z for known-scale families such as Binomial and Poisson), and p-values; a table of approximate significance for each smooth term (effective degrees of freedom, reference degrees of freedom, a chi-squared-type statistic, and a p-value); and overall fit statistics (total EDF, deviance, null deviance, proportion of deviance explained, GCV score, estimated scale, AIC, and BIC). Use this for a quick, human-readable check of which terms are significant and how well the model fits, without extracting individual result objects via parametric_tests() and smooth_tests().

Returns

ModelSummary

Multi-line text summary. Displays cleanly as the last expression in a Jupyter or Quarto cell without needing print().

---

vif()

Compute variance inflation factors for parametric (linear) terms.

Usage

Source

vif()

Returns

list[VIFResult]

One result per parametric term. Empty if fewer than 2 parametric terms.

---

waic()

Compute the Widely Applicable Information Criterion (WAIC).

Usage

Source

waic(
    n_draws=1000,
    *,
    seed=None,
)

WAIC (estimates out-of-sample predictive accuracy from the posterior log-likelihood matrix. It is asymptotically equivalent to LOO-CV but cheaper to compute because it does not require importance-sampling corrections. The result includes the ELPD on the WAIC scale, a standard error, and the effective number of parameters p_WAIC.

Parameters

n_draws: int = 1000

Number of posterior draws used to estimate the log-likelihood matrix. Ignored for MCMC fits (all stored samples are used).

seed: int or None = None

Random seed for drawing posterior samples (VI only).

Returns

WAICResult

Contains elpd_waic, se_elpd_waic, p_waic, waic, and pointwise.

Raises

ValueError

If the model was not fitted with a Bayesian method.

Examples

model = GAM("y ~ s(x)").fit(data, method="VI")
w = model.waic()
print(w)