Generalized Additive Models (GAM)¶

GAMs extend GLMs by replacing the linear predictor with smooth functions:

\[\mathbb{E}[Y] = \beta_0 + f_1(x_1) + f_2(x_2) + \dots + \epsilon\]

This chapter covers spline types, smoothing parameter selection, the formula interface, and visualization.

Univariate GAM¶

Use fit_gam() for a single smooth:

from aurora.models.gam import fit_gam

# Generate nonlinear data
np.random.seed(42)
x = np.linspace(0, 10, 300)
y = 2.0 + np.sin(x) + 0.3 * np.random.randn(300)

# Fit GAM
result = fit_gam(x, y, n_basis=15)

print(f"EDF:    {result.edf:.2f}")      # Effective degrees of freedom
print(f"Lambda: {result.lambda_:.4f}")   # Smoothing parameter
print(f"GCV:   {result.gcv_score:.6f}")

# Predict
x_new = np.linspace(0, 10, 100)
y_pred = result.predict(x_new)

Key parameters for fit_gam:

Parameter	Description
`n_basis`	Number of basis functions (default 10). More = more flexible.
`basis_type`	`"bspline"` (default) or `"cubic"` (natural cubic spline)
`degree`	Spline degree (default 3 = cubic)
`penalty_order`	Difference penalty order (default 2, approximates second derivative)
`lambda_`	Smoothing parameter. `None` (default) = automatic via GCV.
`use_sparse`	Use sparse matrices for large problems (n > 500, k > 20)

Spline types¶

B-splines (default) offer local support and numerical stability via the Cox-de Boor recursion. They are the recommended choice for most problems.

Natural cubic splines enforce f’’(x) = 0 at boundaries and minimize the integrated squared second derivative.

# B-spline (default)
result_bs = fit_gam(x, y, n_basis=15, basis_type="bspline")

# Natural cubic spline
result_cs = fit_gam(x, y, n_basis=12, basis_type="cubic")

Smoothing parameter selection¶

Two methods control the smoothness-accuracy tradeoff:

GCV (Generalized Cross-Validation) – default, fast:

\[GCV(\lambda) = \frac{n \cdot RSS}{(n - edf)^2}\]

REML (Restricted Maximum Likelihood) – more stable for small samples.

from aurora.models.gam import fit_additive_gam

# GCV (default, faster)
result_gcv = fit_additive_gam(X, y, smooth_terms=[...], method="GCV")

# REML (more stable)
result_reml = fit_additive_gam(X, y, smooth_terms=[...], method="REML")

Additive GAM with multiple terms¶

Use fit_additive_gam() for models combining smooth and parametric terms:

from aurora.models.gam import fit_additive_gam
from aurora.models.gam.terms import SmoothTerm, ParametricTerm

# Generate data with two nonlinear effects and one linear effect
np.random.seed(42)
n = 300
X = np.random.randn(n, 3)
y = (np.sin(2 * X[:, 0]) + np.cos(X[:, 1])
      + 0.5 * X[:, 2] + np.random.randn(n) * 0.3)

# Define terms
smooth_terms = [
    SmoothTerm(variable=0, n_basis=12),   # s(x0) - nonlinear
    SmoothTerm(variable=1, n_basis=10),   # s(x1) - nonlinear
]
parametric_terms = [
    ParametricTerm(variable=2),            # x2 - linear
]

result = fit_additive_gam(X, y, smooth_terms, parametric_terms)
print(result.summary())

Formula interface¶

The easiest way to fit a GAM is with fit_gam_formula():

from aurora.models.gam import fit_gam_formula

# Using a dict with named variables
data = {
    "y": y,
    "temperature": X[:, 0],
    "pressure": X[:, 1],
    "humidity": X[:, 2],
}
result = fit_gam_formula("y ~ s(temperature) + s(pressure) + humidity", data)
print(result.summary())

# Using a 2D array with column indices
M = np.column_stack([y, X])
result = fit_gam_formula("0 ~ s(1) + s(2) + 3", M)

See Formula Syntax Reference for the complete formula syntax reference.

Tensor product smooths¶

Note

Tensor product smooths (te(x1, x2) in mgcv) are not currently supported by the formula interface or the public fitting functions. Only univariate smooth terms s(x) can be combined additively.

Visualization¶

Plot smooth functions and diagnostic plots:

from aurora.visualization import plot_smooth, plot_all_smooths

# Plot a single smooth term
fig = plot_smooth(result, term=0)

# Plot all smooth terms
figs = plot_all_smooths(result)

# Diagnostic plots (residuals, QQ, scale-location)
from aurora.visualization import plot_diagnostics
fig = plot_diagnostics(result)

Sparse matrices¶

For large problems (n > 500, n_basis > 20), enable sparse operations:

# Large dataset with many basis functions
result = fit_gam(x, y, n_basis=30, use_sparse=True)

Sparse operations use CSR format internally and dispatch to efficient solvers. Only available for "bspline" basis type.