build_model_matrix()

Assemble the full design matrix and penalty structure from a formula.

Usage

build_model_matrix(
    formula,
    data,
    *,
    apply_constraints=True,
    select=False,
)

A generalized additive model is fit by turning the covariates into a purely numeric model matrix (also called a design matrix) X — an (n, p) array of n observations by p coefficients — such that the linear predictor is simply a matrix-vector product:

\eta = X \beta + \text{offset}

Every term in the formula contributes one or more columns to X. An intercept contributes a single column of ones; a linear term x3 contributes the column x3 itself; an interaction x1:x2 contributes the elementwise product x1 * x2. A smooth term such as s(x), by contrast, expands into several columns: the basis functions of the underlying spline (or other smooth) basis, evaluated at each observation’s value of x. Fitting a GAM is then no different from fitting a (penalized) linear model in these expanded columns — the non-linearity of f(x) lives entirely in how the basis functions are constructed, not in how the coefficients enter the model.

build_model_matrix is the bridge between a parsed ~whittaker.formula.terms.Formula and the numeric matrices the P-IRLS fitting engine (~whittaker.pirls.pirls_fit) actually needs: the design matrix X plus, for each smooth term, a quadratic penalty matrix S_j that penalizes wiggliness via beta.T @ S_j @ beta. Each S_j starts out sized to just that term’s own basis functions but is expanded (“embedded”) to the full (p, p) model dimension, with zeros everywhere outside that term’s column block, so that summing lambda_j * S_j over all terms gives a single penalty matrix that can be added directly to the unpenalized normal equations.

Parameters

formula: Formula: A parsed ~whittaker.formula.terms.Formula, typically produced by ~whittaker.formula.parser.parse.
data: dict[str, numpy.ndarray]: Column-oriented data as {name: 1-D array}. Every column referenced by the formula must be present. All arrays must have the same length.
apply_constraints: bool = True: If True (the default), apply sum-to-zero identifiability constraints to each smooth term so the intercept is identifiable.
select: bool = False: If True, add an extra penalty on each smooth’s null space so that terms can be penalized to zero entirely (double penalty approach, Marra & Wood 2011). This enables automatic smooth selection via GCV or REML. Smooths that already have null_space_dim == 0 (e.g. bs="ts", bs="cs", bs="re", bs="fs") are unaffected.

Returns

ModelMatrix: Bundled design matrix, penalties, and metadata.

Raises

KeyError: If a required column is missing from data.
ValueError: If an unsupported basis type is requested.

Notes

An ordinary smoothing penalty S for a term like s(x) typically has a non-trivial null space — directions in coefficient space that the penalty does not shrink at all (e.g. the linear component of a cubic spline). This means increasing that term’s smoothing parameter lambda_j can flatten the smooth toward a straight line, but never all the way to zero, so ordinary GCV/REML smoothing-parameter selection cannot remove an irrelevant term from the model entirely.

When select=True, build_model_matrix adds a second penalty matrix per term (see _null_space_penalty) built from the eigendecomposition of the term’s own penalty S: the eigenvectors with (numerically) zero eigenvalue span exactly the null space of S, and projecting onto that eigenspace gives a penalty S_null that penalizes only those previously-unpenalized directions. With both S and S_null present (each with its own smoothing parameter), driving both lambda_j and the null-space smoothing parameter to large values shrinks the entire term, including its linear component, to zero — allowing REML- or GCV-based fitting to perform automatic term selection much like a lasso penalty does for linear models. See Marra, G. and Wood, S.N. (2011), “Practical variable selection for generalized additive models”, Computational Statistics & Data Analysis, 55(7), 2372-2387.

Examples

import numpy as np
from whittaker.formula.parser import parse
from whittaker.model_matrix import build_model_matrix

rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 1, 50))
y = np.sin(2 * np.pi * x) + rng.normal(scale=0.2, size=50)

formula = parse("y ~ s(x)")
model_matrix = build_model_matrix(formula, {"x": x, "y": y})

model_matrix.X.shape

(50, 10)

model_matrix.column_names[:5]

['(Intercept)', 's(x)[0]', 's(x)[1]', 's(x)[2]', 's(x)[3]']