OpenTabular · Ananyapam7 · Aug 25, 2026 · Aug 25, 2026
diff --git 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
@@ -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
@@ -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
@@ -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,19 +344,17 @@ 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,)
             )
         )
         if has_explicit_component_lpi
         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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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",
 ]