diff --git a/src/hssm/addm/addm.py b/src/hssm/addm/addm.py
index 4392876a7..52888b390 100644
--- a/src/hssm/addm/addm.py
+++ b/src/hssm/addm/addm.py
@@ -152,6 +152,7 @@ def _make_model_distribution(self) -> type[pm.Distribution]:
             lapse=self.lapse,
             extra_fields=extra_fields_data,
             params_is_trialwise=params_is_trialwise,
+            expected_obs_dim=self._obs_dim,
         )
         # Expose the observed fixations to the RV's generative path so
         # posterior-predictive draws condition on them (see _push_rv_extra_fields).
diff --git a/src/hssm/base.py b/src/hssm/base.py
index 1291c56f2..669a20169 100644
--- a/src/hssm/base.py
+++ b/src/hssm/base.py
@@ -311,14 +311,13 @@ def __init__(
         # region ===== Store the pre-built config =====
         self.model_config: BaseModelConfig = copy(model_config)
         self.model_config.response_domains = deepcopy(model_config.response_domains)
+        configured_response = self.model_config.response
+        assert configured_response is not None
+        self._obs_dim = len(configured_response)
         # endregion
 
         # region ===== Set up shortcuts so old code will work ======
-        self.response: list[str] = (  # type: ignore[assignment]
-            list(self.model_config.response)
-            if self.model_config.response is not None
-            else []
-        )
+        self.response: list[str] = list(configured_response)  # type: ignore[assignment]
         self.list_params = (
             list(self.model_config.list_params)
             if self.model_config.list_params is not None
@@ -391,6 +390,7 @@ def __init__(
             self.list_params,
             self.link,
             self._parent,
+            obs_dim=self._obs_dim,
         )
 
         # Targeted checks against the user's prior dict:
@@ -1277,18 +1277,7 @@ def sample_do(
         )
         do_dt = pm.sample_prior_predictive(model=do_model, draws=draws, **kwargs)
 
-        # clean up `rt,response_mean` to `v`
-        do_dt = self._drop_parent_str_from_datatree(dt=do_dt)
-
-        # rename otherwise inconsistent dims and coords
-        if "rt,response_extra_dim_0" in do_dt["prior_predictive"].dims:
-            do_dt["prior_predictive"] = do_dt["prior_predictive"].ds.rename_dims(
-                {"rt,response_extra_dim_0": "rt,response_dim"}
-            )
-        if "rt,response_extra_dim_0" in do_dt["prior_predictive"].coords:
-            do_dt["prior_predictive"] = do_dt["prior_predictive"].ds.rename_vars(
-                {"rt,response_extra_dim_0": "rt,response_dim"}
-            )
+        do_dt = self._clean_predictive_datatree(dt=do_dt)
 
         if return_model:
             return do_dt, do_model
@@ -1346,18 +1335,7 @@ def sample_prior_predictive(
                     continue
                 self._inference_obj[group_name] = prior_predictive[group_name]
 
-        # clean up `rt,response_mean` to `v`
-        dt = self._drop_parent_str_from_datatree(dt=self._inference_obj)
-
-        # rename otherwise inconsistent dims and coords
-        if "rt,response_extra_dim_0" in dt["prior_predictive"].dims:
-            dt["prior_predictive"] = dt["prior_predictive"].ds.rename_dims(
-                {"rt,response_extra_dim_0": "rt,response_dim"}
-            )
-        if "rt,response_extra_dim_0" in dt["prior_predictive"].coords:
-            dt["prior_predictive"] = dt["prior_predictive"].ds.rename_vars(
-                name_dict={"rt,response_extra_dim_0": "rt,response_dim"}
-            )
+        dt = self._clean_predictive_datatree(dt=self._inference_obj)
 
         # Update self._inference_obj to match the cleaned datatree
         self._inference_obj = dt
@@ -2022,8 +2000,8 @@ def _get_deterministic_var_names(self, dt) -> list[str]:
 
         return var_names
 
-    def _drop_parent_str_from_datatree(self, dt: DataTree | None) -> DataTree:
-        """Drop the parent_str variable from a DataTree object.
+    def _clean_predictive_datatree(self, dt: DataTree | None) -> DataTree:
+        """Normalize generated response names and dimensions.
 
         Parameters
         ----------
@@ -2032,20 +2010,32 @@ def _drop_parent_str_from_datatree(self, dt: DataTree | None) -> DataTree:
 
         Returns
         -------
-        xr.Dataset
+        DataTree
             The modified DataTree object.
         """
         if dt is None:
             raise ValueError("Please provide a DataTree (traces) object.")
-        else:
-            for group in dt.groups:
-                if group == "/":
-                    continue
-                if ("rt,response_mean" in dt[group].data_vars) and (
-                    self._parent not in dt[group].data_vars
-                ):
-                    dt[group] = dt[group].ds.rename({"rt,response_mean": self._parent})
-            return dt
+
+        response_mean = f"{self.response_str}_mean"
+        raw_response_dim = f"{self.response_str}_extra_dim_0"
+        response_dim = f"{self.response_str}_dim"
+        for group in dt.groups:
+            if group == "/":
+                continue
+            dataset = dt[group].ds
+            rename = {}
+            if response_mean in dataset.data_vars and self._parent not in dataset:
+                rename[response_mean] = self._parent
+            if group.rsplit("/", maxsplit=1)[-1] in {
+                "posterior_predictive",
+                "prior_predictive",
+            } and (
+                raw_response_dim in dataset.dims or raw_response_dim in dataset.coords
+            ):
+                rename[raw_response_dim] = response_dim
+            if rename:
+                dt[group] = dataset.rename(rename)
+        return dt
 
     def _postprocess_initvals_deterministic(
         self, initval_settings: dict = INITVAL_SETTINGS
diff --git a/src/hssm/distribution_utils/dist.py b/src/hssm/distribution_utils/dist.py
index 738873020..f62e4d41b 100644
--- a/src/hssm/distribution_utils/dist.py
+++ b/src/hssm/distribution_utils/dist.py
@@ -7,6 +7,7 @@
 
 import logging
 from collections.abc import Callable
+from functools import partial
 from os import PathLike
 from typing import Any, Literal, Protocol, cast, get_args
 
@@ -195,6 +196,7 @@ def make_hssm_rv(
     list_params: list[str],
     lapse: bmb.Prior | float | None = None,
     is_choice_only: bool = False,
+    expected_obs_dim: int | None = None,
 ) -> type[RandomVariable]:
     """Build a RandomVariable Op according to the list of parameters.
 
@@ -208,6 +210,9 @@ def make_hssm_rv(
         A bmb.Prior object representing the lapse distribution.
     is_choice_only : bool
         Whether the model is a choice-only model.
+    expected_obs_dim : int, optional
+        Configured number of physical response columns. When supplied, it must
+        agree with the simulator's effective observation width.
 
     Returns
     -------
@@ -215,7 +220,18 @@ def make_hssm_rv(
         A class of RandomVariable that are to be used in a `pm.Distribution`.
     """
     simulator_fun_internal = get_simulator_fun_internal(simulator_fun)
-    model_name, choices, obs_dim_int = validate_simulator_fun(simulator_fun_internal)
+    model_name, choices, declared_obs_dim = validate_simulator_fun(
+        simulator_fun_internal
+    )
+    effective_obs_dim = 1 if is_choice_only else declared_obs_dim
+    comparison_obs_dim = (
+        declared_obs_dim if callable(simulator_fun) else effective_obs_dim
+    )
+    if expected_obs_dim is not None and comparison_obs_dim != expected_obs_dim:
+        raise ValueError(
+            f"The simulator observation width {comparison_obs_dim} does not match "
+            f"the configured response width {expected_obs_dim}."
+        )
 
     if lapse is not None and list_params[-1] != "p_outlier":
         list_params.append("p_outlier")
@@ -249,7 +265,7 @@ class HSSMRV(RandomVariable):
 
         # Override the output from ssm_simulator based on whether the model is
         # choice-only.
-        output = "()" if is_choice_only else f"({obs_dim_int})"
+        output = "()" if is_choice_only else f"({effective_obs_dim})"
         signature = f"{','.join(['()'] * len(list_params))}->{output}"
 
         dtype = "floatX"
@@ -319,7 +335,7 @@ def rng_fn(
                 size,
                 rng,
                 simulator_fun_internal,
-                obs_dim_int,
+                effective_obs_dim,
                 *args,
                 **kwargs,
             )
@@ -440,6 +456,7 @@ def make_distribution(
     fixed_vector_params: dict[str, np.ndarray] | None = None,
     params_is_trialwise: list[bool] | None = None,
     is_choice_only: bool = False,
+    expected_obs_dim: int | None = None,
 ) -> type[pm.Distribution]:
     """Make a `pymc.Distribution`.
 
@@ -486,6 +503,9 @@ def make_distribution(
         When ``None``, no graph-level broadcasting is applied.
     is_choice_only : optional
         Whether the model is a choice-only model.
+    expected_obs_dim : int, optional
+        Configured number of physical response columns. This is checked against
+        generated RandomVariable metadata.
 
     Returns
     -------
@@ -502,6 +522,7 @@ def make_distribution(
             list_params=list_params,
             lapse=lapse,
             is_choice_only=is_choice_only,
+            expected_obs_dim=expected_obs_dim,
         )
         rv_instance = random_variable()
     elif isinstance(rv, str):
@@ -510,6 +531,7 @@ def make_distribution(
             list_params=list_params,
             lapse=lapse,
             is_choice_only=is_choice_only,
+            expected_obs_dim=expected_obs_dim,
         )
         rv_instance = random_variable()
     else:
@@ -718,6 +740,7 @@ def make_family(
     parent: str = "v",
     likelihood_name: str = "SSM Likelihood",
     family_name="SSM Family",
+    obs_dim: int = 2,
 ) -> bmb.Family:
     """Build a family in bambi.
 
@@ -735,6 +758,8 @@ def make_family(
         the name of the likelihood function. Defaults to "SSM Likelihood".
     family_name
         the name of the family. Defaults to "SSM Family".
+    obs_dim
+        Number of physical response columns in posterior-predictive draws.
 
     Returns
     -------
@@ -745,7 +770,12 @@ def make_family(
         likelihood_name, parent=parent, params=list_params, dist=dist
     )
 
-    family = SSMFamily(family_name, likelihood=likelihood, link=link)
+    family = SSMFamily(
+        family_name,
+        likelihood=likelihood,
+        link=link,
+        obs_dim=obs_dim,
+    )
 
     return family
 
@@ -753,9 +783,10 @@ def make_family(
 class SSMFamily(bmb.Family):
     """Extends bmb.Family to get around the dimensionality mismatch."""
 
-    def create_extra_pps_coord(self):
-        """Create an extra dimension."""
-        return np.arange(2)
+    def __init__(self, name, likelihood, link, obs_dim: int = 2):
+        super().__init__(name, likelihood, link)
+        if obs_dim > 1:
+            self.create_extra_pps_coord = partial(np.arange, obs_dim)
 
 
 def make_likelihood_callable(
diff --git a/src/hssm/hssm.py b/src/hssm/hssm.py
index aad7dcf35..4dae1d7eb 100644
--- a/src/hssm/hssm.py
+++ b/src/hssm/hssm.py
@@ -435,6 +435,7 @@ def dummy_simulator_func(*args, **kwargs):
                 list_params=self.list_params or [],
                 lapse=self.lapse,
                 is_choice_only=True,
+                expected_obs_dim=self._obs_dim,
             )
 
         self.data = typing_cast("pd.DataFrame", _rearrange_data(self.data))
@@ -469,4 +470,5 @@ def dummy_simulator_func(*args, **kwargs):
             params_is_trialwise=params_is_trialwise_base,
             # TODO: add to HSSMBase
             is_choice_only=self.is_choice_only,
+            expected_obs_dim=self._obs_dim,
         )
diff --git a/src/hssm/rl/rlssm.py b/src/hssm/rl/rlssm.py
index 34fdb17e6..c6b33d7fe 100644
--- a/src/hssm/rl/rlssm.py
+++ b/src/hssm/rl/rlssm.py
@@ -307,6 +307,7 @@ def _make_model_distribution(self) -> type[pm.Distribution]:
             extra_fields=extra_fields_data,
             params_is_trialwise=params_is_trialwise,
             is_choice_only=self.model_config.is_choice_only,
+            expected_obs_dim=self._obs_dim,
         )
 
 
diff --git a/src/hssm/utils.py b/src/hssm/utils.py
index 8d4e1c15e..2a87ea4d3 100644
--- a/src/hssm/utils.py
+++ b/src/hssm/utils.py
@@ -576,19 +576,35 @@ def predictive_dt_to_dataframe(
         The DataTree object to convert.
     predictive_group : Literal["posterior_predictive", "prior_predictive"]
         The predictive group to convert.
+    response_str
+        Comma-separated physical response column names in configured order.
+    response_dim
+        Name of the vector response coordinate. Scalar responses have no such
+        coordinate.
 
     Returns
     -------
     pd.DataFrame:
         A dataframe with the predictive samples.
     """
+    predictive = dt[predictive_group].ds[response_str]
     df = dt[predictive_group].ds.to_dataframe().reset_index(drop=False)
+    response_names = response_str.split(",")
+    if response_dim not in df.columns:
+        if any(dim not in {"chain", "draw", "__obs__"} for dim in predictive.dims):
+            raise ValueError(
+                f"Predictive response coordinate {response_dim!r} is missing."
+            )
+        return df.loc[:, ["chain", "draw", "__obs__", response_str]].rename(
+            columns={response_str: response_names[0]}
+        )
+
     df_wide = df.pivot_table(
         index=["chain", "draw", "__obs__"], columns=response_dim, values=response_str
     ).reset_index()
 
     df_wide.columns.name = None
-    df_wide = df_wide.rename(columns={0: "rt", 1: "response"})
+    df_wide = df_wide.rename(columns=dict(enumerate(response_names)))
     return df_wide
 
 
diff --git a/tests/distribution_utils/test_distribution_utils.py b/tests/distribution_utils/test_distribution_utils.py
index 975b7535a..8b45eaf52 100644
--- a/tests/distribution_utils/test_distribution_utils.py
+++ b/tests/distribution_utils/test_distribution_utils.py
@@ -113,8 +113,15 @@ def test_make_distribution_forwards_choice_only_to_generated_rv(monkeypatch, rv)
     """Generated RVs must keep the choice-only support-shape contract."""
     captured = {}
 
-    def fake_make_hssm_rv(simulator_fun, list_params, lapse=None, is_choice_only=False):
+    def fake_make_hssm_rv(
+        simulator_fun,
+        list_params,
+        lapse=None,
+        is_choice_only=False,
+        expected_obs_dim=None,
+    ):
         captured["is_choice_only"] = is_choice_only
+        captured["expected_obs_dim"] = expected_obs_dim
 
         class FakeRV:
             def __call__(self):
@@ -132,6 +139,7 @@ def __call__(self):
     )
 
     assert captured["is_choice_only"] is True
+    assert captured["expected_obs_dim"] is None
 
 
 @pytest.mark.slow
diff --git a/tests/rl/test_choice_only_rl.py b/tests/rl/test_choice_only_rl.py
index 190973cb9..a9cc2f41f 100644
--- a/tests/rl/test_choice_only_rl.py
+++ b/tests/rl/test_choice_only_rl.py
@@ -266,6 +266,7 @@ def _make_shell_rlssm(config, lapse):
     model = _RLSSM.__new__(_RLSSM)
     model.list_params = [*config.list_params, "p_outlier"]
     model.model_config = config
+    model._obs_dim = len(config.response)
     model.bounds = dict(config.bounds)
     model.lapse = lapse
     model.data = pd.DataFrame({"response": [0, 1], "rt": [0.5, 0.6]})
diff --git a/tests/test_hssm.py b/tests/test_hssm.py
index 25768b862..402100cd0 100644
--- a/tests/test_hssm.py
+++ b/tests/test_hssm.py
@@ -756,7 +756,7 @@ def test_vi_idata_rejects_attached_approximation(data_ddm):
         _ = model.vi_idata
 
 
-def test_drop_parent_str_requires_datatree(data_ddm):
+def test_clean_predictive_datatree_requires_datatree(data_ddm):
     """Posterior cleanup rejects a missing trace object."""
     model = HSSM(data=data_ddm)
 
@@ -764,10 +764,10 @@ def test_drop_parent_str_requires_datatree(data_ddm):
         ValueError,
         match=r"Please provide a DataTree \(traces\) object\.",
     ):
-        model._drop_parent_str_from_datatree(None)
+        model._clean_predictive_datatree(None)
 
 
-def test_drop_parent_str_renames_response_mean(data_ddm):
+def test_clean_predictive_datatree_renames_response_mean(data_ddm):
     """Posterior cleanup restores the model's response-parameter name."""
     model = HSSM(data=data_ddm)
     traces = xr.DataTree.from_dict(
@@ -783,7 +783,7 @@ def test_drop_parent_str_renames_response_mean(data_ddm):
         }
     )
 
-    result = model._drop_parent_str_from_datatree(traces)
+    result = model._clean_predictive_datatree(traces)
 
     assert result is traces
     assert model._parent in result["posterior"].data_vars
@@ -800,6 +800,7 @@ def test_is_choice_only_and_deadline(data_ddm):
 
     model = HSSM(data=data_ddm, model="ddm", model_config=config_choice_only)
     assert model.model_config.is_choice_only
+    assert model._obs_dim == 1
     assert model.response_c == "response"
     assert model.response_str == "response"
 
@@ -816,3 +817,4 @@ def test_is_choice_only_and_deadline(data_ddm):
     assert len(model_with_deadline.response) == 2
     assert model_with_deadline.response_c == "c(response, deadline)"
     assert model_with_deadline.response_str == "response,deadline"
+    assert model_with_deadline._obs_dim == 1
diff --git a/tests/test_observation_dimensions.py b/tests/test_observation_dimensions.py
new file mode 100644
index 000000000..a094e43f2
--- /dev/null
+++ b/tests/test_observation_dimensions.py
@@ -0,0 +1,341 @@
+"""Tests for response-order-derived observation dimensions."""
+
+from collections.abc import Callable
+
+import cloudpickle
+import numpy as np
+import pandas as pd
+import pytensor.tensor as pt
+import pytest
+import xarray as xr
+
+import hssm
+from hssm.distribution_utils import make_distribution, make_hssm_rv
+
+
+def _synthetic_simulator(obs_dim: int, *, scalar_output: bool = False) -> Callable:
+    """Return a deterministic-contract simulator with seeded random values."""
+
+    def simulator(theta, random_state, n_replicas, **kwargs):
+        del kwargs
+        theta_array = np.asarray(theta)
+        n_rows = (
+            n_replicas if theta_array.ndim == 1 else theta_array.shape[0] * n_replicas
+        )
+        values = np.random.default_rng(random_state).normal(size=(n_rows, obs_dim))
+        return values[:, 0] if scalar_output else values
+
+    simulator.model_name = f"synthetic_{obs_dim}d"  # type: ignore[attr-defined]
+    simulator.choices = ()  # type: ignore[attr-defined]
+    simulator.obs_dim = obs_dim  # type: ignore[attr-defined]
+    return simulator
+
+
+def _synthetic_logp(data, v):
+    """Return one finite log-likelihood value per observation."""
+    if data.ndim == 1:
+        return -pt.square(data - v)
+    return -pt.sum(pt.square(data - pt.shape_padright(v)), axis=-1)
+
+
+def _synthetic_model(
+    configured_obs_dim: int,
+    *,
+    simulator_obs_dim: int | None = None,
+) -> hssm.HSSM:
+    """Build a small custom model with scalar physical response columns."""
+    simulator_obs_dim = simulator_obs_dim or configured_obs_dim
+    if configured_obs_dim == 1:
+        response = ("response",)
+        data = pd.DataFrame({"response": [0, 1, 0, 1, 0]})
+        response_domains = {"response": {"kind": "categorical", "values": (0, 1)}}
+    else:
+        response = ("rt",) + tuple(
+            f"response_{index}" for index in range(1, configured_obs_dim)
+        )
+        data = pd.DataFrame(
+            {
+                "rt": np.linspace(0.3, 0.7, 5),
+                **{column: np.linspace(-0.5, 0.5, 5) for column in response[1:]},
+            }
+        )
+        response_domains = {
+            column: {"kind": "continuous", "bounds": (-1.0, 1.0)}
+            for column in response[1:]
+        }
+
+    return hssm.HSSM(
+        data=data,
+        model="custom",
+        model_config=hssm.ModelConfig(
+            response=response,
+            response_domains=response_domains,  # type: ignore[arg-type]
+            list_params=["v"],
+            default_priors={"v": {"name": "Normal", "mu": 0.0, "sigma": 1.0}},
+            rv=_synthetic_simulator(
+                simulator_obs_dim,
+                scalar_output=configured_obs_dim == 1,
+            ),  # type: ignore[arg-type]
+        ),
+        loglik=_synthetic_logp,
+        loglik_kind="analytical",
+        p_outlier=None,
+        process_initvals=False,
+    )
+
+
+@pytest.mark.parametrize("obs_dim", [1, 2, 3, 4])
+def test_generated_rv_preserves_width_and_seed_stream(obs_dim):
+    """Generated RVs retain arbitrary widths and deterministic stream order."""
+    is_choice_only = obs_dim == 1
+    simulator = _synthetic_simulator(obs_dim, scalar_output=is_choice_only)
+    rv = make_hssm_rv(
+        simulator,
+        ["v"],
+        is_choice_only=is_choice_only,
+        expected_obs_dim=obs_dim,
+    )
+    legacy_rv = make_hssm_rv(
+        simulator,
+        ["v"],
+        is_choice_only=is_choice_only,
+    )
+    rng_a = np.random.default_rng(42)
+    rng_b = np.random.default_rng(42)
+    legacy_rng = np.random.default_rng(42)
+
+    first_a = rv.rng_fn(rng_a, 0.5, size=5)
+    second_a = rv.rng_fn(rng_a, 0.5, size=5)
+    first_b = rv.rng_fn(rng_b, 0.5, size=5)
+    second_b = rv.rng_fn(rng_b, 0.5, size=5)
+    legacy_first = legacy_rv.rng_fn(legacy_rng, 0.5, size=5)
+    legacy_second = legacy_rv.rng_fn(legacy_rng, 0.5, size=5)
+
+    expected_shape = (5,) if is_choice_only else (5, obs_dim)
+    assert first_a.shape == expected_shape
+    assert second_a.shape == expected_shape
+    np.testing.assert_array_equal(first_a, first_b)
+    np.testing.assert_array_equal(second_a, second_b)
+    np.testing.assert_array_equal(first_a, legacy_first)
+    np.testing.assert_array_equal(second_a, legacy_second)
+    assert not np.array_equal(first_a, second_a)
+
+
+def test_distribution_rejects_callable_width_mismatch():
+    """Configured response width must agree with generated-RV metadata."""
+    with pytest.raises(
+        ValueError,
+        match="simulator observation width 3.*configured response width 4",
+    ):
+        make_distribution(
+            rv=_synthetic_simulator(3),
+            loglik=_synthetic_logp,
+            list_params=["v"],
+            expected_obs_dim=4,
+        )
+
+
+def test_model_rejects_callable_width_mismatch_before_building_distribution():
+    """A custom model fails clearly when its simulator width is inconsistent."""
+    with pytest.raises(
+        ValueError,
+        match="simulator observation width 2.*configured response width 3",
+    ):
+        _synthetic_model(3, simulator_obs_dim=2)
+
+
+@pytest.mark.parametrize("obs_dim", [1, 2, 3, 4])
+def test_predictive_shapes_follow_configured_response_width(obs_dim):
+    """Prior and posterior predictions follow the physical response width."""
+    model = _synthetic_model(obs_dim)
+    prior = model.sample_prior_predictive(
+        draws=3,
+        random_seed=np.random.default_rng(41),
+    )
+    response_name = model.response_str
+    response_dim = f"{response_name}_dim"
+    prior_values = prior["prior_predictive"][response_name]
+    prior["posterior"] = prior["prior"]
+    posterior = model.sample_posterior_predictive(
+        prior,
+        draws=2,
+        safe_mode=False,
+        inplace=False,
+    )
+    assert posterior is not None
+    posterior_values = posterior["posterior_predictive"][response_name]
+
+    if obs_dim == 1:
+        assert not hasattr(model.family, "create_extra_pps_coord")
+        assert prior_values.dims == ("chain", "draw", "__obs__")
+        assert posterior_values.dims == ("chain", "draw", "__obs__")
+    else:
+        np.testing.assert_array_equal(
+            model.family.create_extra_pps_coord(), np.arange(obs_dim)
+        )
+        assert prior_values.dims[-1] == response_dim
+        assert posterior_values.dims[-1] == response_dim
+        np.testing.assert_array_equal(
+            posterior_values.coords[response_dim], np.arange(obs_dim)
+        )
+
+
+def test_safe_mode_preserves_width_four_draw_and_response_coordinates():
+    """Chunked prediction concatenates arbitrary-width draws without drift."""
+    model = _synthetic_model(4)
+    traces = model.sample_prior_predictive(
+        draws=12,
+        random_seed=np.random.default_rng(43),
+    )
+    traces["posterior"] = traces["prior"]
+
+    result = model.sample_posterior_predictive(
+        traces,
+        draws=12,
+        safe_mode=True,
+        inplace=False,
+    )
+
+    assert result is not None
+    response_name = model.response_str
+    response_dim = f"{response_name}_dim"
+    values = result["posterior_predictive"][response_name]
+    assert values.dims == ("chain", "draw", "__obs__", response_dim)
+    np.testing.assert_array_equal(values.coords["draw"], np.arange(12))
+    np.testing.assert_array_equal(values.coords[response_dim], np.arange(4))
+
+
+def test_choice_only_callable_width_mismatch_is_not_masked():
+    """Scalar support does not hide malformed callable width metadata."""
+    with pytest.raises(
+        ValueError,
+        match="simulator observation width 2.*configured response width 1",
+    ):
+        make_hssm_rv(
+            _synthetic_simulator(2, scalar_output=True),
+            ["v"],
+            is_choice_only=True,
+            expected_obs_dim=1,
+        )
+
+
+def test_choice_only_string_keeps_scalar_legacy_width():
+    """Legacy choice-only names retain scalar support despite wrapper metadata."""
+    rv = make_hssm_rv(
+        "choice_only_model",
+        ["beta"],
+        is_choice_only=True,
+        expected_obs_dim=1,
+    )
+
+    assert rv.signature == "()->()"
+
+
+@pytest.mark.parametrize("obs_dim", [1, 4])
+def test_sample_do_dataframe_uses_physical_response_order(obs_dim):
+    """Intervention samples expose each configured physical response column."""
+    model = _synthetic_model(obs_dim)
+    predictive = model.sample_do(params={"v": 0.5}, draws=3)
+    frame = hssm.utils.predictive_dt_to_dataframe(
+        predictive,
+        predictive_group="prior_predictive",
+        response_str=model.response_str,
+        response_dim=f"{model.response_str}_dim",
+    )
+
+    assert list(frame.columns) == ["chain", "draw", "__obs__", *model.response]
+    assert len(frame) == 3 * len(model.data)
+
+
+def test_scalar_dataframe_ignores_deadline_technical_suffix():
+    """Scalar predictions use the physical response name, not its suffix."""
+    response_str = "response,deadline"
+    predictive = xr.DataTree.from_dict(
+        {
+            "posterior_predictive": xr.Dataset(
+                {
+                    response_str: (
+                        ("chain", "draw", "__obs__"),
+                        np.arange(6).reshape(1, 2, 3),
+                    )
+                }
+            )
+        }
+    )
+
+    frame = hssm.utils.predictive_dt_to_dataframe(
+        predictive,
+        response_str=response_str,
+        response_dim=f"{response_str}_dim",
+    )
+
+    assert list(frame.columns) == ["chain", "draw", "__obs__", "response"]
+    np.testing.assert_array_equal(frame["response"], np.arange(6))
+
+
+def test_dataframe_rejects_missing_vector_response_coordinate():
+    """A differently named vector dimension is not mistaken for scalar output."""
+    response_str = "rt,response"
+    predictive = xr.DataTree.from_dict(
+        {
+            "posterior_predictive": xr.Dataset(
+                {
+                    response_str: (
+                        ("chain", "draw", "__obs__", "unexpected_dim"),
+                        np.arange(12).reshape(1, 2, 3, 2),
+                    )
+                }
+            )
+        }
+    )
+
+    with pytest.raises(
+        ValueError,
+        match="Predictive response coordinate 'rt,response_dim' is missing",
+    ):
+        hssm.utils.predictive_dt_to_dataframe(
+            predictive,
+            response_str=response_str,
+            response_dim=f"{response_str}_dim",
+        )
+
+
+def test_predictive_cleanup_uses_custom_response_name():
+    """Parent cleanup derives names from the physical response declaration."""
+    model = _synthetic_model(4)
+    response_mean = f"{model.response_str}_mean"
+    traces = xr.DataTree.from_dict(
+        {
+            "posterior": xr.Dataset(
+                {response_mean: (("chain", "draw"), np.array([[0.5]]))}
+            )
+        }
+    )
+
+    result = model._clean_predictive_datatree(traces)
+
+    assert model._parent in result["posterior"].data_vars
+    assert response_mean not in result["posterior"].data_vars
+
+
+def test_width_four_cloudpickle_round_trip_preserves_predictive_contract():
+    """Existing model reconstruction retains width, metadata, and seeded draws."""
+    model = _synthetic_model(4)
+    restored = cloudpickle.loads(cloudpickle.dumps(model))
+
+    assert restored._obs_dim == model._obs_dim == 4
+    assert restored.response == model.response
+    assert restored.response_domains == model.response_domains
+
+    expected = model.sample_prior_predictive(
+        draws=2,
+        random_seed=np.random.default_rng(91),
+    )
+    actual = restored.sample_prior_predictive(
+        draws=2,
+        random_seed=np.random.default_rng(91),
+    )
+    xr.testing.assert_equal(
+        actual["prior_predictive"].ds,
+        expected["prior_predictive"].ds,
+    )