# build_model_matrix()


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


Usage

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


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


``` python
model_matrix.column_names[:5]
```


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