Generalized Additive Models

Penalized regression spline GAMs via P-IRLS matching R’s mgcv::gam().

Generalized Additive Models (GAM).

Provides smooth term specification and basis construction for penalized regression spline GAMs following Wood (2017).

Usage:

from pystatistics.gam import gam, s, te, ti

sol = gam(
    y,
    smooths=[s('x1', k=15), te('x2', 'x3')],
    smooth_data={'x1': x1, 'x2': x2, 'x3': x3},
)
print(sol.summary())
class pystatistics.gam.SmoothInfo(term_name, var_name, basis_type, k, edf, ref_df, chi_sq, p_value, coef_indices, lambdas, s_scales)[source]

Bases: object

Information about a single smooth term in a fitted GAM.

One SmoothInfo is produced per s() term after fitting. It records the basis metadata and the approximate significance test for the smooth.

Parameters:
term_name

Display name, e.g. 's(x1)'.

Type:

str

var_name

Bare predictor name, e.g. 'x1'.

Type:

str

basis_type

Basis identifier: 'cr' (cubic regression spline) or 'tp' (thin plate regression spline).

Type:

str

k

The basis dimension as declared in s(x, k=...) (the smooth contributes k - 1 coefficients after the sum-to-zero identifiability constraint, same as mgcv).

Type:

int

edf

Effective degrees of freedom for this term.

Type:

float

ref_df

Reference degrees of freedom, tr(2H - HH) over the term’s block (mgcv’s Ref.df convention).

Type:

float

chi_sq

Approximate significance statistic. For fixed-dispersion families this is the chi-squared-referenced quadratic form beta' pinv(V_beta) beta; for estimated-dispersion families it is that form divided by Ref.df, i.e. an F STATISTIC (the summary labels the column accordingly, matching mgcv). A simplified form of mgcv’s Wood (2013) test — agrees with mgcv’s reported statistic to ~0.2% on validated cases, not bit-identical.

Type:

float

p_value

Approximate p-value for the term.

Type:

float

coef_indices

(start, end) slice into the full coefficient vector identifying this term’s (constrained) coefficients.

Type:

tuple[int, int]

lambdas

Selected smoothing parameter(s) for this term — a length-1 tuple for an ordinary s() smooth, one entry per margin for a tensor-product te()/ti() smooth (each margin has its own penalty and smoothing parameter, as in mgcv).

Type:

tuple[float, …]

s_scales

Penalty normalisation factor(s), aligned with lambdas (lambda / s_scale is in function-space units).

Type:

tuple[float, …]

term_name: str
var_name: str
basis_type: str
k: int
edf: float
ref_df: float
chi_sq: float
p_value: float
coef_indices: tuple[int, int]
lambdas: tuple[float, ...]
s_scales: tuple[float, ...]
pystatistics.gam.s(*var_names, k=10, bs=None, by=None, by_type=None)[source]

Convenience constructor for smooth terms, matching mgcv::s().

Usage:

from pystatistics.gam import s, gam
result = gam(y, smooths=[s('x1', k=15), s('x2', 'x3')], ...)
Parameters:
  • *var_names (str) – One predictor name for a univariate smooth, or several for an isotropic multivariate thin-plate smooth s(x, z, ...).

  • k (int) – Basis dimension (default 10).

  • bs (str | None) – Basis type. Univariate: 'cr' (default), 'tp', 'cc', or 'ps'. Multivariate: 'tp' (the default and only choice). For a factor by (by_type='factor') the default is 'tp', matching mgcv::s().

  • by (str | None) – Optional by variable (univariate smooths only).

  • by_type (str | None) –

    How to interpret by:

    • 'continuous' – the varying-coefficient term by * f(x) (mgcv’s s(x, by=z) for a numeric z), keeping the full basis (no centering) since the by-multiplication removes the constant confound.

    • 'factor' – a separate smooth per level of an integer-coded by (mgcv’s s(x, by=factor(g))); the grouping variable’s per-level means are added automatically, so no separate main effect is required.

    • None (default) – treated as continuous, EXCEPT that a factor-looking by (an integer-valued, low-cardinality contiguous coding) fails loud, asking you to pick 'factor' or 'continuous' explicitly rather than silently fitting a meaningless varying coefficient.

Returns:

A SmoothTerm (one variable) or IsotropicSmooth (several variables).

class pystatistics.gam.SmoothTerm(var_name, k=10, bs='cr', by=None, by_type=None)[source]

Bases: object

Specification for a smooth term s(x) in a GAM formula.

A pure, immutable-by-convention specification: variable name, basis type and basis dimension. All fit-derived quantities (basis matrices, penalties, constraint reparameterisations) live inside the fit — a SmoothTerm can safely be reused across multiple gam() calls and datasets without stale-state hazards.

Parameters:
var_name

Name of the predictor variable.

k

Basis dimension (>= 3), exactly as mgcv’s s(x, k=...). The fitted smooth contributes k - 1 coefficients after its sum-to-zero identifiability constraint (same as mgcv).

bs

Basis type: 'cr' (cubic regression spline, default), 'tp' (thin plate regression spline), 'cc' (cyclic cubic regression spline, for periodic/seasonal covariates), or 'ps' (P-spline).

__init__(var_name, k=10, bs='cr', by=None, by_type=None)[source]

Create a smooth term specification.

Parameters:
  • var_name (str) – Name of the predictor variable.

  • k (int) – Basis dimension (default 10, matching mgcv). Must be >= 3.

  • bs (str) – Basis type – 'cr' for cubic regression spline (default) or 'tp' for thin plate regression spline.

  • by (str | None) – Optional by variable name (see s()).

  • by_type (str | None) – How to interpret by'continuous' (varying coefficient by * f(x)), 'factor' (a separate smooth per level of an integer-coded by), or None to let gam decide; a factor-looking by with by_type=None fails loud rather than being silently treated as continuous.

Raises:

ValidationError – If var_name is empty, k is out of range, bs is not a recognised basis type, or by_type is invalid / given without by.

Return type:

None

var_name
k
bs
by
by_type
pystatistics.gam.te(*var_names, k=5, bs='cr')[source]

Tensor-product smooth te(x, z, ...) (mgcv default marginal k=5).

Parameters:
  • *var_names (str) – One or more marginal predictor names (a single name yields a centred 1-D smooth, exactly as mgcv’s te(x) does).

  • k – Marginal basis dimension — a scalar (all margins) or one per margin. Default 5, matching mgcv’s te marginal default.

  • bs (str) – Marginal basis type — a scalar or one per margin.

Returns:

A TensorSmooth (>= 2 margins) or SmoothTerm (one margin).

pystatistics.gam.ti(*var_names, k=5, bs='cr')[source]

Tensor-product interaction ti(x, z, ...) (main effects removed).

Same arguments as te(); used for functional-ANOVA decompositions te(x) + te(z) + ti(x, z) (the te(x)/te(z) main-effect terms are often written ti(x)/ti(z), i.e. single-margin smooths).

Parameters:
class pystatistics.gam.TensorSmooth(var_names, ks, bss, interaction)[source]

Bases: object

Specification for a te() / ti() tensor-product smooth.

Parameters:
  • var_names (Sequence[str])

  • ks (Sequence[int])

  • bss (Sequence[str])

  • interaction (bool)

var_names

Marginal predictor names (>= 2).

ks

Marginal basis dimension per margin (mgcv marginal k).

bss

Marginal basis type per margin ('cr'/'tp'/'cc'/ 'ps').

interaction

True for ti (main effects removed), False for te.

var_names
ks
bss
interaction
property label: str
pystatistics.gam.gam(y, X=None, *, smooths=None, smooth_data=None, family='gaussian', method='GCV', sp=None, tol=1e-08, max_iter=200, names=None)[source]

Fit a Generalized Additive Model.

Matches R’s mgcv::gam() for common cases. A GAM models the response as:

g(E[y]) = beta_0 + f_1(x_1) + f_2(x_2) + ... + X @ beta

where each f_j is a smooth function estimated from the data using penalized regression splines, with per-smooth sum-to-zero identifiability constraints (a smooth declared with k basis functions contributes k - 1 coefficients, exactly as mgcv).

Parameters:
  • y (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Response variable, 1-D array of n observations.

  • X (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – Parametric design matrix (n, p). May include an intercept column. If None, an intercept-only parametric part is used.

  • smooths (list[SmoothTerm] | None) – Smooth term specifications from s().

  • smooth_data (dict[str, _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]] | None) – Mapping from smooth variable names to (n,) arrays.

  • family (str | Family) – Exponential family, e.g. 'gaussian', 'binomial', 'poisson', 'gamma'. method='REML' supports the fixed-dispersion families and the Gaussian-identity model.

  • method (str) – Smoothing-parameter selection criterion: 'GCV' (mgcv GCV.Cp semantics — UBRE is used automatically for known-scale families, as mgcv does) or 'REML' (Laplace, Wood 2011).

  • sp (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – Optional fixed smoothing parameters, one per smooth (mgcv’s sp=). Skips selection; use for reproducing a specific fit.

  • tol (float) – P-IRLS convergence tolerance (relative penalized deviance).

  • max_iter (int) – Maximum P-IRLS iterations.

  • names (list[str] | None) – Optional display names for the parametric coefficients.

Returns:

A GAMSolution.

Raises:
  • ValidationError – invalid inputs, unknown family/method/basis, k exceeding the unique values of a smooth variable, or a sp vector of the wrong length/sign.

  • ConvergenceError – a genuinely divergent P-IRLS (loud, never a silent answer).

Return type:

GAMSolution

class pystatistics.gam.GAMSolution(_result, _smooth_infos, _names=None)[source]

Bases: SolutionReprMixin

User-facing GAM result with convenient properties and R-style output.

Wraps a Result[GAMParams] and provides:

  • Coefficient access.

  • Smooth-term information (EDF, approximate significance tests).

  • R-style summary() matching mgcv::summary.gam().

Parameters:
__init__(_result, _smooth_infos, _names=None)[source]

Initialise the solution wrapper.

Parameters:
  • _result (Result[GAMParams]) – The generic Result envelope carrying GAMParams.

  • _smooth_infos (list[SmoothInfo]) – One SmoothInfo per smooth term.

  • _names (list[str] | None) – Optional names for parametric coefficients.

Return type:

None

property params: GAMParams

The underlying frozen parameter payload.

property coefficients: ndarray[tuple[Any, ...], dtype[_ScalarT]]

Full coefficient vector (parametric + basis).

property fitted_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]

Response-scale predictions.

property residuals: ndarray[tuple[Any, ...], dtype[_ScalarT]]

Working residuals (y - mu_hat).

property edf: ndarray[tuple[Any, ...], dtype[_ScalarT]]

Effective degrees of freedom per smooth term.

property total_edf: float

Total effective degrees of freedom.

property smooth_terms: list[SmoothInfo]

Information about each smooth term.

property deviance: float

Model deviance.

property null_deviance: float

Null-model deviance.

property deviance_explained: float

1 - dev / null_dev.

Type:

Proportion of null deviance explained

property aic: float

Akaike information criterion.

property gcv: float

GCV score.

property scale: float

Estimated dispersion parameter.

property theta: float | None

Negative-binomial dispersion theta (mgcv::getTheta).

The estimated theta for a family='nb' fit, the supplied value for a fixed NegativeBinomial(theta=...) family, or None for every other family. Larger theta means less overdispersion (theta -> inf recovers Poisson).

property ubre: float

UBRE / scaled-AIC score (the criterion mgcv’s GCV.Cp minimises for known-scale families).

property lambdas: ndarray[tuple[Any, ...], dtype[_ScalarT]]

Selected (or user-fixed) smoothing parameters, one per smooth.

Reported in the same penalty scaling mgcv reports sp in; divide by params.s_scales for function-space units.

property covariance: ndarray[tuple[Any, ...], dtype[_ScalarT]]

Bayesian posterior covariance of the coefficients (mgcv Vp).

property standard_errors: ndarray[tuple[Any, ...], dtype[_ScalarT]]

Standard errors for ALL coefficients (posterior, mgcv-style).

Covers the full coefficient vector (parametric block + smooth-basis coefficients), taken from the diagonal of the Bayesian posterior covariance. For the parametric-only slice see coef / conf_int.

property coef: dict[str, float]

Parametric coefficients as a name-to-value dict.

Covers only the parametric block (intercept plus any linear terms) that drives summary()’s parametric-coefficient table; smooth-basis coefficients are excluded. Names come from the model’s parametric names, falling back to B[i] when the block is unnamed.

property conf_int: ndarray[tuple[Any, ...], dtype[_ScalarT]]

95% Wald confidence intervals for the PARAMETRIC coefficients.

Shape (n_param, 2): beta ± q * SE(beta) using the same quantile summary() uses — the normal quantile for known-scale families and the t quantile (residual df) for estimated-scale families. Smooth-basis coefficients are excluded; smooth terms are reported via summary()’s approximate-significance table.

property backend_name: str

Backend identifier (GAM is CPU-only).

property reml_score: float | None

Laplace REML criterion (None unless method='REML').

property converged: bool

Whether the P-IRLS algorithm converged.

property outer_converged: bool

Whether the smoothing-parameter search converged off-boundary.

property n_iter: int

Number of P-IRLS iterations executed.

property r_squared_adj: float

Adjusted R-squared based on deviance explained.

Matches the formula used by mgcv:

R_adj = 1 - (deviance / df_resid) / (null_deviance / df_null)
summary()[source]

R-style summary matching summary(gam(...)).

Includes:

  • Family and link function.

  • Parametric coefficient table with standard errors and tests.

  • Approximate significance tests for each smooth term.

  • Adjusted R-squared, deviance explained, GCV, and scale.

Returns:

Multi-line summary string.

Return type:

str