# MRFBasis Markov random field basis for areal spatial data. Usage ``` python MRFBasis( k=-1, neighborhood=None, ) ``` 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](TPRS.md#whittaker.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. ## Parameters `k: int = ``-1` Maximum 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 = None` The neighborhood structure that defines which regions are considered adjacent. Either a `dict` 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. ## Notes 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](RandomEffectBasis.md#whittaker.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. ## Examples ``` python 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)) ## Attributes | Name | Description | |----|----| | [k](#k) | Requested or actual number of region levels. | | [levels](#levels) | Sorted array of unique region labels retained during `fit()`. | | [n_basis](#n_basis) | Number of basis functions, i.e. the number of retained region levels. | ------------------------------------------------------------------------ ### k 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. ------------------------------------------------------------------------ ### levels Sorted array of unique region labels retained during `fit()`. `levels: NDArray` ------------------------------------------------------------------------ ### n_basis Number of basis functions, i.e. the number of retained region levels. `n_basis: int` #### Raises `RuntimeError` If accessed before `fit()` has been called. ## Methods | Name | Description | |----|----| | [basis_matrix()](#basis_matrix) | Evaluate the MRF indicator basis at `x`. | | [fit()](#fit) | Fit the MRF basis to training data. | | [identifiability_constraints()](#identifiability_constraints) | Return the sum-to-zero constraint row for the intercept. | | [null_space_dimension()](#null_space_dimension) | Return the number of zero eigenvalues of the graph Laplacian. | | [penalty_matrix()](#penalty_matrix) | Return the `(k, k)` graph Laplacian penalty matrix `L = D - A`. | ------------------------------------------------------------------------ ### basis_matrix() Evaluate the MRF indicator basis at `x`. Usage ``` python basis_matrix(x) ``` #### Parameters `x: NDArray` Grouping variable values. Shape `(n,)`. Values not seen during `fit()` produce an all-zero row (no region indicator is set). #### Returns `NDArray` Design matrix of shape `(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() Fit the MRF basis to training data. Usage ``` python fit(x) ``` 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. #### Parameters `x: NDArray` Training grouping variable (region labels). Shape `(n,)`, of any hashable dtype (strings, integers, etc.). #### Returns `MRFBasis` Returns `self` for method chaining. #### Raises `ValueError` If 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. `TypeError` If `neighborhood` is neither a `dict` nor an `ndarray`. ------------------------------------------------------------------------ ### identifiability_constraints() Return the sum-to-zero constraint row for the intercept. Usage ``` python identifiability_constraints() ``` #### Returns `NDArray` A `(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. ------------------------------------------------------------------------ ### null_space_dimension() Return the number of zero eigenvalues of the graph Laplacian. Usage ``` python null_space_dimension() ``` 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. #### Returns `int` Dimension of the penalty null space. ------------------------------------------------------------------------ ### penalty_matrix() Return the `(k, k)` graph Laplacian penalty matrix `L = D - A`. Usage ``` python penalty_matrix() ``` #### Returns `NDArray` Symmetric positive semi-definite matrix of shape `(k, k)`, where `k` is the number of region levels.