diff --git a/README.md b/README.md
index 26ea9c57..5ddb40de 100644
--- a/README.md
+++ b/README.md
@@ -126,8 +126,8 @@ result, and prediction interfaces.
 | Formula surface    | Supported terms                                                                                      |
 | ------------------ | ---------------------------------------------------------------------------------------------------- |
 | Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `tp`, `ts`                              |
-| Structured smooths | random effects `re`, factor smooths `fs`, sum-to-zero factor smooths `sz`                            |
-| Tensor products    | `te(...)` and `ti(...)` over supported numeric marginals                                             |
+| Structured smooths | Markov random fields `mrf`, random effects `re`, factor smooths `fs`, sum-to-zero factor smooths `sz` |
+| Tensor products    | `te(...)` and `ti(...)` over supported numeric or MRF marginals                                      |
 | Parametric terms   | numeric and factor terms, supported interactions, intercept policies, and formula offsets            |
 | Shared smoothing   | supported `id=` groups, fixed/free smoothing parameters, `select=True`, and `pc=` on supported bases |
 
diff --git a/docs/generate_notebooks.py b/docs/generate_notebooks.py
index ea428518..332025d9 100644
--- a/docs/generate_notebooks.py
+++ b/docs/generate_notebooks.py
@@ -1515,6 +1515,7 @@ def gam_notebook() -> dict:
 | `cc` | cyclic cubic spline for periodic covariates |
 | `ps` | P-spline with difference penalties |
 | `gp` | low-rank Gaussian-process smooth with spherical, power-exponential, or Matérn covariance |
+| `mrf` | region effect coupled by a neighbor graph, polygon boundary, or supplied penalty |
 | `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |
 | `te(...)` | scale-invariant tensor product including main-effect directions |
 | `ti(...)` | tensor interaction with marginal main-effect directions removed |
@@ -1539,6 +1540,12 @@ def gam_notebook() -> dict:
     "gaussian_process": GAM(
         formula="demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])"
     ),
+    "markov_random_field": GAM(
+        formula=(
+            "demand ~ s(region, bs='mrf', "
+            "xt={'nb': {'north': ['south'], 'south': ['north']}})"
+        )
+    ),
     "thin_plate": GAM(formula="demand ~ s(temperature, bs='tp', k=10)"),
     "shrinkage_thin_plate": GAM(formula="demand ~ s(temperature, bs='ts', k=10)"),
     "tensor_surface": GAM(
diff --git a/docs/notebooks/01_gam.ipynb b/docs/notebooks/01_gam.ipynb
index 8c93cec3..6280af49 100644
--- a/docs/notebooks/01_gam.ipynb
+++ b/docs/notebooks/01_gam.ipynb
@@ -630,6 +630,7 @@
     "| `cc` | cyclic cubic spline for periodic covariates |\n",
     "| `ps` | P-spline with difference penalties |\n",
     "| `gp` | low-rank Gaussian-process smooth with spherical, power-exponential, or Mat\u00e9rn covariance |\n",
+    "| `mrf` | region effect coupled by a neighbor graph, polygon boundary, or supplied penalty |\n",
     "| `tp`, `ts` | thin-plate regression spline; `ts` adds shrinkage |\n",
     "| `te(...)` | scale-invariant tensor product including main-effect directions |\n",
     "| `ti(...)` | tensor interaction with marginal main-effect directions removed |\n",
@@ -654,6 +655,7 @@
        " 'cyclic',\n",
        " 'factor_smooth',\n",
        " 'gaussian_process',\n",
+       " 'markov_random_field',\n",
        " 'p_spline',\n",
        " 'random_effect',\n",
        " 'shrinkage_cubic',\n",
@@ -680,6 +682,12 @@
     "    \"gaussian_process\": GAM(\n",
     "        formula=\"demand ~ s(temperature, humidity, bs='gp', k=20, m=[3, 1.0])\"\n",
     "    ),\n",
+    "    \"markov_random_field\": GAM(\n",
+    "        formula=(\n",
+    "            \"demand ~ s(region, bs='mrf', \"\n",
+    "            \"xt={'nb': {'north': ['south'], 'south': ['north']}})\"\n",
+    "        )\n",
+    "    ),\n",
     "    \"thin_plate\": GAM(formula=\"demand ~ s(temperature, bs='tp', k=10)\"),\n",
     "    \"shrinkage_thin_plate\": GAM(formula=\"demand ~ s(temperature, bs='ts', k=10)\"),\n",
     "    \"tensor_surface\": GAM(\n",
diff --git a/nampy/gam/compiler/compile_model.py b/nampy/gam/compiler/compile_model.py
index 17bfd7fc..470367ca 100644
--- a/nampy/gam/compiler/compile_model.py
+++ b/nampy/gam/compiler/compile_model.py
@@ -36,10 +36,10 @@ def _apply_overlapping_parametric_identifiability(
     n_obs = int(compiled_predictors[0].design_matrix.shape[0])
     expanded_columns = []
     owners: list[tuple[int, int | None, int | None]] = []
-    keep_by_component: list[dict[int, np.ndarray]] = [
-        {} for _ in compiled_predictors
+    keep_by_component: list[dict[int, np.ndarray]] = [{} for _ in compiled_predictors]
+    keep_intercept = [
+        bool(predictor.has_intercept) for predictor in compiled_predictors
     ]
-    keep_intercept = [bool(predictor.has_intercept) for predictor in compiled_predictors]
 
     def append_column(column, targets, owner):
         expanded = np.zeros(n_obs * n_linear_predictors, dtype=np.float64)
@@ -101,9 +101,7 @@ def append_column(column, targets, owner):
         dropped_local: list[int] = []
         for term_index, term in enumerate(predictor.compiled_terms):
             basis = np.asarray(term.basis_train, dtype=np.float64)
-            keep = component_keep.get(
-                term_index, np.ones(basis.shape[1], dtype=bool)
-            )
+            keep = component_keep.get(term_index, np.ones(basis.shape[1], dtype=bool))
             kept_indices = np.flatnonzero(keep)
             selection = np.eye(basis.shape[1], dtype=np.float64)[:, kept_indices]
             metadata = dict(getattr(term, "metadata", {}) or {})
@@ -199,7 +197,11 @@ def _full_predictor_matrix(predictor, X: np.ndarray) -> tuple[np.ndarray, np.nda
     Z_fit = np.asarray(predictor.design_matrix, dtype=np.float64)
     pred_blocks = []
     for term in predictor.compiled_terms:
-        use_raw = bool(getattr(term, "metadata", {}).get("expose_raw_prediction_basis"))
+        term_metadata = dict(getattr(term, "metadata", {}) or {})
+        use_raw = bool(
+            term_metadata.get("expose_raw_prediction_basis")
+            or term_metadata.get("prediction_basis_map") is not None
+        )
         if use_raw:
             block = np.asarray(
                 term.prediction_parameterization_matrix(X), dtype=np.float64
@@ -250,7 +252,7 @@ def _fit_to_prediction_parameterization_map(
     # qr(Xp, LAPACK=TRUE) -> Rrank(R) -> triangular solve on QtX -> restore pivots.
     Q, R, piv = scipy_qr(X_pred, mode="economic", pivoting=True)
     p_pred = int(R.shape[1])
-    rank = upper_triangular_rrank(R, tol=float(np.finfo(np.float64).eps**0.9))
+    rank = upper_triangular_rrank(R, tol=float(np.finfo(np.float64).eps ** 0.9))
     QtX = np.asarray(Q.T @ X_fit, dtype=np.float64)[:rank, :]
     if rank < p_pred:
         R1 = np.asarray(R[:rank, :], dtype=np.float64)
@@ -342,9 +344,7 @@ def compile_model(
         tuple(
             int(value) - 1
             for value in (
-                (getattr(spec, "metadata", {}) or {}).get(
-                    "lpi", (component_index + 1,)
-                )
+                (getattr(spec, "metadata", {}) or {}).get("lpi", (component_index + 1,))
                 or (component_index + 1,)
             )
         )
@@ -352,9 +352,9 @@ def compile_model(
         else (component_index,)
         for component_index, spec in enumerate(predictor_specs)
     )
-    n_linear_predictors = max(
-        (max(indices) for indices in component_lpi if indices), default=0
-    ) + 1
+    n_linear_predictors = (
+        max((max(indices) for indices in component_lpi if indices), default=0) + 1
+    )
 
     compiled_predictors = compile_predictors(
         X=X,
@@ -452,9 +452,7 @@ def compile_model(
                     (
                         np.zeros(int(predictor.n_coef), dtype=bool)
                         if predictor.positive_coefficient_mask is None
-                        else np.asarray(
-                            predictor.positive_coefficient_mask, dtype=bool
-                        )
+                        else np.asarray(predictor.positive_coefficient_mask, dtype=bool)
                     ),
                 ]
             ),
diff --git a/nampy/gam/compiler/factory.py b/nampy/gam/compiler/factory.py
index 8cfa2bfe..122f2479 100644
--- a/nampy/gam/compiler/factory.py
+++ b/nampy/gam/compiler/factory.py
@@ -10,6 +10,7 @@
     FSmoothInteractionTerm,
     SZSmoothInteractionTerm,
 )
+from ..smooths.categorical.mrf import MarkovRandomFieldTerm
 from ..smooths.categorical.re import RandomEffectTerm
 from ..smooths.parametric import LinearTerm
 from ..smooths.registry import make_smooth_term
@@ -29,6 +30,7 @@
     DuchonSplineSmoothSpec,
     FactorSmoothInteractionSpec,
     GaussianProcessSmoothSpec,
+    MarkovRandomFieldSmoothSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
@@ -228,6 +230,23 @@ def instantiate_term(term_like: TermSpec | Any):
             metadata=metadata,
         )
 
+    if isinstance(smooth_spec, MarkovRandomFieldSmoothSpec):
+        return MarkovRandomFieldTerm(
+            feature=features,
+            k=smooth_spec.k,
+            label=label,
+            term_id=term_like.term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=smooth_spec.sp,
+            select=smooth_spec.select,
+            fixed=smooth_spec.fx,
+            constraint_mode=smooth_spec.constraint_mode,
+            knots=smooth_spec.knots,
+            xt=smooth_spec.xt,
+            metadata=metadata,
+        )
+
     if isinstance(smooth_spec, ShapeConstrainedSmoothSpec):
         if len(features) == 2:
             return BivariateShapePSplineTerm(
diff --git a/nampy/gam/compiler/linked_basis.py b/nampy/gam/compiler/linked_basis.py
index ea9aff0e..ea5da12b 100644
--- a/nampy/gam/compiler/linked_basis.py
+++ b/nampy/gam/compiler/linked_basis.py
@@ -103,6 +103,26 @@ def attach_shared_basis_metadata(predictor_specs, X, feature_names):
                 "pc= is not supported across id-linked s() terms with different "
                 "feature sets; mgcv 1.9-4 fails while constructing the shared basis."
             )
+        has_mrf = any(
+            term.smooth_spec is not None
+            and (
+                str(getattr(term.smooth_spec, "bs", "")).lower() == "mrf"
+                or (
+                    isinstance(getattr(term.smooth_spec, "bs", None), (list, tuple))
+                    and "mrf"
+                    in {
+                        str(value).lower()
+                        for value in getattr(term.smooth_spec, "bs", ())
+                    }
+                )
+            )
+            for term in group_terms
+        )
+        if has_mrf and len(feature_tuples) > 1:
+            raise NotImplementedError(
+                "id= is not supported across MRF terms with different feature "
+                "sets; mgcv 1.9-4 loses the factor topology while pooling them."
+            )
         for term in group_terms[1:]:
             _clone_linked_smooth_spec(base_term, term)
 
diff --git a/nampy/gam/model/api.py b/nampy/gam/model/api.py
index c9e9ba9a..b658c333 100644
--- a/nampy/gam/model/api.py
+++ b/nampy/gam/model/api.py
@@ -174,6 +174,7 @@ def _restore_prediction_na_rows(value, retained_rows, n_rows):
         "covariance",
         "select",
         "knots",
+        "xt",
         "min_sp",
         "drop_intercept",
         "formula",
@@ -240,6 +241,7 @@ def __init__(
         self.covariance = str(self.hparams.get("covariance", "bayes")).lower()
         self.select = bool(self.hparams.get("select", False))
         self.knots = self.hparams.get("knots", None)
+        self.xt = self.hparams.get("xt", None)
         self.min_sp = self.hparams.get("min_sp", None)
         self.drop_intercept = self.hparams.get("drop_intercept", None)
         self.positive_transform = str(
diff --git a/nampy/gam/smooths/__init__.py b/nampy/gam/smooths/__init__.py
index 07ab6dae..65c16552 100644
--- a/nampy/gam/smooths/__init__.py
+++ b/nampy/gam/smooths/__init__.py
@@ -1,4 +1,5 @@
 from .categorical.fs import FSmoothInteractionTerm, SZSmoothInteractionTerm
+from .categorical.mrf import MarkovRandomFieldTerm
 from .categorical.re import RandomEffectTerm
 from .registry import make_smooth_term, register_smooth
 from .shape.scop import ShapeConstrainedPSplineTerm
@@ -35,6 +36,7 @@
 cr = cs = cc = CubicSplineTerm
 ds = DuchonSplineTerm
 gp = GaussianProcessTerm
+mrf = MarkovRandomFieldTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 fs = FSmoothInteractionTerm
@@ -63,6 +65,7 @@
     "CubicSplineTerm",
     "DuchonSplineTerm",
     "GaussianProcessTerm",
+    "MarkovRandomFieldTerm",
     "DerivativeBSplineTerm1D",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
@@ -80,6 +83,7 @@
     "cc",
     "ds",
     "gp",
+    "mrf",
     "cp",
     "ps",
     "tp",
diff --git a/nampy/gam/smooths/categorical/__init__.py b/nampy/gam/smooths/categorical/__init__.py
index a446930f..8479c011 100644
--- a/nampy/gam/smooths/categorical/__init__.py
+++ b/nampy/gam/smooths/categorical/__init__.py
@@ -8,11 +8,13 @@
     try_numeric_1d,
 )
 from .fs import FSmoothInteractionTerm, SZSmoothInteractionTerm
+from .mrf import MarkovRandomFieldTerm
 from .re import RandomEffectTerm
 
 fs = FSmoothInteractionTerm
 sz = SZSmoothInteractionTerm
 re = RandomEffectTerm
+mrf = MarkovRandomFieldTerm
 
 __all__ = [
     "as_object_1d",
@@ -23,9 +25,11 @@
     "factor_indicator_matrix",
     "factor_levels_from_metadata",
     "RandomEffectTerm",
+    "MarkovRandomFieldTerm",
     "FSmoothInteractionTerm",
     "SZSmoothInteractionTerm",
     "fs",
     "sz",
     "re",
+    "mrf",
 ]
diff --git a/nampy/gam/smooths/categorical/fs.py b/nampy/gam/smooths/categorical/fs.py
index 484dcb63..6dd0a113 100644
--- a/nampy/gam/smooths/categorical/fs.py
+++ b/nampy/gam/smooths/categorical/fs.py
@@ -27,6 +27,7 @@
     is_factor_like_vector,
     stable_unique_levels,
 )
+from .mrf import MarkovRandomFieldTerm
 from .re import RandomEffectTerm
 
 
@@ -110,7 +111,7 @@ def _build_base_smooth_term(
     Build the per-level base smooth used inside fs/sz.
 
     Supported base smooth classes in the current codebase:
-    bs, cr, cs, cc, cp, ds, gp, ps, tp, ts
+    bs, cr, cs, cc, cp, ds, gp, mrf, ps, tp, ts
     """
     base_bs = str(base_bs).lower()
     metric_features = list(metric_features)
@@ -132,12 +133,13 @@ def _build_base_smooth_term(
         "cp",
         "ds",
         "gp",
+        "mrf",
         "ps",
         "tp",
         "ts",
     }:
         raise NotImplementedError(
-            "Extra xt options are currently only supported for bs/cp/ds/gp/ps/tp/ts "
+            "Extra xt options are currently only supported for bs/cp/ds/gp/mrf/ps/tp/ts "
             "base smooths, "
             f"got xt={xt_rest!r} with base bs={base_bs!r}."
         )
@@ -240,6 +242,22 @@ def _build_base_smooth_term(
             metadata=metadata,
         )
 
+    if base_bs == "mrf":
+        return MarkovRandomFieldTerm(
+            feature=metric_features[0],
+            k=k,
+            label=label,
+            smoothing_id=None,
+            by=by,
+            sp=None,
+            select=bool(select),
+            fixed=bool(fixed),
+            constraint_mode=str(constraint_mode),
+            knots=knots,
+            xt=xt_rest,
+            metadata=metadata,
+        )
+
     if base_bs in {"tp", "ts"}:
         return make_smooth_term(
             base_bs,
@@ -262,7 +280,7 @@ def _build_base_smooth_term(
 
     raise NotImplementedError(
         f"Current {mode} implementation supports base bs in "
-        f"{{'bs','cr','cs','cc','cp','ds','gp','ps','tp','ts'}}, got {base_bs!r}."
+        f"{{'bs','cr','cs','cc','cp','ds','gp','mrf','ps','tp','ts'}}, got {base_bs!r}."
     )
 
 
@@ -273,6 +291,8 @@ def _penalty_rank_from_base_term(base_term, basis_matrix, penalty_matrix) -> int
         return int(base_term._setup.rank)
     if isinstance(base_term, GaussianProcessTerm):
         return int(base_term._setup.rank)
+    if isinstance(base_term, MarkovRandomFieldTerm):
+        return int(base_term._setup.rank)
     if isinstance(base_term, PSplineTerm1D) and len(base_term.penalties) > 0:
         if str(base_term.basis_name).lower() == "cp":
             return int(base_term._setup.rank)
@@ -668,6 +688,8 @@ def fit(self, X, feature_names):
             metadata=dict(self.metadata),
         )
         base_term.fit(X, feature_names)
+        if isinstance(base_term, MarkovRandomFieldTerm):
+            self.metadata = dict(base_term.metadata)
 
         if len(base_term.penalties) > 1:
             raise NotImplementedError(
@@ -675,6 +697,15 @@ def fit(self, X, feature_names):
             )
 
         self._base_term = base_term
+        if (
+            isinstance(base_term, MarkovRandomFieldTerm)
+            and base_term._setup.used_low_rank
+        ):
+            raise NotImplementedError(
+                "mgcv 1.9-4 cannot predict an fs smooth with a reduced-rank "
+                "MRF base because its factor-smooth P matrix is dimensionally "
+                "incompatible with the full region indicator."
+            )
         B_setup, S0, _ = self._base_constructor_fit_matrices()
         B_setup = np.asarray(B_setup, dtype=np.float64)
         B0 = np.asarray(self._base_constructor_predict_matrix(X), dtype=np.float64)
@@ -722,6 +753,15 @@ def fit(self, X, feature_names):
         null_d = B0.shape[1] - r
 
         self._base_transform = P_coef
+        if isinstance(base_term, MarkovRandomFieldTerm):
+            # The raw prediction callback below deliberately exposes mgcv's
+            # double-P MRF/fs prediction surface. Tell the compiler which
+            # inverse block map recovers the constructor matrix solely for its
+            # fit-vs-prediction parameterization audit; actual prediction does
+            # not apply this metadata map.
+            self.metadata["prediction_basis_map"] = np.kron(
+                np.eye(n_levels, dtype=np.float64), np.linalg.inv(P_coef)
+            )
         self._base_range_penalty_diag = np.concatenate(
             [D, np.zeros(null_d, dtype=np.float64)]
         )
@@ -801,6 +841,13 @@ def transform_new(self, X_new):
         )
         if self._base_transform is not None:
             B0_new = B0_new @ self._base_transform
+            if isinstance(self._base_term, MarkovRandomFieldTerm):
+                # Upstream Predict.matrix.fs.interaction changes the object back
+                # to class mrf.smooth after overwriting object$P with the fs
+                # natural-parameter transform. Predict.matrix.mrf.smooth applies
+                # that P once and the fs wrapper applies it a second time. Preserve
+                # this observable mgcv 1.9-4 prediction parameterization exactly.
+                B0_new = B0_new @ self._base_transform
 
         B_new = rowwise_kronecker([Ifac, B0_new])
         z = by_values_from_new_data(X_new, self._by_state)
@@ -948,6 +995,8 @@ def fit(self, X, feature_names):
             metadata=dict(self.metadata),
         )
         base_term.fit(X, feature_names)
+        if isinstance(base_term, MarkovRandomFieldTerm):
+            self.metadata = dict(base_term.metadata)
 
         if len(base_term.penalties) > 1:
             raise NotImplementedError(
diff --git a/nampy/gam/smooths/categorical/mrf.py b/nampy/gam/smooths/categorical/mrf.py
new file mode 100644
index 00000000..8360838d
--- /dev/null
+++ b/nampy/gam/smooths/categorical/mrf.py
@@ -0,0 +1,532 @@
+"""Markov-random-field smooths matching mgcv ``bs='mrf'``."""
+
+from __future__ import annotations
+
+import warnings
+from dataclasses import dataclass
+
+import numpy as np
+import pandas as pd
+
+from ...constraints.absorption import (
+    fit_single_penalty_with_constraint_policy,
+    fit_single_penalty_with_setup_basis,
+)
+from ...penalties.algebra import penalty_rescale_factor, scale_penalty
+from ...splines.basis.natparam import nat_param_type0
+from ..registry import register_smooth
+from ..smooth_base import BaseSmoothTerm, _resolve_feature, column_as_object
+from .categorical_utils import factor_levels_from_metadata
+
+
+def _canonical_label(value) -> str | None:
+    if pd.isna(value):
+        return None
+    if isinstance(value, (bool, np.bool_)):
+        return "TRUE" if bool(value) else "FALSE"
+    if isinstance(value, (int, np.integer)):
+        return str(int(value))
+    if isinstance(value, (float, np.floating)):
+        number = float(value)
+        if number.is_integer():
+            return str(int(number))
+        return format(number, ".15g")
+    return str(value)
+
+
+def _canonical_labels(values, *, boolean_coded=False) -> np.ndarray:
+    raw = np.asarray(values, dtype=object).ravel()
+    labels = []
+    for value in raw:
+        if (
+            boolean_coded
+            and isinstance(value, (int, np.integer, float, np.floating))
+            and not isinstance(value, (bool, np.bool_))
+            and float(value) in {0.0, 1.0}
+        ):
+            labels.append("TRUE" if float(value) == 1.0 else "FALSE")
+        else:
+            labels.append(_canonical_label(value))
+    return np.asarray(labels, dtype=object)
+
+
+def _sorted_factor_levels(values) -> list[str]:
+    raw = np.asarray(values, dtype=object).ravel()
+    raw = np.asarray([value for value in raw if not pd.isna(value)], dtype=object)
+    try:
+        ordered = np.unique(raw).tolist()
+    except (TypeError, ValueError):
+        ordered = sorted(set(_canonical_labels(raw)))
+    return [str(_canonical_label(value)) for value in ordered]
+
+
+def _knot_levels(knots) -> list[str]:
+    if isinstance(knots, pd.Categorical):
+        return [str(_canonical_label(value)) for value in knots.categories]
+    if isinstance(knots, pd.Series) and isinstance(knots.dtype, pd.CategoricalDtype):
+        return [str(_canonical_label(value)) for value in knots.cat.categories]
+    return _sorted_factor_levels(np.asarray(knots, dtype=object))
+
+
+def _indicator_matrix(values, levels) -> np.ndarray:
+    labels = _canonical_labels(values)
+    index = {str(level): j for j, level in enumerate(levels)}
+    out = np.zeros((labels.size, len(levels)), dtype=np.float64)
+    for row, label in enumerate(labels):
+        if label is not None and label in index:
+            out[row, index[label]] = 1.0
+    return out
+
+
+def _polygons_to_neighbors(polygons) -> dict[str, list[str]]:
+    if not isinstance(polygons, dict):
+        raise TypeError(
+            "MRF xt['polys'] must be a mapping from area names to polygons."
+        )
+    names = [str(_canonical_label(name)) for name in polygons]
+    vertices: dict[str, set[tuple[float, float]]] = {}
+    for raw_name, polygon in polygons.items():
+        array = np.asarray(polygon, dtype=np.float64)
+        if array.ndim != 2 or array.shape[1] != 2:
+            raise ValueError("Each MRF polygon must be a two-column matrix.")
+        finite = array[np.isfinite(array).all(axis=1)]
+        vertices[str(_canonical_label(raw_name))] = {
+            tuple(map(float, row)) for row in finite
+        }
+    return {
+        name: [
+            other
+            for other in names
+            if other != name and vertices[name] & vertices[other]
+        ]
+        for name in names
+    }
+
+
+def _neighbor_penalty(neighbors, levels) -> np.ndarray:
+    if not isinstance(neighbors, dict) or len(neighbors) == 0:
+        raise TypeError("MRF xt['nb'] must be a non-empty named mapping.")
+    area_names = [str(_canonical_label(name)) for name in neighbors]
+    if sorted(area_names) != sorted(levels):
+        raise ValueError(
+            "mismatch between nb/polys supplied area names and data area names"
+        )
+
+    entry_modes = set()
+    for raw_adjacent in neighbors.values():
+        adjacent = np.asarray(raw_adjacent, dtype=object).ravel().tolist()
+        if not adjacent:
+            continue
+        numeric_flags = [
+            isinstance(value, (int, np.integer, float, np.floating))
+            and not isinstance(value, (bool, np.bool_))
+            and float(value).is_integer()
+            for value in adjacent
+        ]
+        if any(numeric_flags) and not all(numeric_flags):
+            raise TypeError(
+                "MRF neighbour entries must use either all area names or all "
+                "one-based numeric indices."
+            )
+        entry_modes.add("numeric" if all(numeric_flags) else "named")
+    if len(entry_modes) > 1:
+        raise TypeError(
+            "MRF neighbour lists must use one uniform representation: either "
+            "area names or one-based numeric indices."
+        )
+    neighbor_mode = next(iter(entry_modes), "named")
+
+    S = np.zeros((len(levels), len(levels)), dtype=np.float64)
+    level_index = {name: i for i, name in enumerate(levels)}
+    for raw_name, raw_adjacent in neighbors.items():
+        name = str(_canonical_label(raw_name))
+        adjacent = np.asarray(raw_adjacent, dtype=object).ravel().tolist()
+        if neighbor_mode == "numeric":
+            adjacent_names = []
+            for value in adjacent:
+                position = int(value) - 1
+                if position < 0 or position >= len(area_names):
+                    raise IndexError("MRF neighbour index is outside the named list.")
+                adjacent_names.append(area_names[position])
+        else:
+            requested = {str(_canonical_label(value)) for value in adjacent}
+            # mgcv maps named neighbours through `which(nb.names %in% entry)`,
+            # so names are set-like: ordering, duplicates, and unknown names
+            # disappear before the degree is computed.
+            adjacent_names = [value for value in area_names if value in requested]
+
+        row = level_index[name]
+        S[row, row] = float(len(adjacent_names))
+        for adjacent_name in adjacent_names:
+            column = level_index[adjacent_name]
+            if column != row:
+                S[row, column] = -1.0
+    if np.any(S != S.T):
+        raise ValueError("Something wrong with auto- penalty construction")
+    return S
+
+
+def _supplied_penalty(penalty, levels) -> np.ndarray:
+    if isinstance(penalty, pd.DataFrame):
+        column_names = [str(_canonical_label(value)) for value in penalty.columns]
+        if sorted(column_names) != sorted(levels):
+            raise ValueError("penalty column names don't match supplied area names!")
+        lookup_rows = {str(_canonical_label(value)): value for value in penalty.index}
+        lookup_cols = {str(_canonical_label(value)): value for value in penalty.columns}
+        penalty = penalty.loc[
+            [lookup_rows[level] for level in levels],
+            [lookup_cols[level] for level in levels],
+        ]
+    matrix = np.asarray(penalty, dtype=np.float64)
+    if matrix.ndim != 2 or matrix.shape[0] != matrix.shape[1]:
+        raise ValueError("supplied penalty not square!")
+    if matrix.shape != (len(levels), len(levels)):
+        raise ValueError("supplied penalty wrong dimension!")
+    return matrix
+
+
+@dataclass
+class MarkovRandomFieldSetup:
+    levels: tuple[str, ...]
+    basis_train: np.ndarray
+    penalty: np.ndarray
+    raw_penalty: np.ndarray
+    P: np.ndarray | None
+    rank: int
+    null_space_dim: int
+    bs_dim: int
+    used_low_rank: bool
+    plot_me: bool
+    boolean_coded: bool
+
+
+def build_markov_random_field_setup(
+    values, *, k=-1, xt=None, knots=None, factor_levels=None
+) -> MarkovRandomFieldSetup:
+    """Port ``smooth.construct.mrf.smooth.spec`` before smoothCon scaling."""
+    raw_values = np.asarray(values, dtype=object).ravel()
+    declared = (
+        []
+        if factor_levels is None
+        else np.asarray(factor_levels, dtype=object).ravel().tolist()
+    )
+    boolean_coded = bool(
+        declared and all(isinstance(value, (bool, np.bool_)) for value in declared)
+    )
+    labels = _canonical_labels(raw_values, boolean_coded=boolean_coded)
+    if any(value is None for value in labels):
+        raise ValueError("MRF smooths do not allow missing regions in fitting.")
+
+    if factor_levels is None:
+        data_levels = _sorted_factor_levels(labels)
+    else:
+        data_levels = sorted(str(_canonical_label(value)) for value in factor_levels)
+    levels = data_levels if knots is None else _knot_levels(knots)
+    if any(value not in levels for value in data_levels):
+        raise ValueError(
+            "data contain regions that are not contained in the knot specification"
+        )
+    observed_levels = {value for value in labels if value is not None}
+    if any(value not in levels for value in observed_levels):
+        raise ValueError(
+            "data contain regions that are not contained in the knot specification"
+        )
+
+    bs_dim = len(levels) if int(k) < 0 else int(k)
+    if bs_dim > len(levels):
+        raise ValueError("MRF basis dimension set too high")
+    if bs_dim <= 2 and bs_dim < len(levels):
+        raise ValueError(
+            "A reduced-rank MRF basis with k<=2 is malformed in mgcv 1.9-4; use k>=3."
+        )
+    if xt is None:
+        raise ValueError(
+            "penalty matrix, boundary polygons and/or neighbours list must be supplied in xt"
+        )
+    if not isinstance(xt, dict):
+        raise TypeError("MRF xt must be a dictionary.")
+
+    plot_me = xt.get("polys") is not None
+    if xt.get("penalty") is not None:
+        raw_penalty = _supplied_penalty(xt["penalty"], levels)
+    else:
+        neighbors = xt.get("nb")
+        if neighbors is None:
+            polygons = xt.get("polys")
+            if polygons is None:
+                raise ValueError("no spatial information provided!")
+            neighbors = _polygons_to_neighbors(polygons)
+        raw_penalty = _neighbor_penalty(neighbors, levels)
+
+    X = _indicator_matrix(labels, levels)
+    P = None
+    if bs_dim < len(levels):
+        missing = np.flatnonzero(np.sum(X, axis=0) == 0.0)
+        X_complete = X
+        if missing.size:
+            dummy = np.zeros((missing.size, X.shape[1]), dtype=np.float64)
+            dummy[np.arange(missing.size), missing] = 1.0
+            X_complete = np.vstack([dummy, X])
+        natural = nat_param_type0(X_complete, raw_penalty)
+        # Retain the final (least penalized) natural-parameter columns.
+        keep = np.arange(len(levels) - bs_dim, len(levels), dtype=int)
+        X_natural = np.asarray(natural["X"], dtype=np.float64)
+        X = X_natural[missing.size :, keep] if missing.size else X_natural[:, keep]
+        P = np.asarray(natural["P"], dtype=np.float64)[:, keep]
+        diagonal = np.asarray(
+            [natural["D"][i] if i < natural["rank"] else 0.0 for i in keep],
+            dtype=np.float64,
+        )
+        penalty = np.diag(diagonal)
+        rank = int(np.sum(keep < natural["rank"]))
+    else:
+        penalty = np.asarray(raw_penalty, dtype=np.float64)
+        eigenvalues = np.linalg.eigvalsh(penalty)
+        largest = float(np.max(eigenvalues)) if eigenvalues.size else 0.0
+        rank = int(np.sum(eigenvalues > np.finfo(float).eps ** 0.8 * largest))
+
+    return MarkovRandomFieldSetup(
+        levels=tuple(levels),
+        basis_train=np.asarray(X, dtype=np.float64),
+        penalty=np.asarray(penalty, dtype=np.float64),
+        raw_penalty=np.asarray(raw_penalty, dtype=np.float64),
+        P=P,
+        rank=rank,
+        null_space_dim=int(bs_dim - rank),
+        bs_dim=int(bs_dim),
+        used_low_rank=bool(P is not None),
+        plot_me=bool(plot_me),
+        boolean_coded=boolean_coded,
+    )
+
+
+def predict_markov_random_field(values, setup: MarkovRandomFieldSetup) -> np.ndarray:
+    labels = _canonical_labels(values, boolean_coded=setup.boolean_coded)
+    unknown = sorted(
+        {
+            str(value)
+            for value in labels
+            if value is not None and value not in setup.levels
+        }
+    )
+    if unknown:
+        raise ValueError(
+            f"MRF prediction data contain unknown regions {unknown}; "
+            f"known regions are {list(setup.levels)}."
+        )
+    basis = _indicator_matrix(labels, setup.levels)
+    if setup.P is not None:
+        basis = basis @ setup.P
+    return np.asarray(basis, dtype=np.float64)
+
+
+@register_smooth("mrf")
+class MarkovRandomFieldTerm(BaseSmoothTerm):
+    term_type = "smooth"
+    basis_name = "mrf"
+    supports_tensor_marginal = True
+
+    def __init__(
+        self,
+        feature,
+        k=-1,
+        label=None,
+        term_id=None,
+        smoothing_id=None,
+        by=None,
+        sp=None,
+        select=False,
+        fixed=False,
+        constraint_mode="auto",
+        knots=None,
+        xt=None,
+        metadata=None,
+    ):
+        features = list(feature) if not isinstance(feature, (str, int)) else [feature]
+        if len(features) != 1:
+            raise ValueError("MRF smooths require exactly one region covariate.")
+        super().__init__(
+            feature=features,
+            label=label or f"s({features[0]})",
+            term_id=term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=sp,
+            metadata=metadata,
+        )
+        self.k = int(k)
+        self.select = bool(select)
+        self.fixed = bool(fixed)
+        self.constraint_mode = str(constraint_mode).lower()
+        self.knots = knots
+        self.xt = xt
+        self._feature_index = None
+        self._feature_name = None
+        self._factor_feature_indices = None
+        self._factor_feature_names = None
+        self._factor_levels = None
+        self._setup = None
+        self._basis_train = None
+        self._penalties = None
+
+    @property
+    def expected_linked_penalty_count(self):
+        return None if self.select else 1
+
+    def fit(self, X, feature_names):
+        if self.fixed and self.sp is not None:
+            sp_values = np.asarray(self.sp, dtype=np.float64).ravel()
+            if np.any(sp_values >= 0.0):
+                raise ValueError(
+                    "incorrect number of smoothing parameters supplied for a smooth term"
+                )
+        index, name = _resolve_feature(self.feature[0], feature_names)
+        values = column_as_object(X, index)
+        declared_levels = factor_levels_from_metadata(self.metadata, name)
+        if declared_levels is None:
+            try:
+                np.asarray(values, dtype=np.float64)
+            except (TypeError, ValueError):
+                pass
+            else:
+                warnings.warn(
+                    "argument of mrf should be a factor variable", stacklevel=2
+                )
+
+        self._set_by_state(X, feature_names)
+        self._feature_index = index
+        self._feature_name = name
+        self._factor_feature_indices = [index]
+        self._factor_feature_names = [name]
+        self._set_resolved_features([name])
+
+        shared_X = self._linked_id_setup_matrix(feature_names)
+        setup_values = values if shared_X is None else column_as_object(shared_X, index)
+        self._setup = build_markov_random_field_setup(
+            setup_values,
+            k=self.k,
+            xt=self.xt,
+            knots=self.knots,
+            factor_levels=(declared_levels if shared_X is None else None),
+        )
+        self._factor_levels = [list(self._setup.levels)]
+        if declared_levels is not None and self.knots is not None:
+            meta = dict(self.metadata or {})
+            factor_meta = dict(meta.get("factor_levels_by_feature", {}))
+            current = dict(factor_meta.get(str(name), {}))
+            allowed = list(current.get("levels", declared_levels))
+            for value in self._setup.levels:
+                if value not in allowed:
+                    allowed.append(value)
+            current["levels"] = allowed
+            factor_meta[str(name)] = current
+            meta["factor_levels_by_feature"] = factor_meta
+            self.metadata = meta
+        setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
+        base = (
+            setup_base
+            if shared_X is None
+            else predict_markov_random_field(values, self._setup)
+        )
+        raw_penalty = np.asarray(self._setup.penalty, dtype=np.float64)
+        penalty = scale_penalty(setup_base, raw_penalty)
+        self._set_penalty_rescale_factors(
+            [penalty_rescale_factor(setup_base, raw_penalty)]
+        )
+
+        auto_constrain = bool(self._by_state.is_constant)
+        if shared_X is None:
+            result = fit_single_penalty_with_constraint_policy(
+                base,
+                penalty,
+                self._by_state,
+                constraint_mode=self.constraint_mode,
+                fixed=self.fixed,
+                auto_constrain_when=auto_constrain,
+            )
+        else:
+            result = fit_single_penalty_with_setup_basis(
+                base,
+                setup_base,
+                penalty,
+                self._by_state,
+                constraint_mode=self.constraint_mode,
+                fixed=self.fixed,
+                auto_constrain_when=auto_constrain,
+            )
+        self._basis_train = result.basis_train
+        self._penalties = result.penalties
+        self._record_constraint_result(
+            result.constraint_kind,
+            result.constraint_transform,
+            absorbed_by="runtime" if result.constraint_transform is not None else None,
+        )
+        return self
+
+    def get_penalty_definitions(self):
+        self._require_fitted()
+        if not self.penalties:
+            return []
+        metadata = self._penalty_metadata_with_scale(
+            {
+                "term_type": self.term_type,
+                "basis_name": self.basis_name,
+                "feature": list(self.feature),
+                "label": self.label,
+                "by": self.by,
+                "by_name": self._by_state.feature_name,
+                "constraint_mode": self.constraint_mode,
+                "constraint_kind": self.constraint_kind,
+                "knots": self.knots,
+                "xt": self.xt,
+                "levels": list(self._setup.levels),
+                "fixed": self.fixed,
+            },
+            penalty_index=0,
+        )
+        return self._build_penalty_block(
+            self.penalties[0],
+            rank=min(int(self._setup.rank), int(self.penalties[0].shape[0])),
+            smooth_metadata=metadata,
+            selection_metadata={**metadata, "is_selection_penalty": True},
+        )
+
+    def transform_new(self, X_new):
+        self._require_fitted()
+        basis = predict_markov_random_field(
+            column_as_object(X_new, self._feature_index), self._setup
+        )
+        return self._apply_constraint_transform_and_by(basis, X_new)
+
+    def tensor_marginal_fit_matrices(
+        self, *, centered=False, apply_np=False, x_train=None
+    ):
+        del apply_np, x_train
+        self._require_fitted()
+        setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
+        setup_penalty = np.asarray(self._setup.penalty, dtype=np.float64)
+        if centered:
+            return super().tensor_marginal_fit_matrices(centered=True)
+        return setup_base, setup_penalty, None
+
+    def tensor_marginal_predict_matrix(
+        self, X_new, *, centered=False, np_transform=None
+    ):
+        basis = (
+            self.transform_new(X_new)
+            if centered
+            else predict_markov_random_field(
+                column_as_object(X_new, self._feature_index), self._setup
+            )
+        )
+        if np_transform is not None:
+            basis = basis @ np.asarray(np_transform, dtype=np.float64)
+        return np.asarray(basis, dtype=np.float64)
+
+
+__all__ = [
+    "MarkovRandomFieldSetup",
+    "MarkovRandomFieldTerm",
+    "build_markov_random_field_setup",
+    "predict_markov_random_field",
+]
diff --git a/nampy/gam/smooths/tensor/marginals.py b/nampy/gam/smooths/tensor/marginals.py
index afdf833d..cec7fb08 100644
--- a/nampy/gam/smooths/tensor/marginals.py
+++ b/nampy/gam/smooths/tensor/marginals.py
@@ -6,6 +6,7 @@
 
 from ...penalties.tensor import normalize_tensor_marginal_penalty
 from ..algebra import rowwise_kronecker
+from ..categorical.mrf import MarkovRandomFieldTerm
 from ..smooth_base import column_as_float
 from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
@@ -15,7 +16,7 @@
 from ..univariate.tp import ThinPlateSplineTerm
 
 TENSOR_MARGINAL_BASES = frozenset(
-    {"bs", "cr", "cs", "cc", "cp", "ds", "gp", "ps", "tp", "ts"}
+    {"bs", "cr", "cs", "cc", "cp", "ds", "gp", "mrf", "ps", "tp", "ts"}
 )
 
 
@@ -46,16 +47,15 @@ def make_tensor_marginal_term(
     knots=None,
     centered=False,
     shared_basis_setup=None,
+    metadata=None,
 ):
     basis = str(basis).lower()
     validate_tensor_marginal_bases([basis])
     constraint_mode = "always" if centered else "never"
     marginal_features = _as_marginal_features(feature)
-    metadata = (
-        None
-        if shared_basis_setup is None
-        else {"shared_basis_setup": shared_basis_setup}
-    )
+    metadata = dict(metadata or {})
+    if shared_basis_setup is not None:
+        metadata["shared_basis_setup"] = shared_basis_setup
     if basis in {"cr", "cs", "cc"}:
         if len(marginal_features) != 1:
             raise ValueError(
@@ -147,6 +147,23 @@ def make_tensor_marginal_term(
             metadata=metadata,
         )
 
+    if basis == "mrf":
+        if len(marginal_features) != 1:
+            raise ValueError("Tensor marginal basis 'mrf' only handles one feature.")
+        return MarkovRandomFieldTerm(
+            feature=marginal_features[0],
+            k=k,
+            xt=xt,
+            label=str(feature),
+            smoothing_id=None,
+            by=None,
+            select=False,
+            fixed=False,
+            constraint_mode=constraint_mode,
+            knots=knots,
+            metadata=metadata,
+        )
+
     if basis in {"tp", "ts"}:
         return ThinPlateSplineTerm(
             feature=marginal_features,
@@ -281,6 +298,7 @@ def build_tensor_marginal_terms(
     knots=None,
     centered=False,
     shared_basis_setups=None,
+    metadata=None,
 ):
     features = list(feature) if not isinstance(feature, (str, int)) else [feature]
     k_list = [int(k)] * len(features) if np.isscalar(k) else [int(v) for v in k]
@@ -343,6 +361,7 @@ def build_tensor_marginal_terms(
             knots=knots_i,
             centered=center_i,
             shared_basis_setup=shared_i,
+            metadata=metadata,
         )
         marginals.append(term)
         feature_ids.append(feat)
@@ -372,6 +391,7 @@ def build_tensor_product_components(
         x_train = (
             column_as_float(X, marginal_indices[0])
             if len(marginal_indices) == 1
+            and str(getattr(m, "basis_name", "")).lower() != "mrf"
             else None
         )
         shared_setup = getattr(m, "shared_basis_setup", None)
@@ -380,7 +400,11 @@ def build_tensor_product_components(
             and str(shared_setup.get("mode", "")).lower() == "linked_id"
             and shared_setup.get("pooled_feature_values")
         )
-        if use_linked_id_predict_path and len(marginal_indices) == 1:
+        if (
+            use_linked_id_predict_path
+            and len(marginal_indices) == 1
+            and str(getattr(m, "basis_name", "")).lower() != "mrf"
+        ):
             x_train = np.asarray(
                 shared_setup["pooled_feature_values"][0], dtype=np.float64
             ).ravel()
@@ -478,7 +502,7 @@ def _tensor_marginal_eval_from_x(term, x, *, centered=False):
 
 
 def _tensor_np_reparameterization(term, x_train, basis_dim, *, centered=False):
-    if str(getattr(term, "basis_name", "")).lower() in {"cr", "cs", "cc"}:
+    if str(getattr(term, "basis_name", "")).lower() in {"cr", "cs", "cc", "mrf"}:
         return None
     if x_train is None:
         return None
diff --git a/nampy/gam/smooths/tensor/te.py b/nampy/gam/smooths/tensor/te.py
index 146dbc56..ec45ebfa 100644
--- a/nampy/gam/smooths/tensor/te.py
+++ b/nampy/gam/smooths/tensor/te.py
@@ -131,9 +131,12 @@ def fit(self, X, feature_names):
             knots=self.knots,
             centered=False,
             shared_basis_setups=marginal_shared_setups,
+            metadata=self.metadata,
         )
         for term in marginals:
             term.fit(X, feature_names)
+            if str(getattr(term, "basis_name", "")).lower() == "mrf":
+                self.metadata = dict(term.metadata)
         feature_indices, feature_names_resolved = resolve_tensor_marginal_features(
             marginals
         )
@@ -223,9 +226,7 @@ def point_basis_fn(point):
         self._set_penalty_rescale_factors(
             [
                 float(scale)
-                for scale, keep in zip(
-                    penalty_scales, keep_penalties, strict=True
-                )
+                for scale, keep in zip(penalty_scales, keep_penalties, strict=True)
                 if keep
             ]
         )
diff --git a/nampy/gam/smooths/tensor/ti.py b/nampy/gam/smooths/tensor/ti.py
index 51877be7..de0ce16f 100644
--- a/nampy/gam/smooths/tensor/ti.py
+++ b/nampy/gam/smooths/tensor/ti.py
@@ -137,9 +137,12 @@ def fit(self, X, feature_names):
             knots=self.knots,
             centered=self._mc,
             shared_basis_setups=marginal_shared_setups,
+            metadata=self.metadata,
         )
         for term in marginals:
             term.fit(X, feature_names)
+            if str(getattr(term, "basis_name", "")).lower() == "mrf":
+                self.metadata = dict(term.metadata)
         feature_indices, feature_names_resolved = resolve_tensor_marginal_features(
             marginals
         )
@@ -209,9 +212,7 @@ def point_basis_fn(point):
         self._set_penalty_rescale_factors(
             [
                 float(scale)
-                for scale, keep in zip(
-                    penalty_scales, keep_penalties, strict=True
-                )
+                for scale, keep in zip(penalty_scales, keep_penalties, strict=True)
                 if keep
             ]
         )
diff --git a/nampy/gam/specs/__init__.py b/nampy/gam/specs/__init__.py
index b8d9d702..0b849c2e 100644
--- a/nampy/gam/specs/__init__.py
+++ b/nampy/gam/specs/__init__.py
@@ -12,6 +12,7 @@
     DuchonSplineSmoothSpec,
     FactorSmoothInteractionSpec,
     GaussianProcessSmoothSpec,
+    MarkovRandomFieldSmoothSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
@@ -36,6 +37,7 @@
     "DuchonSplineSmoothSpec",
     "FactorSmoothInteractionSpec",
     "GaussianProcessSmoothSpec",
+    "MarkovRandomFieldSmoothSpec",
     "PSplineSmoothSpec",
     "RandomEffectSmoothSpec",
     "ShapeConstrainedSmoothSpec",
diff --git a/nampy/gam/specs/build.py b/nampy/gam/specs/build.py
index a93add5c..66d980f1 100644
--- a/nampy/gam/specs/build.py
+++ b/nampy/gam/specs/build.py
@@ -477,7 +477,7 @@ def _fs_by_without_factor_feature_base_spec(smooth_spec: FactorSmoothInteraction
 
     if base_bs == "ps":
         kwargs["m"] = None if xt_rest is None else xt_rest.get("m", None)
-    elif base_bs in {"ds", "gp", "tp", "ts"}:
+    elif base_bs in {"ds", "gp", "mrf", "tp", "ts"}:
         if base_bs in {"ds", "gp"}:
             kwargs["m"] = smooth_spec.m
         kwargs["xt"] = xt_rest
diff --git a/nampy/gam/specs/modeling.py b/nampy/gam/specs/modeling.py
index 371ded20..60251958 100644
--- a/nampy/gam/specs/modeling.py
+++ b/nampy/gam/specs/modeling.py
@@ -21,6 +21,19 @@ def make_predictor_specs(model, feature_names, *, knots=None):
 
     for name in feature_names:
         term_knots = knots_for_feature(model, name, knots=knots)
+        model_xt = getattr(model, "xt", None)
+        is_global_mrf_xt = (
+            basis == "mrf"
+            and isinstance(model_xt, dict)
+            and bool({"penalty", "nb", "polys"}.intersection(model_xt))
+        )
+        term_xt = (
+            model_xt.get(str(name))
+            if isinstance(model_xt, dict)
+            and str(name) in model_xt
+            and not is_global_mrf_xt
+            else model_xt
+        )
 
         if basis in {"cr", "cs", "cc"}:
             main_terms.append(
@@ -42,7 +55,7 @@ def make_predictor_specs(model, feature_names, *, knots=None):
                     metadata={},
                 )
             )
-        elif basis in {"bs", "ds", "gp", "ps", "cp"}:
+        elif basis in {"bs", "ds", "gp", "mrf", "ps", "cp"}:
             main_terms.append(
                 TermSpec(
                     kind="smooth",
@@ -57,6 +70,7 @@ def make_predictor_specs(model, feature_names, *, knots=None):
                         fx=False,
                         select=bool(model.select),
                         knots=term_knots,
+                        xt=term_xt,
                     ),
                     smoothing_id=None,
                     label=name,
@@ -105,7 +119,7 @@ def make_predictor_specs(model, feature_names, *, knots=None):
         else:
             raise NotImplementedError(
                 "Automatic main-effect construction currently supports "
-                "{'bs','cr','cs','cc','cp','ds','gp','ps','tp','ts','re'}, "
+                "{'bs','cr','cs','cc','cp','ds','gp','mrf','ps','tp','ts','re'}, "
                 f"got {model.basis!r}."
             )
 
@@ -150,9 +164,7 @@ def prepare_formula_inputs(
         build_result.predictor_specs, build_result.component_lpi, strict=True
     ):
         spec.metadata["lpi"] = tuple(int(v) for v in lpi)
-        spec.metadata["n_linear_predictors"] = int(
-            build_result.n_linear_predictors
-        )
+        spec.metadata["n_linear_predictors"] = int(build_result.n_linear_predictors)
     return (
         parsed,
         build_result.predictor_specs,
diff --git a/nampy/gam/specs/smooth.py b/nampy/gam/specs/smooth.py
index a8164142..d06281bf 100644
--- a/nampy/gam/specs/smooth.py
+++ b/nampy/gam/specs/smooth.py
@@ -76,6 +76,13 @@ class GaussianProcessSmoothSpec(BaseSmoothSpec):
     pc: Any = None
 
 
+@dataclass(frozen=True)
+class MarkovRandomFieldSmoothSpec(BaseSmoothSpec):
+    bs: str = "mrf"
+    xt: Any = None
+    constraint_mode: str = "auto"
+
+
 @dataclass(frozen=True)
 class ShapeConstrainedSmoothSpec(BaseSmoothSpec):
     """SCAM SCOP-spline specification for a named shape basis code."""
@@ -152,6 +159,7 @@ class TensorInteractionSmoothSpec(BaseSmoothSpec):
     DerivativeBSplineSmoothSpec,
     DuchonSplineSmoothSpec,
     GaussianProcessSmoothSpec,
+    MarkovRandomFieldSmoothSpec,
     CubicShrinkageSmoothSpec,
     PSplineSmoothSpec,
     ShapeConstrainedSmoothSpec,
diff --git a/nampy/gam/specs/smooth_build.py b/nampy/gam/specs/smooth_build.py
index 40420e35..06685d58 100644
--- a/nampy/gam/specs/smooth_build.py
+++ b/nampy/gam/specs/smooth_build.py
@@ -16,6 +16,7 @@
     DuchonSplineSmoothSpec,
     FactorSmoothInteractionSpec,
     GaussianProcessSmoothSpec,
+    MarkovRandomFieldSmoothSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
     ShapeConstrainedSmoothSpec,
@@ -145,6 +146,19 @@ def _build_s_gp(opts) -> GaussianProcessSmoothSpec:
     )
 
 
+def _build_s_mrf(opts) -> MarkovRandomFieldSmoothSpec:
+    return MarkovRandomFieldSmoothSpec(
+        special="s",
+        k=opts["k"],
+        fx=opts["fx"],
+        select=opts["select"],
+        sp=opts["sp"],
+        knots=opts["knots"],
+        xt=opts["xt"],
+        constraint_mode=opts["constraint_mode"],
+    )
+
+
 def _build_s_shape(opts) -> ShapeConstrainedSmoothSpec:
     return ShapeConstrainedSmoothSpec(
         special="s",
@@ -236,6 +250,7 @@ def _build_s_sz(opts) -> SumToZeroFactorSmoothSpec:
     "bs": _build_s_bs,
     "ds": _build_s_ds,
     "gp": _build_s_gp,
+    "mrf": _build_s_mrf,
     "tp": _build_s_tp,
     "ts": _build_s_ts,
     "re": _build_s_re,
@@ -553,7 +568,7 @@ def _default_k_for_smooth(kind, basis, features, default_k):
         # mgcv::te()/ti() default k to 5^d per marginal. The current
         # Python tensor surface supports one feature per marginal, so d = 1.
         return [5] * len(features)
-    if str(basis).lower() in {"bs", "ds", "gp", "tp", "ts"}:
+    if str(basis).lower() in {"bs", "ds", "gp", "mrf", "tp", "ts"}:
         # mgcv/R/smooth.r::s() leaves k = -1. The basis constructor then
         # resolves its dimension-dependent default (TP/TS use M + 8/27/100;
         # DS uses M + 10/30/100; GP uses d + 1 + 10/30/100). A flat default here would be wrong in more
diff --git a/nampy/gam/splines/basis/__init__.py b/nampy/gam/splines/basis/__init__.py
index 18878fab..b668b754 100644
--- a/nampy/gam/splines/basis/__init__.py
+++ b/nampy/gam/splines/basis/__init__.py
@@ -1,7 +1,7 @@
 """Low-level basis algebra and invariant helpers for spline primitives."""
 
 from .cr import cr_exact_null_basis_from_knots, cr_spl, cr_spl_predict
-from .natparam import nat_param_type1
+from .natparam import nat_param_type0, nat_param_type1
 from .tp import eta, tp_T
 
 __all__ = [
@@ -10,5 +10,6 @@
     "cr_exact_null_basis_from_knots",
     "eta",
     "tp_T",
+    "nat_param_type0",
     "nat_param_type1",
 ]
diff --git a/nampy/gam/splines/basis/natparam.py b/nampy/gam/splines/basis/natparam.py
index 54bd7054..49eabd8a 100644
--- a/nampy/gam/splines/basis/natparam.py
+++ b/nampy/gam/splines/basis/natparam.py
@@ -237,4 +237,49 @@ def nat_param_type1(X, S, rank=None, tol=None, unit_fnorm=True):
     }
 
 
-__all__ = ["nat_param_type1"]
+def nat_param_type0(X, S, rank=None, tol=None, unit_fnorm=True):
+    """Python implementation of ``mgcv::nat.param(X, S, type=0)``.
+
+    Unlike type 1, the positive natural-parameter penalty eigenvalues are
+    retained rather than normalized to one.  MRF reduced-rank construction
+    uses this exact parameterization before keeping the least penalized
+    columns.
+    """
+    X = np.asarray(X, dtype=np.float64)
+    S = np.asarray(S, dtype=np.float64)
+    tol = np.finfo(float).eps**0.8 if tol is None else float(tol)
+
+    Q, R = _r_linpack_qr(X, tol)
+    if matrix_is_rank_deficient(R):
+        raise ValueError(
+            "Model matrix is not full rank in natural-parameter construction."
+        )
+
+    tmp = _r_triangular_solve(R.T, S.T, lower=True)
+    RSR = _r_triangular_solve(R.T, tmp.T, lower=True)
+    evals, U = _r_symmetric_eigh_descending(RSR)
+
+    if rank is None or rank < 1 or rank > S.shape[0]:
+        max_eval = np.max(evals) if evals.size else 0.0
+        rank = int(np.sum(evals > max_eval * tol))
+    rank = max(0, min(int(rank), S.shape[0]))
+
+    D = np.asarray(evals[:rank], dtype=np.float64).copy()
+    Xn = np.asarray(Q @ U, dtype=np.float64)
+    P = _r_triangular_solve(R, U, lower=False)
+
+    if unit_fnorm:
+        if rank > 0:
+            scale = 1.0 / np.sqrt(np.mean(Xn[:, :rank] ** 2))
+            Xn[:, :rank] *= scale
+            P[:, :rank] *= scale
+            D *= scale**2
+        if rank < Xn.shape[1]:
+            scalef = 1.0 / np.sqrt(np.mean(Xn[:, rank:] ** 2))
+            Xn[:, rank:] *= scalef
+            P[:, rank:] *= scalef
+
+    return {"X": Xn, "D": D, "P": P, "rank": int(rank)}
+
+
+__all__ = ["nat_param_type0", "nat_param_type1"]
diff --git a/pyproject.toml b/pyproject.toml
index 67e09538..8ebdcccf 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -209,6 +209,7 @@ markers = [
     "smooth_cp: tests covering cyclic P-spline smooths",
     "smooth_ds: tests covering Duchon regression spline smooths",
     "smooth_gp: tests covering Gaussian-process smooths",
+    "smooth_mrf: tests covering Markov-random-field smooths",
     "smooth_ps: tests covering P-spline smooths",
     "smooth_tp: tests covering thin plate smooths",
     "smooth_ts: tests covering shrinkage thin plate smooths",
diff --git a/tests/SUBSYSTEM_COVERAGE.md b/tests/SUBSYSTEM_COVERAGE.md
index 6232deb2..1fbfbb88 100644
--- a/tests/SUBSYSTEM_COVERAGE.md
+++ b/tests/SUBSYSTEM_COVERAGE.md
@@ -15,7 +15,7 @@ local development notes rather than duplicated here.
 | Subsystem | Primary owner(s) | Primary tests | Notes |
 | --- | --- | --- | --- |
 | Formula/spec parsing | `nampy/gam/formula/`, `nampy/gam/specs/` | `tests/parity/test_mgcv_formula_parse_parity.py` | Direct formula parity vs `mgcv`. |
-| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py`, `tests/parity/test_mgcv_ds_combinations_parity.py`, `tests/parity/test_mgcv_gp_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`), integrated-derivative B-splines (`bs`), multivariate Duchon splines (`ds`), and five-family stationary/nonstationary Gaussian-process smooths (`gp`) across prediction, selection, linked bases, tensor/factor-smooth combinations, and the upstream tensor-`m` wrong-length warning and zero fallback. |
+| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py`, `tests/parity/test_mgcv_ds_combinations_parity.py`, `tests/parity/test_mgcv_gp_combinations_parity.py`, `tests/parity/test_mgcv_mrf_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`), integrated-derivative B-splines (`bs`), multivariate Duchon splines (`ds`), five-family stationary/nonstationary Gaussian-process smooths (`gp`), and graph/polygon/direct-penalty Markov random fields (`mrf`) across prediction, selection, linked bases, tensor/factor-smooth combinations, and explicit malformed-upstream boundaries. |
 | `pc=` / linked `id=` routing | smooth metadata + linked basis owners | `tests/smooths/test_mgcv_pc_id_parity.py` | Localizes shared-smoothing and point-constraint issues, including `te`/`ti` point constraints. |
 | Design / pre-fit assembly | `nampy/gam/compiler/`, `nampy/gam/fit/penalized_system.py` | `tests/optimization/test_mgcv_gam_setup_assembly_parity.py`, `tests/optimization/test_mgcv_preoptimization_blocks_parity.py`, `tests/optimization/test_mgcv_preoptimization_reparam_parity.py` | Setup, blocks, and reparameterization parity, including one global shared-component block with overlapping linear-predictor indices. |
 | Term wrapping / by-variable / offset routing | predictor wrapping + compiled term owners | `tests/optimization/test_gam_term_wrapping_owner_contracts.py` | Localizes wrapped predictor blocks, offset routing, and general-family block ownership before broader prediction parity. |
@@ -47,7 +47,7 @@ closes the former seven-stage under-tested backlog.
 | Pipeline stage | Combination coverage now owned |
 | --- | --- |
 | 1. Formula parsing and canonical specs | Formula lists/shared components, transformed covariates, supported numeric/factor interactions, factor-by smooths, formula and fit offsets, intercept policies, weights, knots, `min_sp`, `drop_intercept`, and tensor-`m` warning/fallback behavior; the combined interaction recipe is rebuilt on newdata and compared with a committed `mgcv` reference fixture. |
-| 2. Runtime terms and low-level bases | Train/newdata pairing for every supported basis (`bs/cr/cs/cc/cp/ds/gp/ps/tp/ts/re/fs/sz/te/ti`), including univariate, multivariate Duchon/GP, and structured boundary response/SE parity. Row-permutation contracts cover linked and identified-FS terms, PS/TP/TS bases, `te`/`ti`, SZ, and factor-by smooths. Non-unique constructor representations use column-space projectors and penalized response operators instead of arbitrary coefficient orientation. |
+| 2. Runtime terms and low-level bases | Train/newdata pairing for every supported basis (`bs/cr/cs/cc/cp/ds/gp/mrf/ps/tp/ts/re/fs/sz/te/ti`), including univariate, multivariate Duchon/GP, categorical MRF, and structured boundary response/SE parity. Row-permutation contracts cover linked and identified-FS terms, PS/TP/TS bases, `te`/`ti`, SZ, and factor-by smooths. Non-unique constructor representations use column-space projectors and penalized response operators instead of arbitrary coefficient orientation. |
 | 3. Constructed/wrapped terms | Numeric-by plus linked `id=`, factor-by plus linked `id=`, tensor-by, and mixed fixed/free/select penalty ownership, numerical coefficient-map composition, and fitted response/SE parity. |
 | 4. Predictor/model compilation | Multi-predictor layouts with unequal intercept and offset policies, overlapping shared-component coefficient indices, and three-term linked, fixed/free, select, and rank-deficient assembly; supported `gaulss` and `gammals` layouts are also fitted and compared with `mgcv`. |
 | 5. Side conditions and identifiability | Repeated, nested, reverse-formula-order, tensor/main-effect, three-way, both identified near-rank regimes, zero-width, no-intercept, ordered factor-by, linked, general-family, two-predictor, SZ, and exempt random/factor smooth cases. Deletion rank and behavior replace raw pivot-column identity where QR/eigen choices are non-unique. |
diff --git a/tests/TAXONOMY.md b/tests/TAXONOMY.md
index c86ea87d..d0f528c6 100644
--- a/tests/TAXONOMY.md
+++ b/tests/TAXONOMY.md
@@ -12,7 +12,7 @@ The GAM test suite is intentionally overlapping. The goal is fast subset runs an
 - `tests/`: shared helpers, marker inference, taxonomy registry, static reference fixtures, and parity-generation R scripts
 
 ## Taxonomy Axes
-- `smooth_<name>`: `bs`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `ps`, `tp`, `ts`, `te`, `ti`, `fs`, `sz`, `re`
+- `smooth_<name>`: `bs`, `cr`, `cs`, `cc`, `cp`, `ds`, `gp`, `mrf`, `ps`, `tp`, `ts`, `te`, `ti`, `fs`, `sz`, `re`
 - `family_<name>`: `gaussian`, `binomial`, `poisson`, `gamma`, `negbin`, `gaulss`, `gammals`, `general`
 - `method_<name>`: `fixed`, `reml`, `ml`, `laml`, `gcv`, `ubre`
 - `link_<name>`: `identity`, `log`, `inverse`, `logit`, `probit`, `cloglog`, `cauchit`, `sqrt`
diff --git a/tests/_taxonomy_registry.py b/tests/_taxonomy_registry.py
index ae41327c..0fef056d 100644
--- a/tests/_taxonomy_registry.py
+++ b/tests/_taxonomy_registry.py
@@ -21,6 +21,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
     "cp": "smooth_cp",
     "ds": "smooth_ds",
     "gp": "smooth_gp",
+    "mrf": "smooth_mrf",
     "ps": "smooth_ps",
     "tp": "smooth_tp",
     "ts": "smooth_ts",
@@ -87,6 +88,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
     "test_mgcv_smoothcon_parity.py": {"surface_smoothcon"},
     "test_mgcv_raw_constructor_parity.py": {"surface_smoothcon"},
     "test_mgcv_ds_combinations_parity.py": {"surface_snapshot"},
+    "test_mgcv_mrf_combinations_parity.py": {"surface_snapshot"},
     "test_mgcv_score_hist_trace_parity.py": {"surface_trace"},
     "test_mgcv_optimization_lifecycle_parity.py": {"surface_trace"},
     "test_mgcv_linked_id_trace_parity.py": {"surface_trace"},
@@ -122,6 +124,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
     "test_mgcv_known_gaps.py",
     "test_mgcv_output_parity.py",
     "test_mgcv_ds_combinations_parity.py",
+    "test_mgcv_mrf_combinations_parity.py",
     "test_mgcv_score_hist_trace_parity.py",
     "test_mgcv_linked_id_trace_parity.py",
     "test_mgcv_score_gamma_parity.py",
@@ -149,6 +152,10 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
         "tests/smooths/test_mgcv_raw_constructor_parity.py",
         "tests/smooths/test_mgcv_smoothcon_parity.py",
     ),
+    "smooth_mrf": (
+        "tests/parity/test_mgcv_mrf_combinations_parity.py",
+        "tests/smooths/test_mgcv_raw_constructor_parity.py",
+    ),
     "smooth_cs": (
         "tests/smooths/test_mgcv_pc_id_parity.py",
         "tests/smooths/test_mgcv_raw_constructor_parity.py",
diff --git a/tests/mgcv_invariant_policy.py b/tests/mgcv_invariant_policy.py
index 097555bd..0643ee82 100644
--- a/tests/mgcv_invariant_policy.py
+++ b/tests/mgcv_invariant_policy.py
@@ -320,6 +320,17 @@ def _canonicalize_gp_raw_state(state):
     return state
 
 
+def _canonicalize_mrf_raw_state(state):
+    extra = state["extra"]
+    if extra.get("P") is not None:
+        state["X"] = np.asarray(state["X"], dtype=np.float64) @ np.asarray(
+            state["X"], dtype=np.float64
+        ).T
+        extra["P"] = stable_column_space_projector(extra["P"])
+        state["S"] = [penalty_spectrum(S) for S in state["S"]]
+    return state
+
+
 def _canonicalize_cs_raw_state(state):
     state["S"] = [penalty_spectrum(S) for S in state["S"]]
     return state
@@ -378,6 +389,8 @@ def canonicalize_raw_representation_state(state: dict[str, Any]) -> dict[str, An
         return _canonicalize_duchon_raw_state(state)
     if class_name == "gp.smooth":
         return _canonicalize_gp_raw_state(state)
+    if class_name == "mrf.smooth":
+        return _canonicalize_mrf_raw_state(state)
     if class_name == "fs.interaction":
         return _canonicalize_fs_raw_state(state)
     if class_name == "sz.interaction":
diff --git a/tests/mgcv_parity_utils.py b/tests/mgcv_parity_utils.py
index 85a4c65f..dafc3353 100644
--- a/tests/mgcv_parity_utils.py
+++ b/tests/mgcv_parity_utils.py
@@ -1455,6 +1455,13 @@ def _normalize_raw_constructor_knots(knots):
         if value is None:
             out[str(key)] = None
             continue
+        if isinstance(value, pd.Categorical):
+            out[str(key)] = {
+                "__factor__": True,
+                "values": np.asarray(value, dtype=object).ravel().tolist(),
+                "levels": np.asarray(value.categories, dtype=object).ravel().tolist(),
+            }
+            continue
         arr = np.asarray(value, dtype=object).ravel()
         out[str(key)] = arr.tolist()
     return out
@@ -1521,6 +1528,12 @@ def _run_mgcv_raw_constructor(
   kraw <- fromJSON(args[[4]], simplifyVector = FALSE)
   kn <- lapply(kraw, function(v) {
     if (is.null(v)) return(NULL)
+    if (is.list(v) && isTRUE(v$`__factor__`)) {
+      return(factor(
+        unlist(v$values, recursive = TRUE, use.names = FALSE),
+        levels = unlist(v$levels, recursive = TRUE, use.names = FALSE)
+      ))
+    }
     vals <- unlist(v, recursive = TRUE, use.names = FALSE)
     if (is.numeric(vals)) return(unname(as.numeric(vals)))
     if (is.integer(vals)) return(unname(as.integer(vals)))
@@ -1638,6 +1651,12 @@ def _run_mgcv_raw_constructor(
       shift = pack_vector(sm$shift, "numeric"),
       gp_defn = pack_vector(sm$gp.defn, "numeric")
     ),
+    "mrf.smooth" = list(
+      P = pack_matrix(sm$P),
+      knots = pack_vector(sm$knots, "character"),
+      plot_me = isTRUE(sm$plot.me),
+      noterp = isTRUE(sm$noterp)
+    ),
     "random.effect" = list(
       C = pack_constraint(sm$C),
       random = isTRUE(sm$random),
diff --git a/tests/parity/test_mgcv_mrf_combinations_parity.py b/tests/parity/test_mgcv_mrf_combinations_parity.py
new file mode 100644
index 00000000..e63bcb99
--- /dev/null
+++ b/tests/parity/test_mgcv_mrf_combinations_parity.py
@@ -0,0 +1,330 @@
+"""Integrated parity coverage for Markov-random-field smooths (``bs='mrf'``)."""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from nampy.gam import GAM
+from nampy.gam.smooths.categorical.mrf import build_markov_random_field_setup
+from tests.mgcv_parity_utils import (
+    _assert_basic_mgcv_parity,
+    _fit_nampy_model,
+    _fit_nampy_snapshot,
+    _run_mgcv_snapshot,
+)
+
+_PATH_NB = '{"a":["b"],"b":["a","c"],"c":["b","d"],"d":["c","e"],"e":["d"]}'
+_DISCONNECTED_NB = '{"a":["b"],"b":["a"],"c":["d"],"d":["c"],"e":[]}'
+_SPD_PENALTY = "[[2,-1,0,0,0],[-1,3,-1,0,0],[0,-1,3,-1,0],[0,0,-1,3,-1],[0,0,0,-1,2]]"
+
+
+def _mrf_data(seed=921, n=180):
+    rng = np.random.default_rng(seed)
+    region = np.resize(np.asarray(["a", "b", "c", "d", "e"], dtype=object), n)
+    group = np.resize(np.asarray(["u", "v"], dtype=object), n)
+    x = rng.uniform(-1.5, 1.5, size=n)
+    z = 0.8 + rng.uniform(-0.3, 0.5, size=n)
+    effect = {"a": -0.6, "b": -0.2, "c": 0.15, "d": 0.45, "e": 0.75}
+    y = (
+        np.asarray([effect[value] for value in region])
+        + 0.35 * np.sin(1.3 * x)
+        + 0.2 * (group == "v")
+        + rng.normal(scale=0.12, size=n)
+    )
+    return pd.DataFrame({"y": y, "region": region, "group": group, "x": x, "z": z})
+
+
+def _numeric_mrf_data(seed=930, n=150):
+    data = _mrf_data(seed=seed, n=n)
+    data["region"] = np.resize(np.arange(1, 6, dtype=np.float64), n)
+    return data
+
+
+def _assert_snapshot_fit(actual, expected, *, atol=4e-7):
+    for key in ("response", "link"):
+        np.testing.assert_allclose(
+            actual["predictions"][key],
+            expected["predictions"][key],
+            atol=atol,
+            rtol=atol,
+        )
+    np.testing.assert_allclose(
+        actual["fit"]["edf_total"],
+        expected["fit"]["edf_total"],
+        atol=atol,
+        rtol=atol,
+    )
+
+
+def test_mrf_full_rank_reml_matches_mgcv():
+    data = _mrf_data(seed=922)
+    formula = f'y ~ s(region, bs="mrf", xt={{"nb":{_PATH_NB}}})'
+    _assert_basic_mgcv_parity(
+        _fit_nampy_snapshot(data, formula, "gaussian", "REML"),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+        pred_atol=5e-7,
+        pred_rtol=5e-7,
+        sp_log_atol=8e-6,
+        criterion_atol=5e-7,
+    )
+
+
+def test_mrf_numeric_by_select_adds_null_penalty_and_matches_mgcv():
+    data = _mrf_data(seed=923)
+    formula = f'y ~ s(region, by=z, bs="mrf", k=3, xt={{"nb":{_PATH_NB}}})'
+    model = _fit_nampy_model(data, formula, "gaussian", "REML", select=True)
+    assert len(model.smoothing_params) == 2
+    _assert_basic_mgcv_parity(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML", select=True),
+        pred_atol=7e-7,
+        pred_rtol=7e-7,
+        sp_log_atol=5e-5,
+        criterion_atol=6e-7,
+    )
+
+
+def test_mrf_factor_by_fixed_sp_matches_mgcv():
+    data = _mrf_data(seed=924)
+    formula = f'y ~ s(region, by=group, bs="mrf", k=3, xt={{"nb":{_PATH_NB}}}, sp=.7)'
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, formula, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+        atol=8e-7,
+    )
+
+
+def test_mrf_factor_by_linked_id_shares_one_sp_and_matches_mgcv():
+    data = _mrf_data(seed=931)
+    formula = (
+        f'y ~ s(region, by=group, bs="mrf", k=3, '
+        f'xt={{"nb":{_PATH_NB}}}, id="shared_mrf", sp=.7)'
+    )
+    model = _fit_nampy_model(data, formula, "gaussian", "fixed")
+    assert len(model.smoothing_params) == 1
+    _assert_snapshot_fit(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+        atol=8e-7,
+    )
+
+
+def test_mrf_disconnected_select_retains_one_null_penalty_after_centering():
+    data = _mrf_data(seed=925)
+    formula = f'y ~ s(region, bs="mrf", xt={{"nb":{_DISCONNECTED_NB}}})'
+    model = _fit_nampy_model(data, formula, "gaussian", "REML", select=True)
+    assert len(model.smoothing_params) == 2
+    _assert_basic_mgcv_parity(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML", select=True),
+        pred_atol=7e-7,
+        pred_rtol=7e-7,
+        sp_log_atol=1e-5,
+        criterion_atol=6e-7,
+    )
+
+
+def test_mrf_positive_definite_penalty_and_fx_match_mgcv():
+    data = _mrf_data(seed=926)
+    formula = f'y ~ s(region, bs="mrf", xt={{"penalty":{_SPD_PENALTY}}}, sp=.7)'
+    model = _fit_nampy_model(data, formula, "gaussian", "fixed")
+    penalty = model.gam_result_.compiled_model.compiled_penalties[0]
+    assert penalty.rank == penalty.matrix.shape[0]
+    _assert_snapshot_fit(
+        model.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+    )
+
+    fixed_formula = f'y ~ s(region, bs="mrf", xt={{"nb":{_PATH_NB}}}, fx=True)'
+    fixed = _fit_nampy_model(data, fixed_formula, "gaussian", "fixed")
+    assert fixed.gam_result_.compiled_model.compiled_penalties == ()
+    _assert_snapshot_fit(
+        fixed.parity_snapshot(X=data, include_covariances=True),
+        _run_mgcv_snapshot(data, fixed_formula, "gaussian", "REML"),
+    )
+
+
+@pytest.mark.parametrize(
+    "special",
+    ["te", "ti"],
+)
+def test_mrf_tensor_margin_fixed_sp_matches_mgcv(special):
+    data = _mrf_data(seed=927, n=150)
+    formula = (
+        f'y ~ {special}(region, x, bs=["mrf","cr"], k=[3,5], '
+        f'xt=[{{"nb":{_PATH_NB}}},None], sp=[.6,.8])'
+    )
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, formula, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+        atol=1e-6,
+    )
+
+
+@pytest.mark.parametrize("basis", ["fs", "sz"])
+def test_mrf_numeric_factor_smooth_base_fixed_sp_matches_mgcv(basis):
+    data = _numeric_mrf_data(seed=932)
+    numeric_nb = '{"1":["2"],"2":["1","3"],"3":["2","4"],"4":["3","5"],"5":["4"]}'
+    formula = (
+        f'y ~ s(group, region, bs="{basis}", k=5, '
+        f'xt={{"bs":"mrf","nb":{numeric_nb}}}, sp=[.6,.8])'
+    )
+    _assert_snapshot_fit(
+        _fit_nampy_snapshot(data, formula, "gaussian", "fixed"),
+        _run_mgcv_snapshot(data, formula, "gaussian", "REML"),
+        atol=1e-6,
+    )
+
+
+def test_mrf_array_api_and_persistence(tmp_path):
+    data = _mrf_data(seed=928, n=100)
+    nb = {
+        "a": ["b"],
+        "b": ["a", "c"],
+        "c": ["b", "d"],
+        "d": ["c", "e"],
+        "e": ["d"],
+    }
+    X = data[["region"]]
+    model = GAM(
+        family="gaussian",
+        basis="mrf",
+        k=3,
+        xt={"nb": nb},
+        optimize_smoothing=False,
+        smoothing_params=[0.7],
+    ).fit(X=X, y=data["y"].to_numpy(dtype=np.float64))
+    expected = model.predict(X.iloc[:8])
+    path = tmp_path / "mrf.pkl"
+    model.save_model(path)
+    restored = GAM.load_model(path)
+    np.testing.assert_allclose(restored.predict(X.iloc[:8]), expected)
+
+
+def test_mrf_validation_and_linked_boundaries_are_explicit():
+    values = np.resize(np.asarray(["a", "b", "c"], dtype=object), 30)
+    nb = {"a": ["b"], "b": ["a", "c"], "c": ["b"]}
+    with pytest.raises(ValueError, match="must be supplied in xt"):
+        build_markov_random_field_setup(values)
+    with pytest.raises(ValueError, match="dimension set too high"):
+        build_markov_random_field_setup(values, k=4, xt={"nb": nb})
+    with pytest.raises(ValueError, match="k<=2"):
+        build_markov_random_field_setup(values, k=1, xt={"nb": nb})
+    with pytest.raises(ValueError, match="k<=2"):
+        build_markov_random_field_setup(values, k=2, xt={"nb": nb})
+    with pytest.raises(ValueError, match="auto- penalty construction"):
+        build_markov_random_field_setup(
+            values,
+            xt={"nb": {"a": ["b"], "b": [], "c": []}},
+        )
+    with pytest.raises(TypeError, match="uniform representation"):
+        build_markov_random_field_setup(
+            values,
+            xt={"nb": {"a": [2], "b": ["a", "c"], "c": [2]}},
+        )
+    duplicate_named = build_markov_random_field_setup(
+        values,
+        xt={
+            "nb": {
+                "a": ["b", "b", "unknown"],
+                "b": ["a", "c"],
+                "c": ["b"],
+            }
+        },
+    )
+    assert duplicate_named.raw_penalty[0, 0] == 1.0
+    with pytest.raises(ValueError, match="not contained in the knot specification"):
+        build_markov_random_field_setup(
+            np.asarray(["a", "stale"], dtype=object),
+            xt={"nb": nb},
+            factor_levels=["a", "b", "c"],
+        )
+
+    data = _mrf_data(seed=929, n=80)
+    with pytest.raises(NotImplementedError, match="different feature sets"):
+        GAM(
+            formula=(
+                f'y ~ s(region, bs="mrf", xt={{"nb":{_PATH_NB}}}, id="g")'
+                f' + s(group, bs="mrf", xt={{"nb":{{"u":["v"],"v":["u"]}}}}, id="g")'
+            ),
+            optimize_smoothing=False,
+        ).fit(data=data)
+
+    formula = f'y ~ s(region, bs="mrf", xt={{"nb":{_PATH_NB}}}, sp=.7)'
+    model = _fit_nampy_model(data, formula, "gaussian", "fixed")
+    bad = data.iloc[:3].copy()
+    bad.loc[bad.index[0], "region"] = "unknown"
+    with pytest.raises(ValueError, match="unseen levels|unknown regions"):
+        model.predict(bad)
+
+    with pytest.raises(ValueError, match="incorrect number of smoothing parameters"):
+        GAM(
+            formula=(
+                f'y ~ s(region, bs="mrf", xt={{"nb":{_PATH_NB}}}, fx=True, sp=.7)'
+            ),
+            optimize_smoothing=False,
+        ).fit(data=data)
+    with pytest.raises(NotImplementedError, match="pc"):
+        GAM(
+            formula=(f'y ~ s(region, bs="mrf", xt={{"nb":{_PATH_NB}}}, pc="a")'),
+            optimize_smoothing=False,
+        ).fit(data=data)
+
+
+def test_mrf_boolean_regions_and_array_xt_namespace_are_safe():
+    boolean_data = pd.DataFrame(
+        {
+            "y": np.linspace(-0.2, 0.4, 20),
+            "region": np.resize(np.asarray([False, True]), 20),
+        }
+    )
+    boolean_model = GAM(
+        formula=(
+            'y ~ s(region, bs="mrf", '
+            'xt={"nb":{"FALSE":["TRUE"],"TRUE":["FALSE"]}}, sp=.7)'
+        ),
+        optimize_smoothing=False,
+    ).fit(data=boolean_data)
+    prediction = boolean_model.predict(boolean_data.iloc[:4])
+    assert np.isfinite(prediction).all()
+    assert (
+        np.linalg.norm(
+            boolean_model.gam_result_.compiled_model.compiled_terms[0].basis_train
+        )
+        > 0.0
+    )
+
+    array_data = pd.DataFrame(
+        {"nb": np.resize(np.asarray(["a", "b", "c"], dtype=object), 30)}
+    )
+    array_model = GAM(
+        family="gaussian",
+        basis="mrf",
+        k=3,
+        xt={"nb": {"a": ["b"], "b": ["a", "c"], "c": ["b"]}},
+        optimize_smoothing=False,
+        smoothing_params=[0.7],
+    ).fit(X=array_data, y=np.linspace(-0.4, 0.5, len(array_data)))
+    assert np.isfinite(array_model.predict(array_data.iloc[:5])).all()
+
+
+@pytest.mark.parametrize("special", ["te", "ti"])
+def test_mrf_tensor_knots_promote_unobserved_region_for_prediction(special):
+    data = _mrf_data(seed=934, n=80)
+    data = data.loc[data["region"] != "e"].reset_index(drop=True)
+    knots = pd.Categorical(
+        ["a", "b", "c", "d", "e"],
+        categories=["a", "b", "c", "d", "e"],
+    )
+    formula = (
+        f'y ~ {special}(region, x, bs=["mrf","cr"], k=[3,5], '
+        f'xt=[{{"nb":{_PATH_NB}}},None], sp=[.6,.8])'
+    )
+    model = GAM(formula=formula, optimize_smoothing=False).fit(
+        data=data, knots={"region": knots}
+    )
+    newdata = data.iloc[:2].copy()
+    newdata.loc[newdata.index[0], "region"] = "e"
+    assert np.isfinite(model.predict(newdata)).all()
diff --git a/tests/reference_fixtures/mgcv/24692abbf588cc2353db4619413347ed10a0b664d82ecfc26f5dd2befcfb7f8f.json.gz b/tests/reference_fixtures/mgcv/24692abbf588cc2353db4619413347ed10a0b664d82ecfc26f5dd2befcfb7f8f.json.gz
new file mode 100644
index 00000000..2fe1f7aa
Binary files /dev/null and b/tests/reference_fixtures/mgcv/24692abbf588cc2353db4619413347ed10a0b664d82ecfc26f5dd2befcfb7f8f.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/26b50d46b9c090f7fe718a8a41857b727d0344c9027823c795d20976a4c137c7.json.gz b/tests/reference_fixtures/mgcv/26b50d46b9c090f7fe718a8a41857b727d0344c9027823c795d20976a4c137c7.json.gz
new file mode 100644
index 00000000..74806981
Binary files /dev/null and b/tests/reference_fixtures/mgcv/26b50d46b9c090f7fe718a8a41857b727d0344c9027823c795d20976a4c137c7.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/271fd33286df58a5a765b5050dcb4c66f765717beef59e9003e17ad59329c7e7.json.gz b/tests/reference_fixtures/mgcv/271fd33286df58a5a765b5050dcb4c66f765717beef59e9003e17ad59329c7e7.json.gz
new file mode 100644
index 00000000..3f3b76c4
Binary files /dev/null and b/tests/reference_fixtures/mgcv/271fd33286df58a5a765b5050dcb4c66f765717beef59e9003e17ad59329c7e7.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/39c6a5f6d15955a6b3b004a8fcfa117bdffb98e8b512dc5793135a6960af3e81.json.gz b/tests/reference_fixtures/mgcv/39c6a5f6d15955a6b3b004a8fcfa117bdffb98e8b512dc5793135a6960af3e81.json.gz
new file mode 100644
index 00000000..b68cdff0
Binary files /dev/null and b/tests/reference_fixtures/mgcv/39c6a5f6d15955a6b3b004a8fcfa117bdffb98e8b512dc5793135a6960af3e81.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/3ea99ff49e623cb0e62c8e8083c5fcab5e702066187cc9d137ea913031a0136b.json.gz b/tests/reference_fixtures/mgcv/3ea99ff49e623cb0e62c8e8083c5fcab5e702066187cc9d137ea913031a0136b.json.gz
new file mode 100644
index 00000000..25636c36
Binary files /dev/null and b/tests/reference_fixtures/mgcv/3ea99ff49e623cb0e62c8e8083c5fcab5e702066187cc9d137ea913031a0136b.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/4b9d1ba93528d739e9891d79824baf85a97d0ef90b0395ad95d3de1ed8bcda07.json.gz b/tests/reference_fixtures/mgcv/4b9d1ba93528d739e9891d79824baf85a97d0ef90b0395ad95d3de1ed8bcda07.json.gz
new file mode 100644
index 00000000..84e7142f
Binary files /dev/null and b/tests/reference_fixtures/mgcv/4b9d1ba93528d739e9891d79824baf85a97d0ef90b0395ad95d3de1ed8bcda07.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/7176ffa5c1e63d9ef773d854b9c76db2610a413ca4a2405d4ad23b30c7991668.json.gz b/tests/reference_fixtures/mgcv/7176ffa5c1e63d9ef773d854b9c76db2610a413ca4a2405d4ad23b30c7991668.json.gz
new file mode 100644
index 00000000..744eb041
Binary files /dev/null and b/tests/reference_fixtures/mgcv/7176ffa5c1e63d9ef773d854b9c76db2610a413ca4a2405d4ad23b30c7991668.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/7a7e673de690a99c4c73b203538f061a8049cf8c1009a9b302bfb413ae56b475.json.gz b/tests/reference_fixtures/mgcv/7a7e673de690a99c4c73b203538f061a8049cf8c1009a9b302bfb413ae56b475.json.gz
new file mode 100644
index 00000000..4313120b
Binary files /dev/null and b/tests/reference_fixtures/mgcv/7a7e673de690a99c4c73b203538f061a8049cf8c1009a9b302bfb413ae56b475.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/820d08a40e4007add953053315d2fb81f640b34624fcdd4e335dd5dddf8fd255.json.gz b/tests/reference_fixtures/mgcv/820d08a40e4007add953053315d2fb81f640b34624fcdd4e335dd5dddf8fd255.json.gz
new file mode 100644
index 00000000..2d3e3c64
Binary files /dev/null and b/tests/reference_fixtures/mgcv/820d08a40e4007add953053315d2fb81f640b34624fcdd4e335dd5dddf8fd255.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/87597278c5b988e5fe13233855628ae932b291793f67e6b2a61970828e256ca9.json.gz b/tests/reference_fixtures/mgcv/87597278c5b988e5fe13233855628ae932b291793f67e6b2a61970828e256ca9.json.gz
new file mode 100644
index 00000000..e400c096
Binary files /dev/null and b/tests/reference_fixtures/mgcv/87597278c5b988e5fe13233855628ae932b291793f67e6b2a61970828e256ca9.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/9040616bb3018ab240edecbe723efe7a319fb56aac54554ef381189e3dc4c54b.json.gz b/tests/reference_fixtures/mgcv/9040616bb3018ab240edecbe723efe7a319fb56aac54554ef381189e3dc4c54b.json.gz
new file mode 100644
index 00000000..56334c84
Binary files /dev/null and b/tests/reference_fixtures/mgcv/9040616bb3018ab240edecbe723efe7a319fb56aac54554ef381189e3dc4c54b.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/972cff6f163bd961e0dd7a5888f6807084c37b7505f7256addf4e26f3d9984f4.json.gz b/tests/reference_fixtures/mgcv/972cff6f163bd961e0dd7a5888f6807084c37b7505f7256addf4e26f3d9984f4.json.gz
new file mode 100644
index 00000000..826f15e8
Binary files /dev/null and b/tests/reference_fixtures/mgcv/972cff6f163bd961e0dd7a5888f6807084c37b7505f7256addf4e26f3d9984f4.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/9cc8059d8c86e6989a24bf7cd0a02ea5287aa02f7fc49c6adc65b4617278d40c.json.gz b/tests/reference_fixtures/mgcv/9cc8059d8c86e6989a24bf7cd0a02ea5287aa02f7fc49c6adc65b4617278d40c.json.gz
new file mode 100644
index 00000000..1c16ea4a
Binary files /dev/null and b/tests/reference_fixtures/mgcv/9cc8059d8c86e6989a24bf7cd0a02ea5287aa02f7fc49c6adc65b4617278d40c.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/a86349395966d2d9252b7ae12478b767a4cb228001ddb3378fcff149cadaf09d.json.gz b/tests/reference_fixtures/mgcv/a86349395966d2d9252b7ae12478b767a4cb228001ddb3378fcff149cadaf09d.json.gz
new file mode 100644
index 00000000..42236e2e
Binary files /dev/null and b/tests/reference_fixtures/mgcv/a86349395966d2d9252b7ae12478b767a4cb228001ddb3378fcff149cadaf09d.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/a8eedc37e4c14e7bcb1650c7bcb35602beb51f47043321e82e821a268b7a0ac1.json.gz b/tests/reference_fixtures/mgcv/a8eedc37e4c14e7bcb1650c7bcb35602beb51f47043321e82e821a268b7a0ac1.json.gz
new file mode 100644
index 00000000..9c2f9518
Binary files /dev/null and b/tests/reference_fixtures/mgcv/a8eedc37e4c14e7bcb1650c7bcb35602beb51f47043321e82e821a268b7a0ac1.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/c278dfbf9a39db4338cbc402145129cad15f41e94d79d6a726f6b861c4019a7d.json.gz b/tests/reference_fixtures/mgcv/c278dfbf9a39db4338cbc402145129cad15f41e94d79d6a726f6b861c4019a7d.json.gz
new file mode 100644
index 00000000..2723384f
Binary files /dev/null and b/tests/reference_fixtures/mgcv/c278dfbf9a39db4338cbc402145129cad15f41e94d79d6a726f6b861c4019a7d.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/c5299d99e7fa518f78ec67d4e915ee9f6deda98bcfc5385ed2d8c472e4424c61.json.gz b/tests/reference_fixtures/mgcv/c5299d99e7fa518f78ec67d4e915ee9f6deda98bcfc5385ed2d8c472e4424c61.json.gz
new file mode 100644
index 00000000..042c071f
Binary files /dev/null and b/tests/reference_fixtures/mgcv/c5299d99e7fa518f78ec67d4e915ee9f6deda98bcfc5385ed2d8c472e4424c61.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/e09bf05038a2383dfeafe6b62e75027d7ed732cc3b582ca99fef3e71d8676e0d.json.gz b/tests/reference_fixtures/mgcv/e09bf05038a2383dfeafe6b62e75027d7ed732cc3b582ca99fef3e71d8676e0d.json.gz
new file mode 100644
index 00000000..1dc4ed76
Binary files /dev/null and b/tests/reference_fixtures/mgcv/e09bf05038a2383dfeafe6b62e75027d7ed732cc3b582ca99fef3e71d8676e0d.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/e90433678134533e0bd33782cf7a1953a8b62b5a8ec1ff4466c4e0b9aadc2587.json.gz b/tests/reference_fixtures/mgcv/e90433678134533e0bd33782cf7a1953a8b62b5a8ec1ff4466c4e0b9aadc2587.json.gz
new file mode 100644
index 00000000..37be9cb8
Binary files /dev/null and b/tests/reference_fixtures/mgcv/e90433678134533e0bd33782cf7a1953a8b62b5a8ec1ff4466c4e0b9aadc2587.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/f5304b44e26a8362e75f387a53d64e5a8a2cace168c1eaf35c7c0744b1e7d352.json.gz b/tests/reference_fixtures/mgcv/f5304b44e26a8362e75f387a53d64e5a8a2cace168c1eaf35c7c0744b1e7d352.json.gz
new file mode 100644
index 00000000..904d23f4
Binary files /dev/null and b/tests/reference_fixtures/mgcv/f5304b44e26a8362e75f387a53d64e5a8a2cace168c1eaf35c7c0744b1e7d352.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/fa6672ade967e6993708380f389dc9f98724cc28c79f32fa88539e7245ab08e8.json.gz b/tests/reference_fixtures/mgcv/fa6672ade967e6993708380f389dc9f98724cc28c79f32fa88539e7245ab08e8.json.gz
new file mode 100644
index 00000000..69853fca
Binary files /dev/null and b/tests/reference_fixtures/mgcv/fa6672ade967e6993708380f389dc9f98724cc28c79f32fa88539e7245ab08e8.json.gz differ
diff --git a/tests/smooths/test_mgcv_raw_constructor_parity.py b/tests/smooths/test_mgcv_raw_constructor_parity.py
index f16ced60..262de4db 100644
--- a/tests/smooths/test_mgcv_raw_constructor_parity.py
+++ b/tests/smooths/test_mgcv_raw_constructor_parity.py
@@ -20,6 +20,7 @@
     SZSmoothInteractionTerm,
     _block_penalty_for_group,
 )
+from nampy.gam.smooths.categorical.mrf import MarkovRandomFieldTerm
 from nampy.gam.smooths.categorical.re import RandomEffectTerm
 from nampy.gam.smooths.tensor.marginals import build_tensor_product_components
 from nampy.gam.smooths.tensor.te import TensorProductSplineTerm
@@ -123,6 +124,37 @@ def _make_random_effect_numeric_pair_data():
     )
 
 
+def _make_mrf_data(seed=901, n=90, levels=("a", "b", "c", "d", "e")):
+    rng = np.random.default_rng(seed)
+    region = np.resize(np.asarray(levels, dtype=object), n)
+    effects = dict(
+        zip(levels, np.linspace(-0.5, 0.5, len(levels)), strict=True)
+    )
+    y = np.asarray([effects[value] for value in region]) + rng.normal(
+        scale=0.1, size=n
+    )
+    return pd.DataFrame({"y": y, "region": region})
+
+
+def _make_mrf_tensor_data(seed=905, n=90):
+    data = _make_mrf_data(seed=seed, n=n)
+    data["x"] = np.random.default_rng(seed + 1).uniform(-1.0, 1.0, size=n)
+    return data
+
+
+def _make_numeric_mrf_factor_data(seed=906, n=100):
+    rng = np.random.default_rng(seed)
+    region = np.resize(np.arange(1, 6, dtype=np.float64), n)
+    group = np.resize(np.asarray(["u", "v"], dtype=object), n)
+    y = 0.15 * region + 0.2 * (group == "v") + rng.normal(scale=0.1, size=n)
+    return pd.DataFrame({"y": y, "region": region, "group": group})
+
+
+def _mrf_knots_with_unobserved(_data):
+    levels = ["a", "b", "c", "d", "e"]
+    return {"region": pd.Categorical(levels, categories=levels)}
+
+
 def _factory(fn, **kwargs):
     return lambda fn=fn, kwargs=kwargs: fn(**kwargs)
 
@@ -340,6 +372,81 @@ def _build_cubic_case_matrix():
     return cases
 
 
+def _build_mrf_case_matrix():
+    nb_names = '{"a":["b"],"b":["a","c"],"c":["b","d"],"d":["c","e"],"e":["d"]}'
+    nb_indices = '{"a":[2],"b":[1,3],"c":[2,4],"d":[3,5],"e":[4]}'
+    full_rank_penalty = (
+        "[[2,-1,0,0,0],[-1,3,-1,0,0],[0,-1,3,-1,0],"
+        "[0,0,-1,3,-1],[0,0,0,-1,2]]"
+    )
+    disconnected = '{"a":["b"],"b":["a"],"c":["d"],"d":["c"],"e":[]}'
+    numeric_nb = '{"1":["2"],"2":["1","3"],"3":["2","4"],"4":["3","5"],"5":["4"]}'
+    polygons = (
+        '{"a":[[0,0],[1,0],[1,1],[0,1]],'
+        '"b":[[1,0],[2,0],[2,1],[1,1]],'
+        '"c":[[2,0],[3,0],[3,1],[2,1]]}'
+    )
+    return [
+        _case(
+            "mrf_full_named_nb",
+            _factory(_make_mrf_data, seed=901),
+            f'y ~ s(region, bs="mrf", xt={{"nb":{nb_names}}})',
+        ),
+        _case(
+            "mrf_reduced_numeric_nb",
+            _factory(_make_mrf_data, seed=902),
+            f'y ~ s(region, bs="mrf", k=3, xt={{"nb":{nb_indices}}})',
+            atol=2e-8,
+        ),
+        _case(
+            "mrf_supplied_full_rank_penalty",
+            _factory(_make_mrf_data, seed=903),
+            f'y ~ s(region, bs="mrf", xt={{"penalty":{full_rank_penalty}}})',
+        ),
+        _case(
+            "mrf_unobserved_region_knots",
+            _factory(_make_mrf_data, seed=904, levels=("a", "b", "c", "d")),
+            f'y ~ s(region, bs="mrf", k=3, xt={{"nb":{nb_names}}})',
+            atol=2e-8,
+            knots_factory=_mrf_knots_with_unobserved,
+        ),
+        _case(
+            "mrf_disconnected_graph",
+            _factory(_make_mrf_data, seed=905),
+            f'y ~ s(region, bs="mrf", xt={{"nb":{disconnected}}})',
+        ),
+        _case(
+            "mrf_polygon_adjacency",
+            _factory(_make_mrf_data, seed=906, levels=("a", "b", "c")),
+            f'y ~ s(region, bs="mrf", xt={{"polys":{polygons}}})',
+        ),
+        _case(
+            "mrf_te_margin",
+            _factory(_make_mrf_tensor_data, seed=907),
+            f'y ~ te(region, x, bs=["mrf","cr"], k=[3,5], xt=[{{"nb":{nb_names}}},None])',
+            atol=2e-8,
+        ),
+        _case(
+            "mrf_ti_margin",
+            _factory(_make_mrf_tensor_data, seed=908),
+            f'y ~ ti(region, x, bs=["mrf","cr"], k=[3,5], xt=[{{"nb":{nb_names}}},None])',
+            atol=2e-8,
+        ),
+        _case(
+            "mrf_fs_numeric_base",
+            _factory(_make_numeric_mrf_factor_data, seed=909),
+            f'y ~ s(group, region, bs="fs", k=5, xt={{"bs":"mrf","nb":{numeric_nb}}})',
+            atol=2e-8,
+        ),
+        _case(
+            "mrf_sz_numeric_base",
+            _factory(_make_numeric_mrf_factor_data, seed=910),
+            f'y ~ s(group, region, bs="sz", k=5, xt={{"bs":"mrf","nb":{numeric_nb}}})',
+            atol=2e-8,
+        ),
+    ]
+
+
 def _build_ps_case_matrix():
     return [
         _case(
@@ -1024,6 +1131,7 @@ def _build_tensor_case_matrix():
     *_build_tprs_case_matrix(),
     *_build_duchon_case_matrix(),
     *_build_gp_case_matrix(),
+    *_build_mrf_case_matrix(),
     *_build_re_case_matrix(),
     *_build_factor_smooth_case_matrix(),
     *_build_tensor_case_matrix(),
@@ -1262,6 +1370,23 @@ def _serialize_gp_raw(term):
     )
 
 
+def _serialize_mrf_raw(term):
+    setup = term._setup
+    return _common_raw_state(
+        "mrf.smooth",
+        np.asarray(setup.basis_train, dtype=np.float64),
+        [np.asarray(setup.penalty, dtype=np.float64)],
+        rank=int(setup.rank),
+        null_space_dim=int(setup.null_space_dim),
+        extra={
+            "P": None if setup.P is None else np.asarray(setup.P, dtype=np.float64),
+            "knots": list(setup.levels),
+            "plot_me": bool(setup.plot_me),
+            "noterp": True,
+        },
+    )
+
+
 def _serialize_re_raw(term):
     B = np.asarray(term._basis_train, dtype=np.float64)
     q = int(B.shape[1])
@@ -1491,6 +1616,8 @@ def _serialize_term_raw(term, X):
         return _serialize_duchon_raw(term)
     if isinstance(term, GaussianProcessTerm):
         return _serialize_gp_raw(term)
+    if isinstance(term, MarkovRandomFieldTerm):
+        return _serialize_mrf_raw(term)
     if isinstance(term, RandomEffectTerm):
         return _serialize_re_raw(term)
     if isinstance(term, FSmoothInteractionTerm):