# GAM


Generalized Additive Model with automatic smoothness selection.


Usage

``` python
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](Family.md#whittaker.Family) object.

Use [GAM](GAM.md#whittaker.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](GAM.md#whittaker.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](Formula.md#whittaker.Formula), [SmoothTerm](SmoothTerm.md#whittaker.SmoothTerm), [LinearTerm](LinearTerm.md#whittaker.LinearTerm), [InteractionTerm](InteractionTerm.md#whittaker.InteractionTerm), and [OffsetTerm](OffsetTerm.md#whittaker.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](GAM.md#whittaker.GAM) supports prediction with standard errors and intervals (`predict()`), partial-effect plotting ([plot()](GAM.md#whittaker.GAM.plot)), residual and basis-dimension diagnostics ([check()](check.md#whittaker.check), [gam_check()](GAM.md#whittaker.GAM.gam_check), [k_check()](GAM.md#whittaker.GAM.k_check)), hypothesis tests for parametric and smooth terms ([parametric_tests()](GAM.md#whittaker.GAM.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](Formula.md#whittaker.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()](Gaussian.md#whittaker.Gaussian) (identity link) if not given. Other options include [Binomial](Binomial.md#whittaker.Binomial), [Poisson](Poisson.md#whittaker.Poisson), [Gamma](Gamma.md#whittaker.Gamma), and [Tweedie](Tweedie.md#whittaker.Tweedie)-family classes, each defining the variance function, deviance, and link used during P-IRLS.


## Examples


``` python
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](#aic) | Akaike Information Criterion of the fitted model. |
| [bic](#bic) | Bayesian Information Criterion of the fitted model. |
| [coefficients](#coefficients) | Estimated coefficients \boldsymbol{\beta}. |
| [deviance](#deviance) | Model deviance at convergence. |
| [deviance_explained](#deviance_explained) | Proportion of null deviance explained by the model (analogous to R²). |
| [edf](#edf) | Effective degrees of freedom (EDF) for each smooth term. |
| [edf_total](#edf_total) | Total effective degrees of freedom across all model terms. |
| [family](#family) | The response distribution family used to fit this model. |
| [fitted_values](#fitted_values) | Fitted values \mu on the response scale. |
| [formula](#formula) | The parsed model formula. |
| [gcv_score](#gcv_score) | Generalized Cross-Validation score at the fitted smoothing parameters. |
| [is_fitted](#is_fitted) | Whether the model has been fitted. |
| [mcmc_result](#mcmc_result) | The `MCMCResult` when the model was fitted with `method="MCMC"`, else `None`. |
| [null_deviance](#null_deviance) | Deviance of the intercept-only (null) model. |
| [residuals](#residuals) | Response residuals (y - \mu) on the training data. |
| [scale](#scale) | Estimated scale (dispersion) parameter \phi. |
| [smoothing_params](#smoothing_params) | Selected or fixed smoothing parameters \lambda_j, one per penalty. |
| [vi_result](#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](GAM.md#whittaker.GAM.null_deviance) via [deviance_explained](GAM.md#whittaker.GAM.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()](GAM.md#whittaker.GAM.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()](Gaussian.md#whittaker.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](Formula.md#whittaker.Formula) object produced by parsing the formula string passed to `GAM.__init__()` (or the [Formula](Formula.md#whittaker.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](CausalGAM.md#whittaker.CausalGAM.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](GAM.md#whittaker.GAM.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()](GAM.md#whittaker.GAM.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()](#anova) | Compare this model with one or more other fitted GAMs via deviance-difference tests. |
| [check_data()](#check_data) | Return diagnostic data for custom plotting. |
| [concurvity()](#concurvity) | Compute concurvity diagnostics for all smooth terms. |
| [derivatives()](#derivatives) | Estimate derivatives of smooth terms with respect to a variable. |
| [dispersion_test()](#dispersion_test) | Test for overdispersion in Poisson or Binomial models. |
| [fit()](#fit) | Fit the GAM to data via Penalized Iteratively Reweighted Least Squares (P-IRLS). |
| [gam_check()](#gam_check) | Run all-in-one GAM diagnostics. |
| [get_residuals()](#get_residuals) | Compute residuals of the specified type. |
| [goodness_of_fit()](#goodness_of_fit) | Return all goodness-of-fit statistics in a single object. |
| [influence()](#influence) | Compute hat values and Cook's distance for each observation. |
| [k_check()](#k_check) | Check basis dimension adequacy for each smooth term. |
| [loo()](#loo) | Compute PSIS-LOO cross-validation for a Bayesian fit. |
| [marginal_effects()](#marginal_effects) | Compute marginal (partial) effects of a variable. |
| [pairwise_comparisons()](#pairwise_comparisons) | Compute pairwise contrasts between conditions. |
| [parametric_tests()](#parametric_tests) | Compute Wald tests for parametric (non-smooth) coefficients. |
| [partial_dependence()](#partial_dependence) | Compute partial dependence data for each smooth term. |
| [plot()](#plot) | Plot the estimated partial effect of each smooth term, with a confidence band. |
| [posterior_predict()](#posterior_predict) | Draw from the posterior predictive distribution at new data points. |
| [posterior_samples()](#posterior_samples) | Draw coefficient vectors from the posterior. |
| [ppc()](#ppc) | Run a posterior predictive check against the training data. |
| [predict()](#predict) | Predict from the fitted model on new data. |
| [quantile_residuals()](#quantile_residuals) | Compute randomized quantile residuals (Dunn & Smyth 1996). |
| [simulate()](#simulate) | Draw from the posterior distribution of the fitted model. |
| [simultaneous_ci()](#simultaneous_ci) | Compute simultaneous confidence bands for smooth terms. |
| [smooth_tests()](#smooth_tests) | Compute approximate p-values for all smooth terms. |
| [smoothing_sensitivity()](#smoothing_sensitivity) | Sweep smoothing parameters and record how predictions and fit statistics change. |
| [summary()](#summary) | Return a text summary of the fitted model, analogous to `summary.gam()` in R's mgcv. |
| [vif()](#vif) | Compute variance inflation factors for parametric (linear) terms. |
| [waic()](#waic) | Compute the Widely Applicable Information Criterion (WAIC). |

------------------------------------------------------------------------


### anova()


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


Usage

``` python
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](GAM.md#whittaker.GAM) objects to compare against this model.


#### Returns


`AnovaResult`  
Sequential deviance-comparison table.


------------------------------------------------------------------------


### check_data()


Return diagnostic data for custom plotting.


Usage

``` python
check_data()
```


Provides the same data that [check()](check.md#whittaker.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

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

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

``` python
dispersion_test()
```


#### Returns


`DispersionTestResult`  
Object with `dispersion`, `chi2_stat`, and [p_value](PPCResult.md#whittaker.PPCResult.p_value).


------------------------------------------------------------------------


### fit()


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


Usage

``` python
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](CausalGAM.md#whittaker.CausalGAM.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

``` python
gam_check(
    *,
    n_sim=100,
)
```


Returns a [GamCheckResult](GamCheckResult.md#whittaker.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

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

``` python
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](GoodnessOfFit.md#whittaker.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

``` python
influence()
```


#### Returns


`InfluenceResult`  
Object with `hat_values` and `cooks_distance` arrays.


------------------------------------------------------------------------


### k_check()


Check basis dimension adequacy for each smooth term.


Usage

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

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


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

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

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

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

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


Returns the same underlying data that [partial_effects()](partial_effects.md#whittaker.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

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


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


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

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


------------------------------------------------------------------------


### posterior_predict()


Draw from the posterior predictive distribution at new data points.


Usage

``` python
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()](PosteriorPredictResult.md#whittaker.PosteriorPredictResult.mean), [std()](PosteriorPredictResult.md#whittaker.PosteriorPredictResult.std), [quantile()](PosteriorPredictResult.md#whittaker.PosteriorPredictResult.quantile), and [interval()](PosteriorPredictResult.md#whittaker.PosteriorPredictResult.interval).


#### Examples

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


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


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


``` python
pp.quantile(0.5)
```


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


------------------------------------------------------------------------


### posterior_samples()


Draw coefficient vectors from the posterior.


Usage

``` python
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](Tweedie.md#whittaker.Tweedie.p) is the number of model coefficients.


------------------------------------------------------------------------


### ppc()


Run a posterior predictive check against the training data.


Usage

``` python
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](PPCResult.md#whittaker.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

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

``` python
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](TermsPredictionResult.md#whittaker.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](PredictionResult.md#whittaker.PredictionResult) with [values](TermsPredictionResult.md#whittaker.TermsPredictionResult.values), `se`, `linear_predictor`, `lower`, and `upper`. For `type="terms"`, a [TermsPredictionResult](TermsPredictionResult.md#whittaker.TermsPredictionResult) with per-smooth contributions and standard errors.


#### Examples


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

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

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

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

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

``` python
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](SensitivityResult.md#whittaker.SensitivityResult.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


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

``` python
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()](GAM.md#whittaker.GAM.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

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

``` python
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](GAM.md#whittaker.GAM.waic), and `pointwise`.


#### Raises


`ValueError`  
If the model was not fitted with a Bayesian method.


#### Examples

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