Markov random field basis for areal spatial data.
Usage
A Markov random field (MRF) smooth represents spatial structure over discrete areal units — counties, districts, postcodes, grid cells on a lattice — where the covariate is a categorical label rather than a continuous coordinate, and the only spatial information available is which units are adjacent to which. It is equivalent to mgcv’s bs="mrf" basis. Each unique region gets its own basis function (an indicator column), and smoothness across the map is enforced directly through the neighborhood graph rather than through any distance metric: the penalty discourages the fitted values of neighboring regions from differing, so choose this basis over TPRS-style continuous smooths whenever the domain is a set of discrete areas linked by an adjacency structure (e.g. shared borders) instead of by Euclidean coordinates.
k: int = -1Maximum number of regions to retain. If -1 (the default), all observed levels of the grouping variable are kept as basis functions. If a positive integer smaller than the number of observed levels, only the first k (in sorted level order) are used; this is rarely what a user wants for MRF smooths (unlike continuous bases, reducing k does not give a lower-rank approximation of the same structure — it silently drops regions), so in most workflows the default of -1 should be left alone.
neighborhood: dict | NDArray | None = Nonedict mapping region labels to lists of neighbor labels (only pairs need to be listed once; the adjacency is symmetrized automatically), or a square, symmetric adjacency matrix (ndarray) whose row/column order matches the sorted unique levels of the fitted grouping variable. This argument is required — there is no sensible default neighborhood structure.
Let there be k unique regions after fit(). The basis matrix B is the n x k matrix of region indicators, B[i, j] = 1 if observation i belongs to region j and 0 otherwise — identical in structure to RandomEffectBasis. What distinguishes an MRF smooth is its penalty. Writing A for the symmetric adjacency matrix (A[i, j] = 1 if regions i and j are neighbors) and D = \operatorname{diag}(A \mathbf{1}) for the diagonal matrix of neighbor counts, the penalty matrix is the graph Laplacian
\mathbf{L} = \mathbf{D} - \mathbf{A},
so that the roughness penalty takes the form
\boldsymbol{\beta}^\top \mathbf{L} \boldsymbol{\beta} = \sum_{(i,j) \, \in \, \text{neighbors}} (\beta_i - \beta_j)^2 .
This penalizes exactly the pairwise differences between the fitted level for each region and the fitted levels of its geographic neighbors, pulling adjacent regions toward a common value as lambda grows, while leaving regions that are far apart on the map free to differ. The graph Laplacian of a connected neighborhood graph has exactly one zero eigenvalue, with eigenvector proportional to the all-ones vector; null_space_dimension() therefore returns 1 for a fully connected graph (the penalty cannot shrink a common overall level, only differences between neighbors), but can be larger if the neighborhood graph has multiple disconnected components, since each component then has its own unpenalized constant. Because the raw indicator basis shares an unpenalized constant with the model intercept, a sum-to-zero constraint (returned by identifiability_constraints()) is needed for identifiability when fitting alongside an intercept term. fit() raises if fewer than two regions are present, since a spatial smooth is meaningless with only one area, and raises if no neighborhood is supplied.
import numpy as np
from whittaker.smooths.mrf import MRFBasis
rng = np.random.default_rng(0)
regions = np.array(["A", "B", "C", "D"])
x = rng.choice(regions, size=40)
neighborhood = {
"A": ["B"],
"B": ["A", "C"],
"C": ["B", "D"],
"D": ["C"],
}
basis = MRFBasis(neighborhood=neighborhood).fit(x)
B = basis.basis_matrix(x)
S = basis.penalty_matrix()
B.shape, S.shape((40, 4), (4, 4))| Name | Description |
|---|---|
| k | Requested or actual number of region levels. |
| levels |
Sorted array of unique region labels retained during fit().
|
| n_basis | Number of basis functions, i.e. the number of retained region levels. |
Requested or actual number of region levels.
k: int
Before fit() is called, returns the k value passed at construction (the requested cap on the number of levels, or -1 for “all levels”). After fit(), returns the actual number of region levels retained.
Sorted array of unique region labels retained during fit().
levels: NDArray
Number of basis functions, i.e. the number of retained region levels.
n_basis: int
RuntimeErrorfit() has been called.
| Name | Description |
|---|---|
| basis_matrix() |
Evaluate the MRF indicator basis at x.
|
| fit() | Fit the MRF basis to training data. |
| identifiability_constraints() | Return the sum-to-zero constraint row for the intercept. |
| null_space_dimension() | Return the number of zero eigenvalues of the graph Laplacian. |
| penalty_matrix() |
Return the (k, k) graph Laplacian penalty matrix L = D - A.
|
Evaluate the MRF indicator basis at x.
Usage
x: NDArray(n,). Values not seen during fit() produce an all-zero row (no region indicator is set).
NDArray(n, k) where k is the number of region levels retained during fit(). Each row has at most one entry equal to 1.0.
Fit the MRF basis to training data.
Usage
Determines the set of unique region levels present in x, builds the symmetric adjacency matrix from the supplied neighborhood structure, and computes the graph Laplacian penalty.
x: NDArray(n,), of any hashable dtype (strings, integers, etc.).
MRFBasisself for method chaining.
ValueErrorIf fewer than 2 unique levels are present in x, if neighborhood was not supplied, if an adjacency matrix is provided with the wrong shape, or if k is requested but is less than 2.
TypeErrorneighborhood is neither a dict nor an ndarray.
Return the sum-to-zero constraint row for the intercept.
Usage
NDArray(1, k) matrix of equal weights 1 / k whose product with the coefficient vector forces the mean fitted region effect to be zero, resolving the confound between the MRF smooth’s unpenalized constant and the model intercept.
Return the number of zero eigenvalues of the graph Laplacian.
Usage
This equals the number of connected components of the neighborhood graph: 1 for a fully connected map, and more than 1 if some regions have no path of neighbors linking them to the rest, since each disconnected component then carries its own unpenalized constant.
intReturn the (k, k) graph Laplacian penalty matrix L = D - A.
Usage
NDArray(k, k), where k is the number of region levels.