diff --git a/nampy/gam/diagnostics/concurvity.py b/nampy/gam/diagnostics/concurvity.py
index fd370cb2..63b7b0e8 100644
--- a/nampy/gam/diagnostics/concurvity.py
+++ b/nampy/gam/diagnostics/concurvity.py
@@ -15,6 +15,7 @@
     _term_full_coefficient_indices,
 )
 from ..predict.linear_predictor_matrix import build_lpmatrix
+from ..term_labels import mgcv_term_display_label
 
 
 def _term_indices_for_concurvity(model, n_coef: int):
@@ -37,7 +38,7 @@ def _term_indices_for_concurvity(model, n_coef: int):
         if idx.size == 0:
             continue
         smooth_starts.append(int(np.min(idx)))
-        blocks.append((str(tb.label), idx))
+        blocks.append((mgcv_term_display_label(tb), idx))
 
     if len(blocks) == 0:
         raise ValueError("No smooth or parametric components available for concurvity.")
diff --git a/nampy/gam/diagnostics/derivatives.py b/nampy/gam/diagnostics/derivatives.py
index 4032fbc2..c05ccb4f 100644
--- a/nampy/gam/diagnostics/derivatives.py
+++ b/nampy/gam/diagnostics/derivatives.py
@@ -12,6 +12,7 @@
     _term_blocks_seq,
     _term_full_coefficient_indices,
 )
+from ..term_labels import mgcv_term_display_label
 
 
 @dataclass(frozen=True)
@@ -65,7 +66,7 @@ def smooth_derivative(model, *, X=None, smooth_number: int = 1, deriv: int = 1):
         derivative=values,
         se=np.sqrt(np.maximum(variance, 0.0)),
         derivative_matrix=Xd,
-        term_label=str(term.label),
+        term_label=mgcv_term_display_label(term),
         order=deriv,
     )
 
diff --git a/nampy/gam/diagnostics/k_check.py b/nampy/gam/diagnostics/k_check.py
index 240f6061..fb8cab4c 100644
--- a/nampy/gam/diagnostics/k_check.py
+++ b/nampy/gam/diagnostics/k_check.py
@@ -13,6 +13,7 @@
     _term_blocks_seq,
 )
 from ..predict.linear_predictor_matrix import build_lpmatrix
+from ..term_labels import mgcv_term_display_label
 from .residuals import residuals_gam
 
 
@@ -177,7 +178,7 @@ def k_check(model, subsample: int = 5000, n_rep: int = 400, seed: int | None = N
         edf_by_term = edf_all
     perm_col = 0
     for i, tb in enumerate(term_blocks):
-        label = str(tb.label)
+        label = mgcv_term_display_label(tb)
         X_term = _numeric_feature_block(model, tb, row_idx)
         k_prime = int(tb.coef_slice.stop - tb.coef_slice.start)
         edf = float(edf_by_term[i])
diff --git a/nampy/gam/diagnostics/plots.py b/nampy/gam/diagnostics/plots.py
index e918d05c..d858ed27 100644
--- a/nampy/gam/diagnostics/plots.py
+++ b/nampy/gam/diagnostics/plots.py
@@ -36,7 +36,7 @@
 from ..predict.predictions import _term_has_absorbed_constraint
 from ..predict.terms import _prediction_term_groups
 from ..smooths.categorical import factor_levels_from_metadata
-from ..term_labels import normalize_mgcv_term_label
+from ..term_labels import mgcv_term_display_label
 from .residuals import _prior_weights
 
 __all__ = [
@@ -462,7 +462,8 @@ def prepare_plot_gam_data(
 
     pd_list = []
     for i, tb in enumerate(smooth_blocks):
-        label = _sub_edf(str(tb.label), edf_map.get(id(tb), float("nan")))
+        display_label = mgcv_term_display_label(tb)
+        label = _sub_edf(display_label, edf_map.get(id(tb), float("nan")))
         P = _prepare_smooth(
             model,
             tb,
@@ -535,7 +536,7 @@ def prepare_plot_gam_data(
             P["se"] = False
 
         if partial_resids:
-            normalized = str(normalize_mgcv_term_label(str(tb.label)))
+            normalized = display_label
             group_index = None
             for gi, glabel in enumerate(term_group_labels):
                 if glabel == normalized:
diff --git a/nampy/gam/fit/selection/postfit.py b/nampy/gam/fit/selection/postfit.py
index 2f66a41c..1aaeb2bc 100644
--- a/nampy/gam/fit/selection/postfit.py
+++ b/nampy/gam/fit/selection/postfit.py
@@ -12,7 +12,9 @@
     _n_smoothing_params,
     _penalty_blocks_seq,
     _require_fitted,
+    _term_blocks_seq,
 )
+from ...term_labels import multi_predictor_term_label
 from .criteria.dispatch import criterion_gradient, criterion_hessian, criterion_value
 from .criteria.ml_reml import resolve_ml_reml_scoring_backend
 
@@ -131,6 +133,7 @@ def _gam_vcomp_names(model) -> list[str]:
         return []
 
     names: list[str | None] = [None] * n_sp
+    terms = list(_term_blocks_seq(model))
     for pb in _penalty_blocks_seq(model):
         idx = int(getattr(pb, "smoothing_index", -1))
         if idx < 0 or idx >= n_sp or names[idx] is not None:
@@ -141,7 +144,16 @@ def _gam_vcomp_names(model) -> list[str]:
             label = meta.get("label", None)
         if label is None:
             label = getattr(pb, "label", None)
-        names[idx] = _normalize_vcomp_label(label)
+        normalized = _normalize_vcomp_label(label)
+        term_index = int(getattr(pb, "term_index", -1))
+        if normalized is not None and 0 <= term_index < len(terms):
+            owner = terms[term_index]
+            normalized = multi_predictor_term_label(
+                normalized,
+                predictor_index=int(getattr(owner, "predictor_index", 0)),
+                term_type=str(getattr(owner, "term_type", "")),
+            )
+        names[idx] = normalized
 
     return [name if name is not None else f"sp_{i}" for i, name in enumerate(names)]
 
diff --git a/nampy/gam/inference/anova.py b/nampy/gam/inference/anova.py
index 782c5890..bfa04a26 100644
--- a/nampy/gam/inference/anova.py
+++ b/nampy/gam/inference/anova.py
@@ -34,6 +34,7 @@
     _term_blocks_seq,
     _term_full_coefficient_indices,
 )
+from ..term_labels import mgcv_term_display_label
 from .chi_square_mixtures import psum_chisq
 
 
@@ -45,9 +46,7 @@ def _scale_estimated(model) -> bool:
 
 
 def _formula_term_label(tb) -> str:
-    metadata = dict(getattr(tb, "metadata", {}) or {})
-    formula_term = metadata.get("formula_term", None)
-    return str(getattr(tb, "label", "")) if formula_term is None else str(formula_term)
+    return mgcv_term_display_label(tb, formula_parametric=True)
 
 
 def _parametric_term_groups(model):
@@ -56,7 +55,7 @@ def _parametric_term_groups(model):
         if str(getattr(tb, "term_type", "")) != "parametric":
             continue
         label = _formula_term_label(tb)
-        key = ("parametric", label)
+        key = ("parametric", int(getattr(tb, "predictor_index", 0)), label)
         if groups and groups[-1]["key"] == key:
             groups[-1]["blocks"].append(tb)
             continue
@@ -620,7 +619,7 @@ def _term_table(
         )
         smooth_rows.append(
             {
-                "label": str(tb.label),
+                "label": mgcv_term_display_label(tb),
                 "edf": edf_i,
                 "ref_df": ref_df,
                 "wald_stat": stat_out,
diff --git a/nampy/gam/inference/summary.py b/nampy/gam/inference/summary.py
index 5725eaa9..a9e58dc4 100644
--- a/nampy/gam/inference/summary.py
+++ b/nampy/gam/inference/summary.py
@@ -32,6 +32,7 @@
     _term_blocks_seq,
     _term_full_coefficient_indices,
 )
+from ..term_labels import mgcv_term_display_label
 from .anova import (
     _parametric_term_groups,
     _residual_df,
@@ -116,12 +117,13 @@ def _parametric_coefficient_indices(model) -> tuple[list[int], list[str]]:
     for group in _parametric_term_groups(model):
         for tb in group["blocks"]:
             full_indices = _term_full_coefficient_indices(model, tb)
+            term_label = mgcv_term_display_label(tb)
             for j, full_index in enumerate(full_indices):
                 indices.append(int(full_index))
                 names.append(
-                    str(tb.label)
+                    term_label
                     if full_indices.size == 1
-                    else f"{tb.label}.{j}"
+                    else f"{term_label}.{j}"
                 )
     return indices, names
 
diff --git a/nampy/gam/predict/terms.py b/nampy/gam/predict/terms.py
index 078926fd..1c0b5376 100644
--- a/nampy/gam/predict/terms.py
+++ b/nampy/gam/predict/terms.py
@@ -17,7 +17,7 @@
     _term_blocks_seq,
     _term_full_coefficient_indices,
 )
-from ..term_labels import normalize_mgcv_term_label
+from ..term_labels import multi_predictor_term_label, normalize_mgcv_term_label
 
 
 def _parametric_formula_term(term) -> str | None:
@@ -26,18 +26,6 @@ def _parametric_formula_term(term) -> str | None:
     return None if formula_term is None else str(formula_term)
 
 
-def _multi_predictor_term_label(label: str, *, predictor_index: int, term_type: str):
-    """Apply mgcv's formula-list suffix to later-predictor term labels."""
-    if predictor_index <= 0:
-        return label
-    if term_type == "parametric":
-        return f"{label}.{predictor_index}"
-    open_index = label.find("(")
-    if open_index < 0:
-        return f"{label}.{predictor_index}"
-    return f"{label[:open_index]}.{predictor_index}{label[open_index:]}"
-
-
 def _prediction_term_groups(model):
     """Return ordered mgcv term groups without treating labels as identity."""
     groups = []
@@ -47,7 +35,7 @@ def _prediction_term_groups(model):
         predictor_name = str(getattr(term, "predictor_name", "predictor_0"))
         if term_type == "parametric":
             formula_term = _parametric_formula_term(term)
-            group_label = _multi_predictor_term_label(
+            group_label = multi_predictor_term_label(
                 formula_term or str(getattr(term, "label", "")),
                 predictor_index=predictor_index,
                 term_type=term_type,
@@ -76,7 +64,7 @@ def _prediction_term_groups(model):
             )
             continue
 
-        group_label = _multi_predictor_term_label(
+        group_label = multi_predictor_term_label(
             str(normalize_mgcv_term_label(getattr(term, "label", ""))),
             predictor_index=predictor_index,
             term_type=term_type,
diff --git a/nampy/gam/results/snapshots.py b/nampy/gam/results/snapshots.py
index e74f8e3e..5b276802 100644
--- a/nampy/gam/results/snapshots.py
+++ b/nampy/gam/results/snapshots.py
@@ -37,7 +37,7 @@
     _term_full_coefficient_indices,
 )
 from ..predict.terms import _prediction_term_groups
-from ..term_labels import normalize_mgcv_term_label
+from ..term_labels import mgcv_term_display_label, normalize_mgcv_term_label
 
 
 def _as_pred_or_scalar_array(value):
@@ -430,7 +430,7 @@ def build_parity_snapshot(model, X=None, include_covariances=False):
                 continue
             smooth_blocks.append(tb)
             full_idx = _term_full_coefficient_indices(core, tb)
-            smooth_labels.append(_normalize_reference_term_label(tb.label))
+            smooth_labels.append(mgcv_term_display_label(tb))
             if _cov_bayes(core) is not None:
                 smooth_cov_bayes.append(
                     np.asarray(
diff --git a/nampy/gam/term_labels.py b/nampy/gam/term_labels.py
index 45eb413a..6aa5395d 100644
--- a/nampy/gam/term_labels.py
+++ b/nampy/gam/term_labels.py
@@ -84,4 +84,38 @@ def normalize_mgcv_term_label(label):
     return text
 
 
-__all__ = ["normalize_mgcv_term_label"]
+def multi_predictor_term_label(
+    label: str, *, predictor_index: int, term_type: str
+) -> str:
+    """Apply mgcv's formula-list suffix to later-predictor term labels."""
+    if int(predictor_index) <= 0:
+        return str(label)
+    if str(term_type) == "parametric":
+        return f"{label}.{int(predictor_index)}"
+    open_index = str(label).find("(")
+    if open_index < 0:
+        return f"{label}.{int(predictor_index)}"
+    return f"{label[:open_index]}.{int(predictor_index)}{label[open_index:]}"
+
+
+def mgcv_term_display_label(term, *, formula_parametric: bool = False) -> str:
+    """Return a compiled term's predictor-aware public mgcv label."""
+    term_type = str(getattr(term, "term_type", ""))
+    label = str(getattr(term, "label", ""))
+    if term_type == "parametric" and formula_parametric:
+        metadata = dict(getattr(term, "metadata", {}) or {})
+        label = str(metadata.get("formula_term", label))
+    elif term_type != "parametric":
+        label = str(normalize_mgcv_term_label(label))
+    return multi_predictor_term_label(
+        label,
+        predictor_index=int(getattr(term, "predictor_index", 0)),
+        term_type=term_type,
+    )
+
+
+__all__ = [
+    "mgcv_term_display_label",
+    "multi_predictor_term_label",
+    "normalize_mgcv_term_label",
+]
diff --git a/tests/SUBSYSTEM_COVERAGE.md b/tests/SUBSYSTEM_COVERAGE.md
index aa79c994..717416e9 100644
--- a/tests/SUBSYSTEM_COVERAGE.md
+++ b/tests/SUBSYSTEM_COVERAGE.md
@@ -25,7 +25,7 @@ local development notes rather than duplicated here.
 | ML/REML backend routing | `nampy/gam/fit/selection/criteria/ml_reml.py` | `tests/optimization/test_gam_owner_routing_objective_contracts.py`, `tests/optimization/test_mgcv_gaussian_backend_selection.py` | Exact vs dynamic vs PIRLS vs general-family selection. |
 | Objective wrappers / optimizer wiring | `nampy/gam/fit/selection/optimize/objectives.py`, `.../driver.py` | `tests/optimization/test_gam_owner_routing_objective_contracts.py`, `tests/optimization/test_mgcv_parametric_only_parity.py`, `tests/optimization/test_mgcv_outer_optimization_parity.py`, `tests/optimization/test_mgcv_optimization_lifecycle_parity.py` | Owner contracts first, direct empty-smoothing-vector parity, mgcv trace parity, then lifecycle parity. |
 | Postfit smoothing diagnostics | `nampy/gam/fit/selection/postfit.py` | `tests/optimization/test_gam_postfit_owner_contracts.py`, `tests/optimization/test_mgcv_vcomp_parity.py`, `tests/optimization/test_mgcv_sp_vcov_stage_parity.py` | Endpoint diagnostics, Hessian sourcing, smoothing covariance surfaces, and stage-local `sp.vcov` / unconditional-covariance checkpoints. |
-| General-family fixed-smoothing / postprocess | `nampy/gam/fit/solvers/general_family/fixed_smoothing.py`, `.../newton.py` | `tests/families/test_gam_general_family_owner_contracts.py`, `tests/optimization/test_mgcv_fixed_inner_fit_parity.py`, `tests/optimization/test_mgcv_general_family_preoptimization_parity.py`, `tests/families/test_general_family_mgcv_parity.py`, `tests/parity/test_mgcv_under_tested_supported_combinations.py` | Owner precedence plus mgcv `gam.fit5` parity; reparameterized and original-coordinate singleton/multi-penalty `Sl` blocks are covered through setup, derivatives, roots/totals, and full fits. Structured `re`, `fs`, and linked-`sz` cases include fixed and optimized `fs` behavior and an `fs` block in linear predictor two. |
+| General-family fixed-smoothing / postprocess | `nampy/gam/fit/solvers/general_family/fixed_smoothing.py`, `.../newton.py` | `tests/families/test_gam_general_family_owner_contracts.py`, `tests/optimization/test_mgcv_fixed_inner_fit_parity.py`, `tests/optimization/test_mgcv_general_family_preoptimization_parity.py`, `tests/families/test_general_family_mgcv_parity.py`, `tests/parity/test_mgcv_snapshot_extended_matrix.py`, `tests/parity/test_mgcv_under_tested_supported_combinations.py` | Owner precedence plus mgcv `gam.fit5` parity; reparameterized and original-coordinate singleton/multi-penalty `Sl` blocks are covered through setup, derivatives, roots/totals, and full fits. A four-smooth/two-predictor gaulss case covers wrapped blocks, `sp.vcov`, inference, and diagnostics. Structured `re`, `fs`, and linked-`sz` cases include fixed and optimized `fs` behavior and an `fs` block in linear predictor two. |
 | Diagnostics owners | `nampy/gam/diagnostics/residuals.py`, `concurvity.py`, `summary.py`, `plots.py` | `tests/diagnostics/test_gam_diagnostics_owner_contracts.py`, `tests/diagnostics/test_gam_plot_and_public_api_contracts.py`, `tests/parity/test_mgcv_secondary_diagnostics_parity.py` | Owner-level residual/summary/plot contracts plus direct secondary-diagnostics parity. |
 | Prediction / inference / diagnostics | `nampy/gam/predict/`, `nampy/gam/inference/`, `nampy/gam/diagnostics/` | `tests/parity/test_mgcv_output_parity.py`, `tests/parity/test_mgcv_prediction_arguments_parity.py`, `tests/parity/test_mgcv_prediction_inference_diagnostics_parity.py`, `tests/parity/test_mgcv_general_family_lpmatrix_stage_parity.py`, `tests/parity/test_mgcv_general_family_prediction_stage_parity.py`, `tests/parity/test_mgcv_inference_stage_parity.py`, `tests/diagnostics/test_mgcv_general_family_secondary_diagnostics_parity.py` | Public-surface parity plus direct `block.size`, `newdata.guaranteed`, `na.action`, `unconditional`, `iterms.type`, and stage-local general-family checkpoints. |
 | Parity snapshot / trace tooling | `nampy/gam/parity/` | `tests/parity/test_gam_parity_owner_contracts.py`, `tests/parity/test_gam_results_api_stage_owner_contracts.py`, `tests/optimization/test_mgcv_score_hist_trace_parity.py`, `tests/optimization/test_mgcv_outer_optimization_parity.py`, `tests/optimization/test_mgcv_optimization_lifecycle_parity.py`, `tests/optimization/test_mgcv_inner_trace_parity.py`, `tests/optimization/test_mgcv_joint_branch_trace_parity.py` | Localizes serialization, criterion-view logic, outer-object trace schemas, lifecycle branch parity, and inner/joint trace branches. |
diff --git a/tests/families/test_general_family_mgcv_parity.py b/tests/families/test_general_family_mgcv_parity.py
index cf94bb64..eb7e4e89 100644
--- a/tests/families/test_general_family_mgcv_parity.py
+++ b/tests/families/test_general_family_mgcv_parity.py
@@ -53,6 +53,22 @@ def _gaulss_tensor_data(seed=22, n=160):
     return pd.DataFrame({"y": y, "x0": x0, "x1": x1})
 
 
+GENERAL_MULTISMOOTH_FORMULA = [
+    'y ~ s(x0, bs="cr", k=8) + s(x1, bs="cr", k=8)',
+    '~ s(x0, bs="cr", k=7) + s(x1, bs="cr", k=7)',
+]
+
+
+def _gaulss_multismooth_data(seed=1501, n=220):
+    rng = np.random.default_rng(seed)
+    x0 = rng.uniform(-1.5, 1.5, size=n)
+    x1 = rng.uniform(-1.5, 1.5, size=n)
+    mu = 0.3 + 1.8 * np.sin(2.5 * x0) + 1.1 * np.cos(3.0 * x1)
+    sigma = np.exp(-1.0 + 0.45 * np.sin(2.0 * x0) - 0.35 * np.cos(2.5 * x1))
+    y = rng.normal(mu, sigma, size=n)
+    return pd.DataFrame({"y": y, "x0": x0, "x1": x1})
+
+
 def _gammals_data(n=100, seed=2):
     rng = np.random.default_rng(seed)
     x = rng.uniform(-1.0, 1.0, n)
diff --git a/tests/optimization/test_mgcv_general_family_preoptimization_parity.py b/tests/optimization/test_mgcv_general_family_preoptimization_parity.py
index 708c1677..ac790b9c 100644
--- a/tests/optimization/test_mgcv_general_family_preoptimization_parity.py
+++ b/tests/optimization/test_mgcv_general_family_preoptimization_parity.py
@@ -20,10 +20,12 @@
 from tests._paths import PARITY_DIR, REPO_ROOT
 from tests.families.test_general_family_mgcv_parity import (
     GAULSS_FORMULA,
+    GENERAL_MULTISMOOTH_FORMULA,
     _gammals_by_data,
     _gammals_data,
     _gaulss_by_data,
     _gaulss_data,
+    _gaulss_multismooth_data,
 )
 from tests.mgcv_parity_utils import _family_specs, _fit_nampy_model_fixed_sp
 from tests.reference_fixtures import (
@@ -404,6 +406,15 @@ def _assert_general_fit5_setup_parity(
 
 GENERAL_PREOPT_CASES = [
     ("gaulss_cr", "gaulss", GAULSS_FORMULA, _gaulss_data, "ML", False, True),
+    (
+        "gaulss_multi_smooth_both_predictors",
+        "gaulss",
+        GENERAL_MULTISMOOTH_FORMULA,
+        _gaulss_multismooth_data,
+        "ML",
+        False,
+        True,
+    ),
     (
         "gaulss_fs",
         "gaulss",
diff --git a/tests/parity/test_mgcv_general_family_prediction_filters_parity.py b/tests/parity/test_mgcv_general_family_prediction_filters_parity.py
index 09f28466..75d84120 100644
--- a/tests/parity/test_mgcv_general_family_prediction_filters_parity.py
+++ b/tests/parity/test_mgcv_general_family_prediction_filters_parity.py
@@ -3,6 +3,7 @@
 import numpy as np
 import pytest
 
+from nampy.gam.inference.summary import summary_gam
 from tests.families.test_general_family_mgcv_parity import (
     _gaulss_by_data,
     _general_newdata,
@@ -17,6 +18,19 @@
 ]
 
 
+def test_general_family_inference_uses_predictor_aware_term_labels():
+    """Formula-list parametric and smooth rows use mgcv's later-LP suffix."""
+    data = _gaulss_by_data(seed=270, n=120)
+    gam = _fit_nampy_model(data, _FORMULA, "gaulss", "fixed")
+
+    summary = summary_gam(gam)
+    assert list(summary.pterms_table["label"]) == ["x", "z.1"]
+    assert list(summary.s_table["label"]) == ["s(z)", "s.1(x)"]
+    anova = gam.anova()
+    assert list(anova.parametric_table["label"]) == ["x", "z.1"]
+    assert list(anova.smooth_table["label"]) == ["s(z)", "s.1(x)"]
+
+
 def test_general_family_terms_filter_values_labels_and_se_match_mgcv():
     data = _gaulss_by_data(seed=271, n=120)
     newdata = _general_newdata(data, n=19)
diff --git a/tests/parity/test_mgcv_snapshot_extended_matrix.py b/tests/parity/test_mgcv_snapshot_extended_matrix.py
index b762fcf9..1e0604f3 100644
--- a/tests/parity/test_mgcv_snapshot_extended_matrix.py
+++ b/tests/parity/test_mgcv_snapshot_extended_matrix.py
@@ -6,9 +6,14 @@
 import pandas as pd
 import pytest
 
+from nampy.gam.inference.summary import summary_gam
 from tests._mgcv_snapshot_parity_shared import (
     TestAdditionalScenarioParity as _SharedTestAdditionalScenarioParity,
 )
+from tests.families.test_general_family_mgcv_parity import (
+    GENERAL_MULTISMOOTH_FORMULA,
+    _gaulss_multismooth_data,
+)
 from tests.mgcv_parity_utils import (
     _assert_basic_mgcv_parity,
     _assert_exact_mgcv_snapshot_parity,
@@ -804,6 +809,118 @@ def test_general_family_gaulss_and_gammals_tensor_multi_smooth_predictions_match
             )
 
 
+def test_general_family_multismooth_fit_inference_and_diagnostics_match_mgcv():
+    """Close the wrapped-block surface across both gaulss predictors."""
+    data = _gaulss_multismooth_data()
+    formula = GENERAL_MULTISMOOTH_FORMULA
+    expected = _run_mgcv_snapshot(data, formula, "gaulss", "ML")
+    model = _fit_nampy_model(data, formula, "gaulss", "ML")
+    core = model.gam_result_.fit_summary.core
+    expected_fit = expected["fit"]
+    diagnostics = expected["parity"]["diagnostics"]
+
+    np.testing.assert_allclose(
+        model.smoothing_params,
+        expected_fit["smoothing_params"],
+        atol=2e-5,
+        rtol=2e-6,
+    )
+    np.testing.assert_allclose(
+        core.coef_full, expected_fit["coef_full"], atol=5e-7, rtol=5e-7
+    )
+    for name in ("cov_bayes", "cov_freq", "cov_unconditional"):
+        np.testing.assert_allclose(
+            np.asarray(getattr(core, name), dtype=np.float64),
+            np.asarray(expected_fit[name], dtype=np.float64),
+            atol=5e-7,
+            rtol=5e-7,
+        )
+    assert model.smoothing_score_ == pytest.approx(
+        float(expected_fit["criterion_value"]), abs=1e-6
+    )
+    assert model.loglik() == pytest.approx(float(expected_fit["loglik"]), abs=1e-6)
+    assert model.aic() == pytest.approx(float(expected_fit["aic"]), abs=2e-6)
+    assert model.bic() == pytest.approx(float(diagnostics["bic"]), abs=3e-6)
+
+    np.testing.assert_allclose(
+        model.sp_vcov(edge_correct=False),
+        diagnostics["sp_vcov"],
+        atol=2e-6,
+        rtol=2e-6,
+    )
+    np.testing.assert_allclose(
+        model.one_se_rule(),
+        diagnostics["one_se_rule"],
+        atol=5e-5,
+        rtol=2e-6,
+    )
+    vcomp = model.gam_vcomp(rescale=False)
+    assert list(vcomp["names"]) == list(diagnostics["gam_vcomp_names"])
+    np.testing.assert_allclose(
+        vcomp["vc"], diagnostics["gam_vcomp"], atol=2e-6, rtol=2e-6
+    )
+
+    summary = summary_gam(model)
+    expected_summary = diagnostics["summary"]
+    assert list(summary.p_table.index) == list(expected_summary["p_table"]["labels"])
+    np.testing.assert_allclose(
+        summary.p_table.to_numpy(dtype=np.float64),
+        expected_summary["p_table"]["values"],
+        atol=2e-6,
+        rtol=2e-6,
+    )
+    expected_smooth = diagnostics["anova_smooth"]
+    assert list(summary.s_table["label"]) == list(expected_smooth["labels"])
+    np.testing.assert_allclose(
+        summary.s_table[["edf", "ref_df", "wald_stat", "p_value"]].to_numpy(
+            dtype=np.float64
+        ),
+        expected_smooth["values"],
+        atol=2e-5,
+        rtol=2e-6,
+    )
+    anova = model.anova()
+    assert list(anova.smooth_table["label"]) == list(expected_smooth["labels"])
+
+    expected_residuals = diagnostics["residuals"]
+    for residual_type in ("response", "pearson", "deviance"):
+        np.testing.assert_allclose(
+            model.residuals(type=residual_type),
+            expected_residuals[residual_type],
+            atol=2e-6,
+            rtol=2e-6,
+        )
+
+    full = model.concurvity(full=True)
+    assert list(full["labels"]) == list(diagnostics["concurvity_labels"])
+    # Two predictors contain smooths of the same covariates, making the
+    # function-space estimate nearly singular; the observed/worst rows remain
+    # tight while the estimate row is stable to about 0.7% across QR stacks.
+    np.testing.assert_allclose(
+        full["values"], diagnostics["concurvity_full"], atol=7e-3, rtol=0.0
+    )
+    pairwise = model.concurvity(full=False)
+    expected_pairwise = diagnostics["concurvity_pairwise"]
+    assert list(pairwise["labels"]) == list(expected_pairwise["labels"])
+    for name, values in pairwise["values"].items():
+        np.testing.assert_allclose(
+            values, expected_pairwise[name], atol=2e-7, rtol=2e-7
+        )
+
+    k_table = model.k_check(subsample=120, n_rep=8, seed=0)
+    expected_k = diagnostics["k_check"]
+    assert list(k_table.index) == list(expected_k["labels"])
+    actual_k = k_table[["k_prime", "edf", "k_index", "p_value"]].to_numpy(
+        dtype=np.float64
+    )
+    expected_k_values = np.asarray(expected_k["values"], dtype=np.float64)
+    np.testing.assert_allclose(actual_k[:, :2], expected_k_values[:, :2], atol=2e-6)
+    np.testing.assert_allclose(
+        actual_k[:, 2], expected_k_values[:, 2], atol=0.4, rtol=0.5
+    )
+    assert np.all((actual_k[:, 3] >= 0.0) & (actual_k[:, 3] <= 1.0))
+
+
 def _snapshot_matrix_assert(actual, expected, *, atol=1e-5):
     for key in ("response", "link"):
         np.testing.assert_allclose(
diff --git a/tests/reference_fixtures/mgcv/dcd66fd99ca234c55a4089e59090e447fc3d0f2616f08bad7b9206f457cd3c93.json.gz b/tests/reference_fixtures/mgcv/dcd66fd99ca234c55a4089e59090e447fc3d0f2616f08bad7b9206f457cd3c93.json.gz
new file mode 100644
index 00000000..b99b884a
Binary files /dev/null and b/tests/reference_fixtures/mgcv/dcd66fd99ca234c55a4089e59090e447fc3d0f2616f08bad7b9206f457cd3c93.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/f35b0e352ce1af0906a9764c8b12ad3e79317d5dbb9924c59ea49662b01ac8a6.json.gz b/tests/reference_fixtures/mgcv/f35b0e352ce1af0906a9764c8b12ad3e79317d5dbb9924c59ea49662b01ac8a6.json.gz
new file mode 100644
index 00000000..3f09c585
Binary files /dev/null and b/tests/reference_fixtures/mgcv/f35b0e352ce1af0906a9764c8b12ad3e79317d5dbb9924c59ea49662b01ac8a6.json.gz differ