From 57f499c735b7e3dcc99f18b6cc89ddcdf94fb581 Mon Sep 17 00:00:00 2001
From: Chris Nicholas <4948774+cnicholas@users.noreply.github.com>
Date: Sun, 18 Oct 2026 12:58:28 -0400
Subject: [PATCH 1/4] refactor(xbar-s): one per-group estimator for the global
 and stratified paths (#144)

The global and stratified Xbar/S builders carried the same aggregation, fit and
limit code twice. It now lives in five helpers (_aggregate_subgroups,
_fit_xbar_group, _apply_xbar_limits, _apply_s_limits, _subgroup_label_column),
the lane-divider label rule is its own function (_lane_labels), and
calculate_limits_vectorized accepts a sigma-hat per row. Phased Xbar/S (#144)
is N per-stage charts on one axis, each one exactly the per-stratum estimator;
with the estimator extracted the feature adds a loop, not arithmetic.

No behaviour change. validation/xbar_s_sweep.py (new) fingerprints 504 Xbar/S
requests on T100 PM SDS 1 and 3 and prints identical output before and after;
golden masters untouched; gate 1547 passed, 0 failed, 4 pending.
---
 processbehavior/analysis.py      | 352 +++++++++++++++----------------
 processbehavior/spc_constants.py |   7 +-
 tests/test_xbar_s_helpers.py     | 153 ++++++++++++++
 validation/README.txt            |   7 +
 validation/xbar_s_sweep.py       |  92 ++++++++
 5 files changed, 431 insertions(+), 180 deletions(-)
 create mode 100644 tests/test_xbar_s_helpers.py
 create mode 100644 validation/xbar_s_sweep.py

diff --git a/processbehavior/analysis.py b/processbehavior/analysis.py
index b04585e..f229664 100644
--- a/processbehavior/analysis.py
+++ b/processbehavior/analysis.py
@@ -211,6 +211,38 @@ def _stratum_row_positions(df: pd.DataFrame, stratify_col: str) -> dict:
     return df.groupby(stratify_col, sort=False, observed=True).indices
 
 
+@dataclass(frozen=True)
+class _GroupFit:
+    """What one group of Xbar/S subgroups (a chart, a stratum or a phase) needs for its limits.
+
+    Built by :meth:`Analysis._fit_xbar_group` from an aggregated subgroup frame; consumed by
+    :meth:`Analysis._apply_xbar_limits` and, for the companion S chart, by
+    :meth:`Analysis._apply_s_limits`.
+    """
+
+    out: pd.DataFrame  # the charted subgroups (n >= 1); N overwritten under n_mode='average'
+    fit: pd.Series  # mask over ``out``: the subgroups with n >= 2 that sigma-hat is estimated from
+    s_bar: float  # mean S over the fitted subgroups
+    n_to_use: str  # 'N' when the Xbar chart's subgroup sizes are constant, else 'n'
+    s_n_to_use: str  # the same decision for the S chart, which has no one-reading subgroups
+    n_max: int
+    n_avg: float | None  # mean fitted n under n_mode='average', else None
+
+
+def _lane_labels(df: pd.DataFrame, collapsed_vars: list[str]) -> np.ndarray:
+    """One display label per row of ``df``: the collapsed columns as ``astype(str)``, joined with ``'_'``.
+
+    Column-wise rather than ``.agg('_'.join, axis=1)``, which ran once per row (1.05s at 1M).
+    Deliberately NOT encode_rsg_series: that renders a datetime as str(Timestamp)
+    ('2024-01-01 00:00:00') where DataFrame.astype(str) gives '2024-01-01'. These are display
+    labels for a chart divider or a phase, not stratum identity, so they keep the astype
+    rendering — do not "unify" the two.
+    """
+    parts = [df[col].astype(str) for col in collapsed_vars]
+    combined = parts[0] if len(parts) == 1 else reduce(lambda left, right: left + '_' + right, parts)
+    return combined.to_numpy()
+
+
 def _stratum_rows(df: pd.DataFrame, positions: dict, stratum) -> pd.DataFrame:
     """The rows of `df` in `stratum`, in `df`'s order (a view; copy before writing)."""
     pos = positions.get(stratum)
@@ -1063,19 +1095,11 @@ def _calculate_lane_boundaries(self, df: pd.DataFrame, collapsed_vars: list[str]
         if not collapsed_vars or df.empty:
             return []
 
-        # Combine column-wise rather than with `.agg('_'.join, axis=1)`, which ran once per
-        # row (1.05s at 1M). Deliberately NOT encode_rsg_series: that renders a datetime as
-        # str(Timestamp) ('2024-01-01 00:00:00') where DataFrame.astype(str) gives
-        # '2024-01-01'. These are display labels for a chart divider, not stratum identity,
-        # so they keep the astype rendering — do not "unify" the two.
-        parts = [df[col].astype(str) for col in collapsed_vars]
-        combined = parts[0] if len(parts) == 1 else reduce(lambda left, right: left + '_' + right, parts)
-
         # Positional comparison against the previous row. The old form built a boolean mask,
         # took df.index[mask], dropped the first entry, then called df.index.get_loc() per
         # boundary — and get_loc is not O(1) on these frames (~1.3ms a call, ~1.8s total).
         # Comparing the array to itself shifted gives the 0-based positions directly.
-        labels = combined.to_numpy()
+        labels = _lane_labels(df, collapsed_vars)
         changed = np.empty(len(labels), dtype=bool)
         changed[0] = False  # the first row is a start, not a boundary
         changed[1:] = labels[1:] != labels[:-1]
@@ -1187,6 +1211,117 @@ def _xbar_limits(self, out: pd.DataFrame, s_bar: float, n_to_use: str) -> pd.Dat
             sigma_multiplier=self.request.n_sigma,
         )
 
+    def _aggregate_subgroups(
+        self,
+        df: pd.DataFrame,
+        groupby_cols: list[str],
+        value_col: str,
+        limits_col: str,
+        ybar_col: str | None = None,
+    ) -> pd.DataFrame:
+        """Aggregate ``df`` into Xbar/S subgroups, one row per ``groupby_cols`` combination.
+
+        Columns: the grouping columns; ``s`` and ``n`` from ``limits_col`` (for an effect
+        residual that is R2, so a reading without one leaves the subgroup one smaller, as VAS
+        treats it); ``xbar`` from ``value_col``; ``s_value`` when the two columns differ; and
+        ``N``, the largest ``n``. With ``ybar_col`` the residual calculator's cached cell mean
+        stands in for the subgroup mean.
+        """
+        agg_dict = {
+            's': pd.NamedAgg(column=limits_col, aggfunc='std'),
+            'mean': pd.NamedAgg(column=value_col, aggfunc='mean'),
+            # N_r counts the values S_r is computed from: a residual missing on a row (R2's
+            # first value) leaves that subgroup one smaller, as VAS treats it.
+            'n': pd.NamedAgg(column=limits_col, aggfunc='count'),
+        }
+        if limits_col != value_col:
+            agg_dict['s_value'] = pd.NamedAgg(column=value_col, aggfunc='std')
+        out = df.groupby(groupby_cols, as_index=False, observed=True).agg(**agg_dict)
+
+        if ybar_col and ybar_col in df.columns:
+            # Get unique Ybar values per group
+            ybar_by_group = df.groupby(groupby_cols, observed=True)[ybar_col].first().reset_index()
+            out = out.merge(ybar_by_group, on=groupby_cols, how='left')
+            out['xbar'] = out[ybar_col]
+            out = out.drop(columns=[ybar_col, 'mean'], errors='ignore')
+        else:
+            out = out.rename(columns={'mean': 'xbar'})
+
+        out['N'] = out['n'].max()
+        return out
+
+    def _fit_xbar_group(self, out: pd.DataFrame) -> _GroupFit | None:
+        """Decide which subgroups of ``out`` are charted and which estimate sigma-hat.
+
+        Every subgroup with a value is charted. A one-reading subgroup has no S, so sigma-hat
+        comes from the subgroups with n >= 2, and its limits are center ± 3·sigma-hat (10-12
+        manual Eqs 11-16, 11-17 with N_r = 1), as VAS draws them; the companion S chart keeps
+        only n >= 2. n counts the limits column, so a reading without one (R2's first value) is
+        not charted, as in VAS. Returns None when no subgroup has n >= 2; the caller decides
+        whether that raises (one chart) or is recorded (a stratum, a phase). Under
+        ``n_mode='average'`` the ``N`` column becomes the mean fitted n and both charts take
+        the constant-N path (Eq 11-20).
+        """
+        out = out[out['n'].ge(1)].copy()  # n=0: every charted value missing
+        fit = out['n'].ge(2)
+        if not fit.any():
+            return None
+
+        s_bar = out.loc[fit, 's'].mean()
+        # Constant or variable subgroup sizes: on the Xbar chart, and on the S chart, which has
+        # no one-reading subgroups
+        n_to_use, n_max = self._determine_n_to_use(out)
+        s_n_to_use, _ = self._determine_n_to_use(out[fit])
+
+        n_avg = None
+        if self.request.n_mode == 'average':
+            n_avg = out.loc[fit, 'n'].mean()
+            out['N'] = n_avg  # overwrite N column with average
+            n_to_use = s_n_to_use = 'N'  # force constant-N path
+        return _GroupFit(
+            out=out, fit=fit, s_bar=s_bar, n_to_use=n_to_use, s_n_to_use=s_n_to_use, n_max=n_max, n_avg=n_avg
+        )
+
+    def _apply_xbar_limits(self, out: pd.DataFrame, center: float, group: _GroupFit) -> pd.DataFrame:
+        """Centre, natural process limits, beyond-limits flag and rounding for one Xbar group."""
+        out['center'] = center
+        out[['lpl', 'upl']] = self._xbar_limits(out, group.s_bar, group.n_to_use)
+        out = self._add_beyond_limits_flag(out, value_col='xbar')
+        return out.round(self.spec.round_to)
+
+    def _apply_s_limits(self, out: pd.DataFrame, s_bar: float, n_to_use: str) -> pd.DataFrame:
+        """Centre (S-bar), B3/B4 limits per row, beyond-limits flag and rounding for one S group."""
+        out['center'] = s_bar
+        out[['lpl', 'upl']] = calculate_limits_vectorized(
+            'S',
+            sd=out['center'],
+            N=out[n_to_use],
+            sigma_multiplier=self.request.n_sigma,
+        )
+        out = self._add_beyond_limits_flag(out, value_col='s')
+        return out.round(self.spec.round_to)
+
+    def _subgroup_label_column(self, out: pd.DataFrame, groupby_cols: list[str]) -> tuple[pd.DataFrame, str | None]:
+        """Name the subgroup label column of an Xbar/S frame, building it when it spans several columns.
+
+        An explicit single ``by`` keeps the factor's own column (cast to str so Plotly treats
+        coded factors as categories); several ``by`` columns join into ``subgroup``; the
+        default grouping joins into ``group``; a single default column is used as is.
+        """
+        by = self.request.by
+        if by is not None and len(by) > 0:
+            if len(by) == 1:
+                out[groupby_cols[0]] = out[groupby_cols[0]].astype(str)
+                return out, groupby_cols[0]
+            out['subgroup'] = out[groupby_cols].astype(str).agg('_'.join, axis=1)
+            return out, 'subgroup'
+        if len(groupby_cols) > 1:
+            out['group'] = out[groupby_cols].astype(str).agg('_'.join, axis=1)
+            return out, 'group'
+        if groupby_cols:
+            return out, groupby_cols[0]
+        return out, None
+
     def _calculate_xbar(  # noqa: C901
         self, value_col: str | None = None, _return_intermediates: bool = False
     ) -> dict:
@@ -1269,81 +1404,37 @@ def _calculate_xbar(  # noqa: C901
                 out_dict['s_value'] = [df[value_col].std()]
             out = pd.DataFrame(out_dict)
         else:
-            # Group by specified columns
-            agg_dict = {
-                's': pd.NamedAgg(column=_limits_col, aggfunc='std'),
-                'mean': pd.NamedAgg(column=value_col, aggfunc='mean'),
-                # N_r counts the values S_r is computed from: a residual missing on a row (R2's
-                # first value) leaves that subgroup one smaller, as VAS treats it.
-                'n': pd.NamedAgg(column=_limits_col, aggfunc='count'),
-            }
-            if _limits_col != value_col:
-                agg_dict['s_value'] = pd.NamedAgg(column=value_col, aggfunc='std')
-            out = df.groupby(groupby_cols, as_index=False, observed=True).agg(**agg_dict)
-
-            # If we have pre-calculated Ybar, use it for xbar values
-            if ybar_col and ybar_col in df.columns:
-                # Get unique Ybar values per group
-                ybar_by_group = df.groupby(groupby_cols, observed=True)[ybar_col].first().reset_index()
-                out = out.merge(ybar_by_group, on=groupby_cols, how='left')
-                out['xbar'] = out[ybar_col]
-                out = out.drop(columns=[ybar_col, 'mean'], errors='ignore')
-            else:
-                # Rename columns for consistency: mean→xbar (values)
-                out = out.rename(columns={'mean': 'xbar'})
-
-            _N = out['n'].max()
-            out['N'] = _N
-
-        # Every subgroup with a value is charted. A one-reading subgroup has no S, so sigma-hat
-        # comes from the subgroups with n >= 2, and its limits are center ± 3·sigma-hat (10-1
-        # manual Eqs 11-16/17 with N_r = 1), as VAS draws them; the companion S chart keeps only
-        # n >= 2. n counts the limits column, so a reading without one (R2's first value) is not
-        # charted, as in VAS. Xbar stays in valid_charts on ADS 2 — the chart family is valid
-        # there (Xbar of a residual, or of a pooled subgroup, computes fine); it is only charting
-        # *the response* that has nothing to subgroup. See Study._response_pair_problem, which
-        # reports that in advance.
-        out = out[out['n'].ge(1)].copy()  # n=0: every charted value missing
-        fit = out['n'].ge(2)
-        if not fit.any():
+            out = self._aggregate_subgroups(df, groupby_cols, value_col, _limits_col, ybar_col)
+
+        # Which subgroups are charted and which estimate sigma-hat: _fit_xbar_group. Xbar stays
+        # in valid_charts on ADS 2 — the chart family is valid there (Xbar of a residual, or of
+        # a pooled subgroup, computes fine); it is only charting *the response* that has nothing
+        # to subgroup. See Study._response_pair_problem, which reports that in advance.
+        group = self._fit_xbar_group(out)
+        if group is None:
             sds = self.ads._ads_result.sds if self.ads._ads_result else '?'
             _raise_no_replicated_subgroups('Xbar', sds)
+        out, fit = group.out, group.fit
+        _S, n_to_use, s_n_to_use = group.s_bar, group.n_to_use, group.s_n_to_use
+        n_max, n_avg = group.n_max, group.n_avg
+        n_mode = self.request.n_mode
 
         # Use Bishop VAS grand mean (mean of cell means on value_col) as center
         _Xbar = _Ybar
-        _S = out.loc[fit, 's'].mean()
         _N = out.loc[fit, 'n'].max()
-        if 'N' not in out.columns:
-            out['N'] = _N
-
-        # Determine if subgroup sizes are constant or variable: on the Xbar chart, and on the
-        # S chart, which has no one-reading subgroups
-        n_to_use, n_max = self._determine_n_to_use(out)
-        s_n_to_use, _ = self._determine_n_to_use(out[fit])
-
-        # Override n_to_use if n_mode="average"
-        n_mode = self.request.n_mode
-        n_avg = None
-        if n_mode == 'average':
-            n_avg = out.loc[fit, 'n'].mean()
-            out['N'] = n_avg  # overwrite N column with average
-            n_to_use = s_n_to_use = 'N'  # force constant-N path
 
         # CALCULATE XBAR
-        xbar = out.copy()
         if self.calibration is not None:
             # Standards-given Xbar: center=mean, mean ± n_sigma·sigma/√N (per-row N).
+            xbar = out.copy()
             xbar[['center', 'lpl', 'upl']] = xbar.apply(
                 lambda row: self._calibrated_row('Xbar', N=row[n_to_use]),
                 axis=1,
             )
+            xbar = self._add_beyond_limits_flag(xbar, value_col='xbar')
+            xbar = xbar.round(spec.round_to)
         else:
-            xbar['center'] = _Xbar  # Add center column for Xbar chart
-            xbar[['lpl', 'upl']] = self._xbar_limits(xbar, _S, n_to_use)
-
-        # Detect beyond limits signals
-        xbar = self._add_beyond_limits_flag(xbar, value_col='xbar')
-        xbar = xbar.round(spec.round_to)
+            xbar = self._apply_xbar_limits(out.copy(), _Xbar, group)
 
         _center_stat = xbar['center'].iloc[0] if self.calibration is not None else _Xbar
         statistics['center'] = round(_center_stat, spec.round_to)
@@ -1357,27 +1448,7 @@ def _calculate_xbar(  # noqa: C901
             statistics['upl'] = None
             statistics['limits_vary'] = True
 
-        # Determine the grouping column for output
-        by = self.request.by
-        if by is not None and len(by) > 0:
-            # Explicit by list
-            if len(by) == 1:
-                # Single factor: preserve original column name for axis labeling.
-                # Cast to string so Plotly treats coded factors (e.g., 1,2,3) as categories.
-                xbar[groupby_cols[0]] = xbar[groupby_cols[0]].astype(str)
-                group_col = groupby_cols[0]
-            else:
-                xbar['subgroup'] = xbar[groupby_cols].astype(str).agg('_'.join, axis=1)
-                group_col = 'subgroup'
-        elif len(groupby_cols) > 1:
-            # by=None - create 'group' from multiple columns (original behavior)
-            xbar['group'] = xbar[groupby_cols].astype(str).agg('_'.join, axis=1)
-            group_col = 'group'
-        elif groupby_cols:
-            # Single groupby column
-            group_col = groupby_cols[0]
-        else:
-            group_col = None
+        xbar, group_col = self._subgroup_label_column(xbar, groupby_cols)
 
         # n is each subgroup's own count from the aggregation above: the size its limits
         # used. A pooled subgroup spans several cells, so no per-cell size could stand in (#143).
@@ -1470,46 +1541,21 @@ def _calculate_xbar_stratified(
             # Group by time within this stratum
             # Note: we always compute mean directly from filtered data, NOT using
             # Ybar_t (which is the marginal time mean across ALL factors).
-            agg_dict = {
-                's': pd.NamedAgg(column=_limits_col, aggfunc='std'),
-                'xbar': pd.NamedAgg(column=value_col, aggfunc='mean'),
-                'n': pd.NamedAgg(column=_limits_col, aggfunc='count'),
-            }
-            if _limits_col != value_col:
-                agg_dict['s_value'] = pd.NamedAgg(column=value_col, aggfunc='std')
-            out = sdf.groupby(groupby_cols, as_index=False, observed=True).agg(**agg_dict)
-
-            out['N'] = out['n'].max()
+            out = self._aggregate_subgroups(sdf, groupby_cols, value_col, _limits_col)
 
             # Centre first: Bishop's unweighted mean of the subgroup means, every
             # subgroup counted once — a one-observation subgroup is still a cell.
             _Xbar = out['xbar'].mean()
 
             # Every subgroup with a value is charted; sigma-hat comes from those with n >= 2
-            # (see _calculate_xbar)
-            out = out[out['n'].ge(1)].copy()
-            fit = out['n'].ge(2)
-            if not fit.any():
+            # (see _fit_xbar_group)
+            group = self._fit_xbar_group(out)
+            if group is None:
                 insufficient_strata.append(stratum)
                 continue
-
-            # Per-stratum statistics
-            _S = out.loc[fit, 's'].mean()
-            n_to_use, n_max = self._determine_n_to_use(out)
-            s_n_to_use, _ = self._determine_n_to_use(out[fit])
-
-            n_mode = self.request.n_mode
-            n_avg = None
-            if n_mode == 'average':
-                n_avg = out.loc[fit, 'n'].mean()
-                out['N'] = n_avg
-                n_to_use = s_n_to_use = 'N'
-
-            out['center'] = _Xbar
-            out[['lpl', 'upl']] = self._xbar_limits(out, _S, n_to_use)
-
-            out = self._add_beyond_limits_flag(out, value_col='xbar')
-            out = out.round(spec.round_to)
+            _S, n_to_use, s_n_to_use = group.s_bar, group.n_to_use, group.s_n_to_use
+            n_max, n_avg = group.n_max, group.n_avg
+            out = self._apply_xbar_limits(group.out, _Xbar, group)
 
             # Add stratum identifier
             out[stratify_col] = stratum
@@ -1678,16 +1724,7 @@ def _calculate_s_stratified(
                 if groupby_cols:
                     out[groupby_cols[0]] = out[groupby_cols[0]].astype(str)
 
-            out['center'] = _S
-            out[['lpl', 'upl']] = calculate_limits_vectorized(
-                'S',
-                sd=out['center'],
-                N=out[n_to_use],
-                sigma_multiplier=self.request.n_sigma,
-            )
-
-            out = self._add_beyond_limits_flag(out, value_col='s')
-            out = out.round(spec.round_to)
+            out = self._apply_s_limits(out, _S, n_to_use)
             out[stratify_col] = stratum
 
             all_s_frames.append(out)
@@ -1833,22 +1870,9 @@ def _calculate_s(  # noqa: C901
                 out['N'] = n_avg
                 n_to_use = 'N'
 
-            # Determine group_col
-            by = self.request.by
-            if by is not None and len(by) > 0:
-                if len(by) == 1:
-                    out[groupby_cols[0]] = out[groupby_cols[0]].astype(str)
-                    group_col = groupby_cols[0]
-                else:
-                    out['subgroup'] = out[groupby_cols].astype(str).agg('_'.join, axis=1)
-                    group_col = 'subgroup'
-            elif len(groupby_cols) > 1:
-                out['group'] = out[groupby_cols].astype(str).agg('_'.join, axis=1)
-                group_col = 'group'
-            elif groupby_cols:
-                group_col = groupby_cols[0]
-            else:
-                group_col = None
+        # Label the subgroups. Xbar labelled its own copy; the shared intermediate frame is
+        # still unlabelled, so the companion path labels it here too.
+        out, group_col = self._subgroup_label_column(out, groupby_cols)
 
         # Calculate S chart (common path for both precomputed and independent)
         out['groups'] = out['n'].count()
@@ -1860,41 +1884,15 @@ def _calculate_s(  # noqa: C901
                 lambda row: self._calibrated_row('S', N=row[n_to_use]),
                 axis=1,
             )
+            out = self._add_beyond_limits_flag(out, value_col='s')
+            out = out.round(spec.round_to)
         else:
-            out['center'] = _S
-            # Add limits columns
-            out[['lpl', 'upl']] = calculate_limits_vectorized(
-                'S',
-                sd=out['center'],
-                N=out[n_to_use],
-                sigma_multiplier=self.request.n_sigma,
-            )
-
-        # Detect beyond limits signals
-        out = self._add_beyond_limits_flag(out, value_col='s')
-
-        # Apply grouping column transformation if needed and not already done
-        if _precomputed is None:
-            # Already handled in independent path above
-            pass
-        else:
-            # For precomputed path, apply group_col transformation
-            by = self.request.by
-            if by is not None and len(by) > 0 and 'subgroup' not in out.columns:
-                if len(by) == 1:
-                    out[groupby_cols[0]] = out[groupby_cols[0]].astype(str)
-                    group_col = groupby_cols[0]
-                else:
-                    out['subgroup'] = out[groupby_cols].astype(str).agg('_'.join, axis=1)
-                    group_col = 'subgroup'
-            elif len(groupby_cols) > 1 and group_col == 'group' and 'group' not in out.columns:
-                out['group'] = out[groupby_cols].astype(str).agg('_'.join, axis=1)
+            out = self._apply_s_limits(out, _S, n_to_use)
 
         # n: each subgroup's own count, as on Xbar (#143).
         cols_to_keep = [group_col, 'n', 's', 'center', 'lpl', 'upl', 'beyond_limits']
         cols_to_keep = [c for c in cols_to_keep if c in out.columns]
         out = out[cols_to_keep]
-        out = out.round(spec.round_to)
 
         # Build statistics
         _S_stat = out['center'].iloc[0] if len(out) > 0 else None
diff --git a/processbehavior/spc_constants.py b/processbehavior/spc_constants.py
index adfe7af..3eb50ec 100644
--- a/processbehavior/spc_constants.py
+++ b/processbehavior/spc_constants.py
@@ -374,7 +374,8 @@ def calculate_limits_vectorized(
     ``DataFrame.apply(axis=1)``, each call constructing a ``pd.Series`` to carry two
     numbers. At 1M rows (~50,000 subgroups) that was ~2.5s per chart; this is ~0.001s.
 
-    ``mean``, ``sd``, ``N`` and ``mR`` accept scalars or array-likes and broadcast.
+    ``mean``, ``sd``, ``N``, ``mR`` and ``sigma`` accept scalars or array-likes and broadcast,
+    so one chart can carry a sigma-hat per row (a phased chart: one per phase).
     Subgroup-size constants (c4/b3/b4) are looked up once per *distinct* N and mapped, so
     an N of 1 still raises ``ValueError`` from :func:`c4` exactly as the scalar path does.
 
@@ -396,7 +397,7 @@ def calculate_limits_vectorized(
     not subgroups or strata. Do not assume the absence of a call here means a site was
     checked and rejected on correctness grounds.
     """
-    index = next((arg.index for arg in (mean, sd, N, mR) if isinstance(arg, pd.Series)), None)
+    index = next((arg.index for arg in (mean, sd, N, mR, sigma) if isinstance(arg, pd.Series)), None)
 
     def _per_n(func, sizes, *args):
         """Evaluate a subgroup-size constant once per distinct N, then broadcast."""
@@ -413,7 +414,7 @@ def _per_n(func, sizes, *args):
         sizes = np.asarray(N)
         if sigma is not None:
             # Eqs 11-16/11-17: half-width = multiplier * sigma_hat / sqrt(N_r)
-            half = (sigma_multiplier * float(sigma)) / np.sqrt(sizes)
+            half = (sigma_multiplier * np.asarray(sigma, dtype=float)) / np.sqrt(sizes)
         else:
             # Wd = S / c4(n), then half-width = (multiplier * Wd) / sqrt(n)
             half = (sigma_multiplier * (np.asarray(sd) / _per_n(c4, sizes))) / np.sqrt(sizes)
diff --git a/tests/test_xbar_s_helpers.py b/tests/test_xbar_s_helpers.py
new file mode 100644
index 0000000..a47f13e
--- /dev/null
+++ b/tests/test_xbar_s_helpers.py
@@ -0,0 +1,153 @@
+"""The Xbar/S helpers extracted for issue #144 (its first PR), tested on their own.
+
+The behaviour they carry is pinned end to end by the golden masters and by
+``validation/xbar_s_sweep.py`` (identical fingerprints before and after the extraction). These
+tests pin the pieces the phased Xbar/S chart will reuse: one label rule for lane dividers and
+phases, a per-row sigma-hat in ``calculate_limits_vectorized``, and the per-group fit.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+from conftest import detect_sds_for_test, make_spec
+
+from processbehavior.analysis import Analysis, _lane_labels
+from processbehavior.datasets.synthetic import make_design
+from processbehavior.formulation_spec import ChartRequest
+from processbehavior.spc_constants import b3, b4, calculate_limits_vectorized
+
+SPEC = {'rsg_vars': ['factor 1', 'factor 2'], 'response_var': 'y', 'time_var': 'time', 'round_to': 3}
+
+
+def _analysis(**request) -> Analysis:
+    df = make_design(1, K1=2, K2=2, T=3, n_min=2, n_max=3, seed=1)
+    sds = detect_sds_for_test(df, SPEC)
+    return Analysis(spec=make_spec(SPEC), request=ChartRequest(chart='Xbar', **request), sds=sds, df=df)
+
+
+class TestLaneLabels:
+    def test_single_column_keeps_astype_rendering(self):
+        df = pd.DataFrame({'d': pd.to_datetime(['2024-01-01', '2024-01-02', '2024-01-02'])})
+        # astype(str) renders a date without the time part; str(Timestamp) would not.
+        assert list(_lane_labels(df, ['d'])) == ['2024-01-01', '2024-01-02', '2024-01-02']
+
+    def test_several_columns_join_with_underscore(self):
+        df = pd.DataFrame({'a': [1, 1, 2], 'b': ['x', 'y', 'y']})
+        labels = _lane_labels(df, ['a', 'b'])
+        assert isinstance(labels, np.ndarray)
+        assert list(labels) == ['1_x', '1_y', '2_y']
+
+
+class TestVectorizedSigmaPerRow:
+    def test_series_sigma_matches_the_formula_row_by_row(self):
+        mean = pd.Series([10.0, 10.0, 10.0], index=[5, 6, 7])
+        sigma = pd.Series([1.0, 2.0, 0.5], index=[5, 6, 7])
+        n = pd.Series([4, 4, 1], index=[5, 6, 7])
+        limits = calculate_limits_vectorized('Xbar', mean=mean, sigma=sigma, N=n)
+        expected_half = 3 * sigma / np.sqrt(n)
+        pd.testing.assert_series_equal(limits['lpl'], mean - expected_half, check_names=False)
+        pd.testing.assert_series_equal(limits['upl'], mean + expected_half, check_names=False)
+        assert list(limits.index) == [5, 6, 7]
+
+    def test_scalar_sigma_equals_a_constant_series(self):
+        n = pd.Series([2, 3, 5, 1])
+        scalar = calculate_limits_vectorized('Xbar', mean=4.5, sigma=0.7, N=n)
+        series = calculate_limits_vectorized('Xbar', mean=4.5, sigma=pd.Series([0.7] * 4), N=n)
+        assert np.array_equal(scalar.to_numpy(), series.to_numpy())
+
+    def test_series_sigma_supplies_the_index_when_nothing_else_does(self):
+        sigma = pd.Series([1.0, 2.0], index=[30, 40])
+        limits = calculate_limits_vectorized('Xbar', mean=0.0, sigma=sigma, N=np.array([4, 4]))
+        assert list(limits.index) == [30, 40]
+
+
+class TestFitXbarGroup:
+    def test_none_when_no_subgroup_is_replicated(self):
+        out = pd.DataFrame({'s': [np.nan, np.nan], 'xbar': [1.0, 2.0], 'n': [1, 1], 'N': [1, 1]})
+        assert _analysis()._fit_xbar_group(out) is None
+
+    def test_one_reading_subgroups_are_charted_but_not_fitted(self):
+        out = pd.DataFrame(
+            {'s': [np.nan, 1.0, 3.0, 2.0], 'xbar': [9.0, 10.0, 11.0, 12.0], 'n': [1, 3, 3, 0], 'N': [3] * 4}
+        )
+        group = _analysis()._fit_xbar_group(out)
+        assert group is not None
+        assert list(group.out['n']) == [1, 3, 3]  # n = 0 dropped, n = 1 kept
+        assert list(group.fit) == [False, True, True]
+        assert group.s_bar == pytest.approx(2.0)
+        assert group.n_to_use == 'n'  # sizes vary on the Xbar chart (1, 3, 3) ...
+        assert group.s_n_to_use == 'N'  # ... but not on the S chart (3, 3)
+        assert group.n_max == 3
+        assert group.n_avg is None
+        assert list(group.out['N']) == [3, 3, 3]
+
+    def test_average_n_mode_forces_the_constant_n_path(self):
+        out = pd.DataFrame({'s': [1.0, 3.0, 2.0], 'xbar': [10.0, 11.0, 12.0], 'n': [2, 4, 3], 'N': [4] * 3})
+        group = _analysis(n_mode='average')._fit_xbar_group(out)
+        assert group is not None
+        assert group.n_avg == pytest.approx(3.0)
+        assert group.n_to_use == group.s_n_to_use == 'N'
+        assert list(group.out['N']) == [3.0, 3.0, 3.0]
+        assert list(group.out['n']) == [2, 4, 3]  # actual sizes survive on n
+
+
+class TestApplySLimits:
+    def test_b3_b4_per_row_around_s_bar(self):
+        out = pd.DataFrame({'s': [0.5, 4.0, 1.2], 'n': [3, 5, 5], 'N': [5] * 3})
+        result = _analysis()._apply_s_limits(out.copy(), s_bar=1.5, n_to_use='n')
+        assert list(result['center']) == [1.5, 1.5, 1.5]
+        assert list(result['lpl']) == [round(b3(3) * 1.5, 3), round(b3(5) * 1.5, 3), round(b3(5) * 1.5, 3)]
+        assert list(result['upl']) == [round(b4(3) * 1.5, 3), round(b4(5) * 1.5, 3), round(b4(5) * 1.5, 3)]
+        assert list(result['beyond_limits']) == [0, 1, 0]
+
+
+class TestAggregateSubgroups:
+    def test_cached_cell_mean_replaces_the_subgroup_mean(self):
+        df = pd.DataFrame(
+            {
+                'k': ['A', 'A', 'A', 'B', 'B'],
+                't': [1, 1, 2, 1, 1],
+                'y': [1.0, 3.0, 5.0, 2.0, 4.0],
+                'r2': [np.nan, 0.5, -0.5, 0.2, -0.2],
+                'Ybar_kt': [2.0, 2.0, 5.0, 3.0, 3.0],
+            }
+        )
+        analysis = _analysis()
+        plain = analysis._aggregate_subgroups(df, ['k', 't'], 'y', 'y')
+        assert list(plain.columns) == ['k', 't', 's', 'xbar', 'n', 'N']
+        assert list(plain['xbar']) == [2.0, 5.0, 3.0]
+        assert list(plain['n']) == [2, 1, 2]
+        assert list(plain['N']) == [2, 2, 2]
+
+        cached = analysis._aggregate_subgroups(df, ['k', 't'], 'y', 'y', ybar_col='Ybar_kt')
+        assert list(cached['xbar']) == [2.0, 5.0, 3.0]
+        assert 'Ybar_kt' not in cached.columns and 'mean' not in cached.columns
+
+        residual = analysis._aggregate_subgroups(df, ['k', 't'], 'y', 'r2')
+        # n and s count the limits column: the subgroup missing its R2 is one smaller
+        assert list(residual['n']) == [1, 1, 2]
+        assert 's_value' in residual.columns
+
+
+class TestSubgroupLabelColumn:
+    def test_default_grouping_joins_into_group(self):
+        out = pd.DataFrame({'factor 1': [1, 2], 'time': [10, 10], 'xbar': [0.0, 0.0]})
+        labelled, col = _analysis(by=None)._subgroup_label_column(out, ['factor 1', 'time'])
+        assert col == 'group' and list(labelled['group']) == ['1_10', '2_10']
+
+    def test_single_by_keeps_the_column_as_strings(self):
+        out = pd.DataFrame({'factor 1': [1, 2], 'xbar': [0.0, 0.0]})
+        labelled, col = _analysis(by=('factor 1',))._subgroup_label_column(out, ['factor 1'])
+        assert col == 'factor 1' and list(labelled['factor 1']) == ['1', '2']
+
+    def test_several_by_columns_join_into_subgroup(self):
+        out = pd.DataFrame({'factor 1': [1, 2], 'time': [10, 10], 'xbar': [0.0, 0.0]})
+        _, col = _analysis(by=('factor 1', 'time'))._subgroup_label_column(out, ['factor 1', 'time'])
+        assert col == 'subgroup'
+
+    def test_no_grouping_has_no_label_column(self):
+        out = pd.DataFrame({'group': ['All'], 'xbar': [0.0]})
+        _, col = _analysis(by=None)._subgroup_label_column(out, [])
+        assert col is None
diff --git a/validation/README.txt b/validation/README.txt
index 5496519..0013389 100644
--- a/validation/README.txt
+++ b/validation/README.txt
@@ -51,6 +51,13 @@ What lives here and why
   mR signals) and which do not (the X chart's "not one process" conclusion).
   Run from the repo root:
       python validation/mr_permutation_invariance.py validation/aco_per_capita_expenditure.csv
+- xbar_s_sweep.py - fingerprints every Xbar/S request on PM SDS 1 and PM SDS 3
+  of the T100 database (chart x companion x seven by= shapes x value x n_mode,
+  504 requests): one SHA-256 per request over the chart frames, statistics and
+  metadata, or the exception a refused request raises. Run it on two revisions
+  and diff the outputs; a refactor that changes nothing prints identical files
+  (issue #144, PR 1). Run from the repo root:
+      python validation/xbar_s_sweep.py > sweep.txt
 
 Removed in the 2026-08 cleanup (recoverable from git history):
 sds1-6_analysis.xlsx (Excel exports generated by the package itself - not
diff --git a/validation/xbar_s_sweep.py b/validation/xbar_s_sweep.py
new file mode 100644
index 0000000..cf1559e
--- /dev/null
+++ b/validation/xbar_s_sweep.py
@@ -0,0 +1,92 @@
+#!/usr/bin/env python3
+"""
+Fingerprint every Xbar/S request on Tom Bishop's reference data.
+
+One line per request: the request, then a SHA-256 over the chart frames (every value at 17
+significant digits), the statistics and the metadata, or the exception a refused request
+raises. Run it on two revisions and diff the outputs: a refactor that changes nothing prints
+identical files. Written for issue #144, whose first PR extracts the per-group Xbar/S
+estimator and must leave every existing chart bit-identical.
+
+Grid: response PM SDS 1 (ADS 1) and PM SDS 3 (ADS 3) x chart {Xbar, S} x companion {off, on}
+x by {None, [], one factor, both factors, time, one factor + time, both factors + time}
+x value {response, R5, R2} x n_mode {actual, average}.
+
+Usage, from the repo root:
+    python validation/xbar_s_sweep.py > sweep.txt
+"""
+
+from __future__ import annotations
+
+import hashlib
+import json
+import sys
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+import processbehavior as pb
+from processbehavior.exceptions import ProcessBehaviorError
+
+ROOT = Path(__file__).resolve().parent
+REFERENCE = ROOT / 'PBTESTDATABASE_T100.csv'
+FACTORS = ['FACTOR 1', 'FACTOR 2']
+TIME = 'PRODUCTION TIME'
+RESPONSES = ('PM SDS 1', 'PM SDS 3')
+BY_SHAPES = (None, [], ['FACTOR 1'], FACTORS, [TIME], ['FACTOR 1', TIME], FACTORS + [TIME])
+VALUES = (None, 'R5', 'R2')
+N_MODES = ('actual', 'average')
+
+
+def _normalise(obj):
+    """JSON-ready copy: numpy scalars to Python, NaN to the string 'nan', dict keys to str."""
+    if isinstance(obj, dict):
+        return {str(k): _normalise(v) for k, v in obj.items()}
+    if isinstance(obj, (list, tuple)):
+        return [_normalise(v) for v in obj]
+    if isinstance(obj, np.generic):
+        obj = obj.item()
+    if isinstance(obj, float) and np.isnan(obj):
+        return 'nan'
+    return obj
+
+
+def fingerprint(result) -> str:
+    digest = hashlib.sha256()
+    for name in sorted(result.charts):
+        payload = result.charts[name]
+        digest.update(name.encode())
+        digest.update(payload['data'].to_csv(index=True, float_format='%.17g').encode())
+        digest.update(json.dumps(_normalise(payload['statistics']), sort_keys=True).encode())
+        digest.update(json.dumps(_normalise(payload['metadata']), sort_keys=True).encode())
+        digest.update(json.dumps(_normalise(payload.get('strata')), sort_keys=True).encode())
+    return digest.hexdigest()
+
+
+def main() -> int:
+    df = pd.read_csv(REFERENCE, na_values=['*'])
+    for response in RESPONSES:
+        study = pb.formulate(df, response=response, factors=FACTORS, time=TIME)
+        for chart in ('Xbar', 'S'):
+            for companion in (False, True):
+                for by in BY_SHAPES:
+                    for value in VALUES:
+                        for n_mode in N_MODES:
+                            request = (
+                                f'{response} | chart={chart} companion={companion} by={by} '
+                                f'value={value} n_mode={n_mode}'
+                            )
+                            try:
+                                result = study.execute(
+                                    chart=chart, by=by, value=value, companion=companion, n_mode=n_mode
+                                )
+                            except ProcessBehaviorError as exc:
+                                print(f'{request} -> refused {type(exc).__name__}: {exc}')
+                                continue
+                            print(f'{request} -> {fingerprint(result)}')
+    return 0
+
+
+if __name__ == '__main__':
+    sys.exit(main())

From fb8f10deddee3be77b7e4ec32cafc225356fb000 Mon Sep 17 00:00:00 2001
From: Chris Nicholas <4948774+cnicholas@users.noreply.github.com>
Date: Sun, 18 Oct 2026 13:41:08 -0400
Subject: [PATCH 2/4] feat(charts): phased limits on Xbar and S charts (#144)

execute(chart='Xbar' | 'S', phased=True) gives each factor combination its own
centre line and limits on one chart, drawn as stepped lines with a labelled
divider at each change, as phased X/mR does for its lanes. Each phase is
computed as the per-condition chart of by=[time] computes one stratum: the
mean of its subgroup means as centre, its own S-bar and sigma-hat, its own
limits (10-12 manual Eqs 11-5, 11-8, 11-9, 11-16, 11-17). With other factors
pooled into the subgroup (by=['stage', time]) the centre is the mean of the
pooled subgroup averages, as Dr. Bishop specified on 2026-10-18. On Wheeler's
four-cavity thickness data the per-cavity centres and limits match Minitab's
staged Xbar-S chart to the digit.

New result.phase_statistics(chart); chart data and chart_table gain a phase
column; metadata carries phase_by, lane_boundaries, single_point_phases and
insufficient_phases. by=[time], residual values and calibrations are refused
with messages that say what to use instead. phased=False output is unchanged
(validation/xbar_s_sweep.py: 0 differing lines); gate 1547 / 0 / 4; full
suite 2628 passed.
---
 CHANGELOG.md                                  |  21 +
 docs/reference/api.md                         |  26 +-
 docs/user-guide/chart-types.md                |  51 ++
 processbehavior/analysis.py                   | 170 +++++-
 processbehavior/analysis_result.py            |  37 ++
 processbehavior/exceptions.py                 |   8 +-
 processbehavior/plotting/plotter.py           |   1 +
 processbehavior/plotting/renderers.py         |   5 +-
 processbehavior/result_tabular.py             |   5 +
 processbehavior/study.py                      |  94 ++-
 processbehavior/types.py                      |  19 +-
 .../data/wheeler_cavity_thickness.csv         | 401 +++++++++++++
 .../phased_xbar_s/S_data.parquet              | Bin 0 -> 6089 bytes
 .../phased_xbar_s/S_metadata.json             |  53 ++
 .../phased_xbar_s/S_statistics.json           |   7 +
 .../phased_xbar_s/Xbar_data.parquet           | Bin 0 -> 6294 bytes
 .../phased_xbar_s/Xbar_metadata.json          |  58 ++
 .../phased_xbar_s/Xbar_statistics.json        |   7 +
 tests/test_analysis_golden_master.py          |  17 +
 tests/test_phased_limits.py                   |  10 +-
 tests/test_phased_xbar_s.py                   | 536 ++++++++++++++++++
 tests/test_result_tabular.py                  |   2 +-
 tests/test_supports_calibration.py            |  15 +-
 23 files changed, 1509 insertions(+), 34 deletions(-)
 create mode 100644 tests/fixtures/data/wheeler_cavity_thickness.csv
 create mode 100644 tests/fixtures/golden_masters/phased_xbar_s/S_data.parquet
 create mode 100644 tests/fixtures/golden_masters/phased_xbar_s/S_metadata.json
 create mode 100644 tests/fixtures/golden_masters/phased_xbar_s/S_statistics.json
 create mode 100644 tests/fixtures/golden_masters/phased_xbar_s/Xbar_data.parquet
 create mode 100644 tests/fixtures/golden_masters/phased_xbar_s/Xbar_metadata.json
 create mode 100644 tests/fixtures/golden_masters/phased_xbar_s/Xbar_statistics.json
 create mode 100644 tests/test_phased_xbar_s.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5bd8d52..edd254d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,6 +14,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   single stream, ADS 2, X/mR, sigma from N − 1 moving ranges. An NA response keeps its place as an
   empty cell. The design report adds `Time: received order (no time column given)`. If the data
   already has an `Observation` column, formulate asks you to pass it as `time=` or rename it (#129).
+- **Phased limits on Xbar and S charts.** `execute(chart='Xbar', phased=True)` gives each factor
+  combination its own centre line and limits on one chart, drawn as stepped lines with a divider at
+  each change, the way phased X/mR already does for its lanes. Each phase is computed as a separate
+  Xbar/S analysis, exactly as `by=[time]` computes one chart per condition: the centre is the mean of
+  the phase's subgroup means, with its own S-bar and sigma-hat (10-12 manual Eqs 11-5, 11-8, 11-9,
+  11-16, 11-17). `by=['stage', time]` pools the other factors into the subgroup and phases on the
+  stage; the centre is then the mean of the pooled subgroup averages, as Dr. Bishop specified
+  (2026-10-18). On Wheeler's four-cavity thickness data the per-cavity centres and limits match
+  Minitab's staged Xbar-S chart. Residual values, calibrations and `by=[time]` are refused with
+  messages that say what to use instead; a phase with no replicated subgroup keeps its points without
+  limits and is named in `metadata['insufficient_phases']`. (#144)
+- **`result.phase_statistics(chart)`.** One row per phase of a phased Xbar/S chart: `phase`,
+  `subgroups`, `N`, `center`, `s_bar`, `sigma_hat`, `lpl`, `upl`. (#144)
 
 ### Changed
 - **`formulate(response=..., factors=[])` with no time is the same study.** It used to treat the
@@ -23,6 +36,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - **Xbar and S chart data carry `n`.** `get_chart('Xbar')` / `get_chart('S')`, and the Excel chart
   sheets, gain an `n` column after the subgroup label. Unstratified Xbar/S metadata names the label
   column as `group_col` (#143).
+- **Phased Xbar/S chart data and `chart_table` carry a `phase` column** after the subgroup label, and
+  the metadata carries `phase_by`, `lane_boundaries`, `single_point_phases` and `insufficient_phases`.
+  `n_avg` is not reported on a phased chart, whose `n_mode='average'` is applied per phase; the
+  per-phase mean size is `phase_statistics(chart)['N']`. (#144)
+- **Lane dividers on an Xbar or S chart are labelled with the phase.** These charts have no time
+  column and so no cell band under the axis, which is where X charts name their lanes. (#144)
 
 ### Fixed
 - **Xbar/S `chart_table` keeps the subgroup label and `n` on every `by=`.** With a `by=` other than
@@ -33,6 +52,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   (R1, R3, R4, R5) in a design with singletons: the stream's first cell, where R2 is undefined at
   the first value, showed its cell size, not the one-fewer values its S and limits use. Limits,
   centres and statistics are unchanged (#143).
+- **`execute()` no longer documents `phased=True` as requiring `by=[]`.** Stratified phasing of X/mR
+  (`by=['factor 1'], phased=True`) has been accepted since 75bcdd1; the docstring never caught up. (#144)
 
 ## [0.4.0] - 2026-10-14
 
diff --git a/docs/reference/api.md b/docs/reference/api.md
index 78456c3..06ababb 100644
--- a/docs/reference/api.md
+++ b/docs/reference/api.md
@@ -113,7 +113,7 @@ result = study.execute(
     recentered=False,    # bool — recentered residual (RCR*) form
     bins=None,           # int | None — Histogram bins
     companion=False,     # bool — also compute the companion chart (see below)
-    phased=False,        # bool — phase-wise limits
+    phased=False,        # bool — a center line and limits per phase (X/mR and Xbar/S)
     n_sigma=3.0,         # float — limit width
     n_mode='actual',     # 'actual' | 'average' — subgroup-size handling
     calibration=None,    # Calibration | str | None — standards-given limits
@@ -276,7 +276,8 @@ Every chart's statistics dict carries exactly these keys:
 
 When limits vary point-to-point (e.g. varying subgroup sizes with
 `n_mode='actual'`), `N`, `lpl` and `upl` are `None` and the dict additionally has
-`'limits_vary': True` — read per-point limits from `get_chart()`. Histogram adds
+`'limits_vary': True` — read per-point limits from `get_chart()`. On a phased chart
+`center` is `None` too. Histogram adds
 `{'mean', 'std', 'n'}` and has `lpl`/`upl` of `None`. For a stratified result,
 `get_statistics()` returns `{stratum: {...}}`; `statistics(chart, stratum=...)`
 picks one.
@@ -326,6 +327,27 @@ result.focus(stratum)          # -> FocusedAnalysisResult
 `False`. (`FocusedAnalysisResult` is created by `focus()`, not constructed
 directly, and is not a top-level export.)
 
+### Phases
+
+```python
+result.phase_statistics(chart)   # -> DataFrame, one row per phase (phased Xbar/S)
+```
+
+A chart executed with `phased=True` gives each phase — one combination of the grouping
+factors, a contiguous run of subgroups in time order — its own center line and limits.
+On Xbar/S the chart data gains a `phase` column (the lane label) next to the subgroup
+label, and `chart_table` shows it. `phase_statistics(chart)` returns the numbers behind
+the stepped lines: `phase`, `subgroups`, `N` (the subgroup size the phase's limits assume;
+`NaN` when sizes vary within it), `center`, `s_bar`, `sigma_hat`, `lpl`, `upl` (`NaN` when
+the limits step within the phase). It raises `ValidationError` on a chart that is not
+phased, and on phased X/mR, whose per-phase values are on each row of `get_chart()`.
+
+Phased metadata carries `phased=True`, `run_rules_applicable=False`, `phase_by` (the
+factors a phase is one combination of), `lane_boundaries` (the divider positions),
+`single_point_phases` and `insufficient_phases` (phases with no replicated subgroup,
+charted without limits). Phased Xbar/S refuses residual values, calibrations and
+`by=[time]`; see [Staged Xbar/S](../user-guide/chart-types.md#staged-xbars).
+
 ### Plotting, reporting, export
 
 ```python
diff --git a/docs/user-guide/chart-types.md b/docs/user-guide/chart-types.md
index d4cea28..4901d7f 100644
--- a/docs/user-guide/chart-types.md
+++ b/docs/user-guide/chart-types.md
@@ -240,6 +240,51 @@ Charted](residuals.md#how-effect-residuals-are-charted).
 2. **Then interpret Xbar** - Only meaningful if S is stable
 3. **Investigate signals** - What makes that subgroup different?
 
+### Staged Xbar/S
+
+A staged study records the stage as a factor, and the default Xbar chart already lays the
+subgroups out stage by stage: every subgroup of stage A, then every subgroup of stage B,
+each in time order. `phased=True` gives each stage its own center line and limits on that
+one chart, drawn as stepped lines with a divider at each stage change.
+
+```python
+study = pb.formulate(df, response='Thickness', time='Hour', factors=['Cavity'])
+staged = study.execute(chart='Xbar', phased=True, companion=True)
+staged.plot(chart='Xbar')
+staged.phase_statistics('Xbar')   # one row per stage: center, S-bar, sigma-hat, limits
+```
+
+Each stage is computed as a separate analysis, exactly as `by=['Hour']` computes one chart
+per stage: its center is the unweighted mean of its subgroup means, its S̄ is the mean of its
+subgroup standard deviations, and its σ̂ and limits follow from those (10-12 manual Eqs 11-5,
+11-8, 11-9, 11-16, 11-17). Nothing is pooled across stages. On Wheeler's four-cavity
+thickness data the per-cavity centers and limits match Minitab's staged Xbar-S chart to the
+digit; the per-cavity charts of `by=['Hour']` give the same numbers, one chart at a time.
+
+With more factors than the stage, `by=['Cavity', 'Hour']` pools the others into the
+subgroup and phases on `Cavity`. The center is then the mean of the pooled subgroup
+averages, per Dr. Bishop: ignore the other factor completely, group by stage and time, and
+chart per the usual rules. That equals charting a study formulated without the pooled
+factors.
+
+What the phased chart reports:
+
+- `get_chart('Xbar')` has a `phase` column and each row's own `center`, `lpl` and `upl`;
+  `chart_table` shows `phase` after `subgroup`.
+- `get_statistics` returns `None` for `center`, `N`, `lpl` and `upl` with
+  `limits_vary=True`, since no single value describes the chart; `phase_statistics(chart)`
+  carries the per-stage numbers.
+- `detect_signals` applies Rule 1 against each stage's own limits; the run rules are off,
+  as on phased X/mR, because the zones differ from stage to stage.
+- A stage with a single replicated subgroup takes its σ̂ from that subgroup. A stage with
+  no replicated subgroup keeps its points but has no limits, and is named in
+  `metadata['insufficient_phases']`.
+
+Stages are ordered as the factor's levels sort, so label them to sort (`1 before`,
+`2 after`, or an ordered pandas Categorical). `by=[time]` with `phased=True` is refused:
+that view is already one chart per stage. Residual values and calibrations are refused on
+phased Xbar/S.
+
 ## Understanding X and mR Charts
 
 ### The X (Individual) Chart
@@ -314,6 +359,10 @@ phased = study.execute(chart='X', by=[], phased=True, companion=True)
 phased.plot(chart='X')
 ```
 
+The same view exists for subgrouped data: `study.execute(chart='Xbar', phased=True,
+companion=True)` gives each factor combination its own center and limits on one Xbar/S
+chart (see [Staged Xbar/S](#staged-xbars)).
+
 **3. One subgroup at a time.** Full stratification: one chart per subgroup, each with
 its own limits.
 
@@ -350,6 +399,7 @@ only the lane order, and reports which results move.
 | "This transition is a signal" | `chart='mR', by=[]` | No. Which transitions are flagged depends on which lanes are adjacent. |
 | "Subgroup *k* differs in level" | `chart='Xbar', by=['lane'], value='R5', recentered=True` | The points, yes: each is a subgroup mean. The limits, not entirely in DS 2 and 3: they are set from R2, which differences the condition-and-period-adjusted series along the same lane-major sequence. The level difference between lanes no longer enters R2, but which observations are paired still depends on lane order. In DS 1, R2 is within-cell and carries no order. |
 | "Subgroup *k* is stable on its own" | `phased=True`, or `by=['lane']` | Yes. Only that subgroup's own ranges enter. |
+| "Stage *k* is stable on its own" (subgrouped data) | `chart='Xbar', phased=True` | Yes. Only that stage's own subgroups enter its center, S̄ and limits. |
 
 So a boundary signal should be read as evidence along the path, and claims about the
 subgroups themselves should come from the charts that have no path: the R5 chart for
@@ -401,6 +451,7 @@ VAS residuals (R1–R6) need both factors and time.
 |----------|-------|----------------|
 | Is this one process? | `chart='X', by=[], companion=True` | Subgroups do not join into one stable stream |
 | Each subgroup on its own? | `chart='X', by=[], phased=True` | Special cause within that subgroup |
+| Each stage on its own, with subgroups? | `chart='Xbar', phased=True, companion=True` | Special cause within that stage |
 | Are groups different? | `chart='Xbar'` | Group deviates from average |
 | Is variation stable? | `chart='S'` | Group has unusual variation |
 | Process over time? | `chart='X', by=[...]` | Special cause detected |
diff --git a/processbehavior/analysis.py b/processbehavior/analysis.py
index f229664..c8a22b8 100644
--- a/processbehavior/analysis.py
+++ b/processbehavior/analysis.py
@@ -42,6 +42,7 @@
 from .residual_calculator import calculate_r6_residuals, resolve_r6_groupby
 from .spc_constants import (
     REQUEST_RESIDUALS,
+    c4,
     calculate_limits,
     calculate_limits_vectorized,
     calibrated_limits,
@@ -705,7 +706,7 @@ def calibration_rejection_reason(self) -> str | None:
         """The path description a calibration would be refused for, or None if it applies.
 
         Calibration is applied on global and grouped (single-chart) Xbar/S/X/mR paths.
-        Stratified paths and phased X/mR are a documented follow-up.
+        Stratified paths and phased charts (X/mR and Xbar/S) are a documented follow-up.
 
         This is the predicate behind :meth:`Study.supports_calibration`. It asks the
         same two resolvers the chart code asks — :meth:`_resolve_by_grouping` for
@@ -724,7 +725,11 @@ def calibration_rejection_reason(self) -> str | None:
 
         if chart_type in ('Xbar', 'S'):
             _, _, stratify_by = self._resolve_by_grouping(value_col)
-            return f'stratified {chart_type}' if stratify_by else None
+            if stratify_by:
+                return f'stratified {chart_type}'
+            if self.request.phased:
+                return f'phased {chart_type}'
+            return None
 
         if chart_type in ('X', 'mR'):
             stratify_by, collapsed_factors = self._resolve_xmr_stratification()
@@ -1361,6 +1366,9 @@ def _calculate_xbar(  # noqa: C901
         # Resolve grouping based on `by` parameter
         groupby_cols, ybar_col, stratify_by = self._resolve_by_grouping(value_col)
 
+        if self.request.phased:
+            return {'Xbar': self._calculate_xbar_s_phased(value_col)['Xbar']}
+
         # Stratified path: separate chart per factor combination
         if stratify_by:
             self._reject_calibration('stratified Xbar')
@@ -1821,6 +1829,9 @@ def _calculate_s(  # noqa: C901
             # Resolve grouping based on `by` parameter
             groupby_cols, _, stratify_by = self._resolve_by_grouping(value_col)
 
+            if self.request.phased:
+                return {'S': self._calculate_xbar_s_phased(value_col)['S']}
+
             # Stratified path: separate chart per factor combination
             if stratify_by:
                 self._reject_calibration('stratified S')
@@ -1942,6 +1953,9 @@ def _calculate_xbar_s(self, value_col: str | None = None) -> dict[str, ChartPayl
         dict
             Combined chart data: {'Xbar': {...}, 'S': {...}}
         """
+        if self.request.phased:
+            return self._calculate_xbar_s_phased(value_col)
+
         # Calculate Xbar with intermediates for S reuse
         xbar_result = self._calculate_xbar(value_col, _return_intermediates=True)
 
@@ -1964,6 +1978,158 @@ def _calculate_xbar_s(self, value_col: str | None = None) -> dict[str, ChartPayl
         # Combine results
         return {**xbar_result, **s_result}
 
+    def _calculate_xbar_s_phased(self, value_col: str | None = None) -> dict[str, ChartPayload]:
+        """Xbar and S with a centre line and limits per phase (issue #144).
+
+        A phase is one combination of the grouping factors: a contiguous run of subgroups in
+        time order on the stage-by-stage layout that ``by=None`` already draws. Each phase is
+        computed as the per-condition chart of ``by=[time]`` computes one stratum — the
+        unweighted mean of its subgroup means as centre (10-12 manual Eq 10-9 within the
+        phase), its own S-bar and sigma-hat (Eq 11-5), its own limits (Eqs 11-8, 11-9, 11-16,
+        11-17) — and the phases are concatenated in display order with the stepped columns.
+        Nothing is pooled across phases: Dr. Bishop (2026-10-18), "treat each stage separately
+        as if it were a single analysis".
+
+        With factors besides the phase factors pooled into the subgroup (``by=['stage', time]``
+        in a study with more factors) the subgroups are the pooled ones and the centre is the
+        mean of their averages, per Bishop: "ignore the [other factor] completely and simply
+        group the data by [stage] and Time". That equals re-formulating the study without the
+        pooled factors and charting its per-condition charts.
+
+        Solo Xbar, solo S and the companion pair all come through here, so they agree by
+        construction. Residual values are refused in validation; a calibration here.
+        """
+        spec = self.spec
+        if value_col is None:
+            value_col = spec.response_var
+        self._reject_calibration(f'phased {self.analysis_type}')
+
+        df = self._source_frame().copy()
+        groupby_cols, _, _ = self._resolve_by_grouping(value_col)
+        time_var = spec.time_var
+        phase_by = [c for c in groupby_cols if c != time_var]
+        limits_col = self._resolve_limits_column(value_col, df)
+
+        # The subgroup means are computed from the data, as the per-condition charts do, not
+        # read from the cached cell means: a pooled subgroup has no cached mean anyway.
+        out = self._aggregate_subgroups(df, groupby_cols, value_col, limits_col)
+        out = out[out['n'].ge(1)].reset_index(drop=True)
+        lane_boundaries = self._calculate_lane_boundaries(out, phase_by)
+        labels = _lane_labels(out, phase_by)
+        starts = [0] + [b['position'] for b in lane_boundaries]
+        ends = starts[1:] + [len(out)]
+
+        xbar_frames: list[pd.DataFrame] = []
+        s_frames: list[pd.DataFrame] = []
+        xbar_rows: list[dict] = []
+        s_rows: list[dict] = []
+        insufficient: list[str] = []
+        single_point = 0
+        for start, end in zip(starts, ends, strict=True):
+            label = str(labels[start])
+            segment = out.iloc[start:end].copy()
+            segment['N'] = segment['n'].max()  # this phase's own N, as a stratum's would be
+            center = segment['xbar'].mean()
+            group = self._fit_xbar_group(segment)
+            if group is None:
+                # No replicated subgroup: the points stay, the limits are NaN (so nothing is
+                # flagged) and the phase is named in the metadata.
+                segment[['center', 'lpl', 'upl']] = np.nan
+                segment['beyond_limits'] = 0
+                segment['phase'] = label
+                xbar_frames.append(segment.round(spec.round_to))
+                insufficient.append(label)
+                xbar_rows.append({'phase': label, 'subgroups': len(segment)})
+                continue
+
+            xbar_seg = self._apply_xbar_limits(group.out, center, group)
+            xbar_seg['phase'] = label
+            xbar_frames.append(xbar_seg)
+            single_point += int(len(xbar_seg) == 1)
+
+            # The S chart works from the rounded Xbar slice with n >= 2, as the stratified
+            # companion path does, so the two S charts agree to the last digit.
+            s_seg = xbar_seg[xbar_seg['n'].ge(2)].drop(columns=['s_value'], errors='ignore').copy()
+            s_seg = self._apply_s_limits(s_seg, group.s_bar, group.s_n_to_use)
+            s_seg['phase'] = label
+            s_frames.append(s_seg)
+
+            n_used = group.out['N'].iloc[0]  # n_avg under n_mode='average', else the phase's N
+            sigma_hat = (
+                subgroup_sigma_hat(group.out.loc[group.fit, 's'], group.out.loc[group.fit, 'n'])
+                if group.n_to_use == 'n'
+                else group.s_bar / c4(int(n_used))
+            )
+            common = {'center': center, 's_bar': group.s_bar, 'sigma_hat': sigma_hat}
+            xbar_rows.append(
+                {
+                    'phase': label,
+                    'subgroups': len(xbar_seg),
+                    **self._phase_limits(xbar_seg, group.n_to_use, group.n_avg),
+                    **common,
+                }
+            )
+            s_rows.append(
+                {
+                    'phase': label,
+                    'subgroups': len(s_seg),
+                    **self._phase_limits(s_seg, group.s_n_to_use, group.n_avg),
+                    **{**common, 'center': group.s_bar},
+                }
+            )
+
+        if not s_frames:
+            sds = self.ads._ads_result.sds if self.ads._ads_result else '?'
+            _raise_no_replicated_subgroups(self.analysis_type, sds)
+
+        xbar_out = pd.concat(xbar_frames, ignore_index=True)
+        xbar_out, group_col = self._subgroup_label_column(xbar_out, groupby_cols)
+        s_out = pd.concat(s_frames, ignore_index=True)
+        s_out, _ = self._subgroup_label_column(s_out, groupby_cols)
+        s_lane_boundaries = self._calculate_lane_boundaries(s_out, ['phase'])
+
+        phase_cols = ['phase', 'subgroups', 'N', 'center', 's_bar', 'sigma_hat', 'lpl', 'upl']
+        statistics = {'N': None, 'center': None, 'lpl': None, 'upl': None, 'limits_vary': True}
+
+        def payload(chart: str, value_name: str, frame: pd.DataFrame, lanes: list[dict], rows: list[dict]):
+            cols = [group_col, 'phase', 'n', value_name, 'center', 'lpl', 'upl', 'beyond_limits']
+            phases = pd.DataFrame(rows, columns=phase_cols).round(spec.round_to)
+            return {
+                'data': frame[[c for c in cols if c in frame.columns]],
+                'statistics': dict(statistics),
+                'metadata': {
+                    'chart_type': chart,
+                    'value_col': value_name,
+                    'center_col': 'center',
+                    'n_sigma': self.request.n_sigma,
+                    'n_mode': self.request.n_mode,
+                    **({'group_col': group_col} if group_col else {}),
+                    'phased': True,
+                    'run_rules_applicable': False,
+                    'phase_by': list(phase_by),
+                    'lane_boundaries': lanes or None,
+                    'single_point_phases': single_point,
+                    'insufficient_phases': insufficient or None,
+                },
+                'phases': phases,
+            }
+
+        return {
+            'Xbar': payload('Xbar', 'xbar', xbar_out, lane_boundaries, xbar_rows),
+            'S': payload('S', 's', s_out, s_lane_boundaries, s_rows),
+        }
+
+    @staticmethod
+    def _phase_limits(frame: pd.DataFrame, n_to_use: str, n_avg: float | None) -> dict:
+        """The N, lpl and upl a phase reports: scalars when its limits are one pair, NaN when they step."""
+        if n_to_use != 'N':
+            return {'N': np.nan, 'lpl': np.nan, 'upl': np.nan}
+        return {
+            'N': n_avg if n_avg is not None else frame['N'].iloc[0],
+            'lpl': frame['lpl'].iloc[0],
+            'upl': frame['upl'].iloc[0],
+        }
+
     @staticmethod
     def _resolve_mr_source_column(value_col: str) -> str:
         """Map recentered column (RCR*) to its non-recentered residual (R*).
diff --git a/processbehavior/analysis_result.py b/processbehavior/analysis_result.py
index d0d7267..3a3f2b8 100644
--- a/processbehavior/analysis_result.py
+++ b/processbehavior/analysis_result.py
@@ -656,6 +656,43 @@ def statistics(self, chart: str, stratum: str | None = None) -> dict:
             )
         return dict(stats[stratum])
 
+    def phase_statistics(self, chart: str) -> pd.DataFrame:
+        """One row per phase of a phased chart: the numbers behind its stepped lines.
+
+        Columns: ``phase`` (the lane label, as in ``chart_table``), ``subgroups`` (points
+        charted in it), ``N`` (the subgroup size its limits assume: the common size, the mean
+        size under ``n_mode='average'``, or NaN when sizes vary within the phase), ``center``,
+        ``s_bar``, ``sigma_hat``, ``lpl`` and ``upl`` (NaN when the limits step within the
+        phase; ``get_chart`` carries them per row). A phase with no replicated subgroup keeps
+        its ``subgroups`` and is NaN elsewhere; ``metadata['insufficient_phases']`` names it.
+
+        Raises
+        ------
+        ValidationError
+            On a chart that is not phased, and on phased X/mR, whose per-phase values are on
+            each row of ``get_chart`` (a later release adds this table there; issue #144).
+        ChartNotAvailableError
+            If the chart is not in this result.
+        """
+        name = self._resolve_chart_name(chart)
+        if name not in self.charts:
+            raise ChartNotAvailableError(
+                f"Chart '{name}' not found. Available charts: {self.all_charts}", chart=name, available=self.all_charts
+            )
+        payload = self.charts[name]
+        if not payload.get('metadata', {}).get('phased'):
+            raise ValidationError(
+                f'Chart {chart!r} is not phased, so it has no per-phase statistics. '
+                'execute(..., phased=True) gives each phase its own centre line and limits.'
+            )
+        phases = payload.get('phases')
+        if phases is None:
+            raise ValidationError(
+                'phase_statistics is available on phased Xbar and S charts; a phased X or mR chart '
+                f"carries each phase's centre and limits on its rows (get_chart({chart!r}))."
+            )
+        return phases.copy()
+
     def get_residual(self, residual_type: str) -> pd.Series:
         """
         Get a specific residual series.
diff --git a/processbehavior/exceptions.py b/processbehavior/exceptions.py
index a1238e0..b16f996 100644
--- a/processbehavior/exceptions.py
+++ b/processbehavior/exceptions.py
@@ -132,9 +132,9 @@ class CalibrationNotSupportedError(ValidationError):
     """A calibration was supplied on a chart path that cannot apply it.
 
     Calibration is applied on global and grouped (single-chart) Xbar/S/X/mR paths.
-    Stratified paths (one chart per stratum) and phased X/mR are a documented
-    follow-up, so a calibration passed to one of those is refused rather than
-    silently ignored.
+    Stratified paths (one chart per stratum) and phased charts (X/mR and Xbar/S)
+    are a documented follow-up, so a calibration passed to one of those is refused
+    rather than silently ignored.
 
     Ask *before* executing with :meth:`Study.supports_calibration`, which answers
     the same question without raising — preferable to catching this, and far
@@ -146,7 +146,7 @@ class CalibrationNotSupportedError(ValidationError):
     Attributes
     ----------
     context : str | None
-        The unsupported path, e.g. ``'stratified Xbar'`` or ``'phased X/mR'``.
+        The unsupported path, e.g. ``'stratified Xbar'``, ``'phased X/mR'`` or ``'phased Xbar'``.
 
     Examples
     --------
diff --git a/processbehavior/plotting/plotter.py b/processbehavior/plotting/plotter.py
index 425ac93..7623279 100644
--- a/processbehavior/plotting/plotter.py
+++ b/processbehavior/plotting/plotter.py
@@ -1040,6 +1040,7 @@ def _get_x_column(self, data: pd.DataFrame) -> str | None:
             'groups',
             'n',
             'N',
+            'phase',
         }
         for col in data.columns:
             if col in metric_cols:
diff --git a/processbehavior/plotting/renderers.py b/processbehavior/plotting/renderers.py
index 54902f9..12bf9e1 100644
--- a/processbehavior/plotting/renderers.py
+++ b/processbehavior/plotting/renderers.py
@@ -254,7 +254,8 @@ def render_control_chart(
             # Single-chart mode: cell labels are rendered as a band below the
             # axis title by Plotter._apply_time_tick_labels, so suppress the
             # top-of-plot lane labels here. Faceted mode keeps top labels
-            # per panel.
+            # per panel. Phased Xbar/S has no time column and so no band: its
+            # divider labels are the only place the phase is named.
             add_lane_boundaries(
                 fig,
                 lane_boundaries,
@@ -262,7 +263,7 @@ def render_control_chart(
                 theme,
                 row=row,
                 col=col,
-                show_labels=ctx.is_faceted,
+                show_labels=ctx.is_faceted or spec.chart_type in ('Xbar', 'S'),
             )
 
 
diff --git a/processbehavior/result_tabular.py b/processbehavior/result_tabular.py
index 44ac81d..3729fb3 100644
--- a/processbehavior/result_tabular.py
+++ b/processbehavior/result_tabular.py
@@ -34,6 +34,7 @@
         "beyond_limits",
         "n",
         "N",
+        "phase",
         "obs_id",
         "x",
         "pull",
@@ -250,6 +251,10 @@ def _output_column_plan(
         output_cols.append(subgroup_col)
         renames[subgroup_col] = "subgroup"
 
+    # Phased Xbar/S: which phase's centre and limits the row uses.
+    if "phase" in chart_data.columns:
+        output_cols.append("phase")
+
     if "n" in chart_data.columns:
         output_cols.append("n")
 
diff --git a/processbehavior/study.py b/processbehavior/study.py
index 01ff909..48d460c 100644
--- a/processbehavior/study.py
+++ b/processbehavior/study.py
@@ -1919,16 +1919,28 @@ def execute(
             Not applicable to Histogram.
 
         phased : bool, default False
-            When True, computes per-phase center lines and control limits.
-            Each phase is a contiguous run of the same collapsed factor
-            combination (the groups demarcated by vertical lane boundaries).
-            Within each phase, the moving range is computed independently.
-
-            - phased=False (default): Global limits across entire chart
-            - phased=True: Per-phase limits (requires by=[])
-
-            Single-point phases yield zero-width limits because mR is
-            imputed as 0; count available in metadata['single_point_phases'].
+            When True, each phase gets its own center line and limits, drawn
+            as stepped lines on one chart. A phase is a contiguous run of one
+            factor combination (the lanes between the vertical dividers).
+
+            - X/mR: the collapsed factors define the phases (``by=[]``
+              collapses all of them; ``by=['factor 1']`` stratifies on it and
+              phases on the rest). Each phase's moving range is computed on
+              its own.
+            - Xbar/S: the grouping factors define the phases. ``by=None`` (or
+              all factors + time) phases on every factor combination;
+              ``by=['stage', time]`` pools the other factors into the
+              subgroup and phases on ``stage``. Each phase is computed as a
+              separate Xbar/S analysis: its center is the mean of its
+              subgroup means, with its own S-bar and sigma-hat (Dr. Bishop,
+              2026-10-18). ``by=[time]`` and residual values are refused.
+              Per-phase numbers: ``result.phase_statistics(chart)``.
+
+            Statistics report ``center``, ``N``, ``lpl`` and ``upl`` as None
+            with ``limits_vary=True``; run rules are off (Rule 1 only).
+            Single-point phases: X/mR impute mR as 0 (zero-width limits);
+            Xbar/S take sigma-hat from the one subgroup's S. Count in
+            metadata['single_point_phases'].
 
         n_sigma : float, default 3.0
             Sigma multiplier for Xbar/S chart control limits. Standard SPC
@@ -2261,7 +2273,7 @@ def _validate_execute_request(  # noqa: C901
         if recentered:
             self._validate_recentered(value)
         if phased:
-            self._validate_phased(base_chart, by_validated)
+            self._validate_phased(base_chart, by_validated, value)
 
         if companion and base_chart == 'Histogram':
             raise ValidationError('companion is not applicable to Histogram charts.')
@@ -2423,10 +2435,13 @@ def _validate_recentered(self, value: str | None) -> None:
         if value.upper() not in recenterable:
             raise ValidationError(f"recentered=True requires a residual value (R1-R6), got '{value}'")
 
-    def _validate_phased(self, base_chart: str, by_validated: list[str] | None) -> None:
-        """Validate phased=True requirements."""
+    def _validate_phased(self, base_chart: str, by_validated: list[str] | None, value: str | None = None) -> None:
+        """Validate phased=True: X/mR phase on the collapsed factors, Xbar/S on the grouping factors."""
+        if base_chart in ('Xbar', 'S'):
+            self._validate_phased_xbar_s(base_chart, by_validated, value)
+            return
         if base_chart not in ('X', 'mR'):
-            raise ValidationError(f"phased=True is only valid for X or mR charts, got '{base_chart}'.")
+            raise ValidationError(f"phased=True is only valid for X, mR, Xbar or S charts, got '{base_chart}'.")
         if by_validated is None:
             raise ValidationError('phased=True requires an explicit by= parameter.')
         if not self._spec.rsg_vars_list:
@@ -2445,6 +2460,57 @@ def _validate_phased(self, base_chart: str, by_validated: list[str] | None) -> N
                     'Remove factors from by= or set phased=False.'
                 )
 
+    def _validate_phased_xbar_s(self, base_chart: str, by_validated: list[str] | None, value: str | None) -> None:
+        """phased=True on Xbar/S: one centre line and one set of limits per factor combination.
+
+        A phase is a contiguous run of subgroups sharing the grouping factors, in time order, on
+        the stage-by-stage layout ``by=None`` already draws. So the request must chart the
+        response, have factors and a time variable, and group by some factors with time last;
+        ``by=None`` / ``[]`` phases on every factor, ``by=['stage', time]`` pools the others into
+        the subgroup. ``by=[time]`` is refused because that view is already one chart per
+        combination (issue #144).
+        """
+        factors = self._spec.rsg_vars_list
+        time_var = self._spec.time_var
+        if value is not None and value.upper().startswith('R'):
+            raise ValidationError(
+                f'phased=True on {base_chart} charts the response; value={value!r} is not supported. '
+                "Residual charts repeat the response's phases or are flat within a phase; for a phased "
+                f"residual use chart='X', by=[], value={value!r}, phased=True."
+            )
+        if not factors:
+            raise ValidationError(
+                'phased=True requires factors to define phases (rsg_key). '
+                'This study has no factors; phased limits do not apply.'
+            )
+        if not time_var:
+            raise ValidationError(
+                f'phased=True on {base_chart} requires a time variable: a phase is one factor '
+                "combination's subgroups in time order, and this study has none."
+            )
+        if not by_validated:
+            return  # None or []: phases are the factor combinations
+        by_list = list(by_validated)
+        if set(by_list) == set(factors) | {time_var}:
+            return  # all factors + time, any order: the same subgroups as by=None
+        if by_list == [time_var]:
+            raise ValidationError(
+                f'phased=True with by={by_list} on {base_chart}: that view already gives each factor '
+                'combination its own chart, centre line and limits (result.focus(k)). For the same '
+                'limits side by side on one chart use by=None with phased=True.'
+            )
+        if time_var not in by_list:
+            raise ValidationError(
+                f'phased=True on {base_chart} requires the time variable in by=: by={by_list} gives one '
+                f'point per level with no time order, so there are no phases. Use by={by_list + [time_var]}.'
+            )
+        if by_list[-1] != time_var:
+            ordered = [c for c in by_list if c != time_var] + [time_var]
+            raise ValidationError(
+                f'phased=True requires the time variable last in by= (got by={by_list}): phasing changes '
+                f'the limits, not the order of the points. Use by={ordered}.'
+            )
+
     def _parse_chart_request(self, chart: str) -> str:
         """
         Parse chart string to validate base chart type.
diff --git a/processbehavior/types.py b/processbehavior/types.py
index cc84d80..0325768 100644
--- a/processbehavior/types.py
+++ b/processbehavior/types.py
@@ -61,7 +61,9 @@ class ChartStatistics(TypedDict):
     When the control limits and/or subgroup size vary across subgroups
     (e.g. Xbar/S with ``n_mode='actual'`` on unbalanced data),
     ``N``, ``lpl``, and ``upl`` are ``None`` and the optional
-    ``limits_vary`` flag is ``True``. ``center`` remains a single scalar.
+    ``limits_vary`` flag is ``True``; ``center`` stays a scalar. On a phased
+    chart ``center`` is ``None`` too, and each row's ``center`` column
+    carries its phase's value.
     """
 
     N: int | None
@@ -95,10 +97,15 @@ class ChartMetadata(TypedDict, total=False):
       charts only.
     - ``group_col``: unstratified Xbar/S — the chart-data column that labels
       each subgroup (``group``, ``subgroup``, or the factor's own name).
-    - ``lane_boundaries``: stratified X/mR charts with collapsed factors.
+    - ``lane_boundaries``: X/mR charts with collapsed factors and phased
+      Xbar/S — a flat list of ``{'position', 'label', 'variables'}`` dicts,
+      or one such list per stratum on a stratified chart.
     - ``insufficient_strata``: mR/r charts that dropped first-per-stratum.
     - ``phased`` / ``single_point_phases``: charts run with
       ``phased=True``.
+    - ``phase_by`` / ``insufficient_phases``: phased Xbar/S — the factors a
+      phase is one combination of, and the phases with no replicated
+      subgroup (charted without limits).
     - ``run_rules_applicable``: present on charts where Western Electric
       rules apply (or explicitly don't).
     - ``limits_source`` / ``calibration``: reserved for the standards-
@@ -112,10 +119,12 @@ class ChartMetadata(TypedDict, total=False):
     stratify_col: str | None
     stratify_by: list[str]
     group_col: str  # Unstratified Xbar/S: the subgroup label column
-    lane_boundaries: dict[str, list[float]] | None
+    lane_boundaries: list[dict] | dict[str, list[dict]] | None
     insufficient_strata: list[str] | None
     phased: bool
     single_point_phases: int
+    phase_by: list[str]  # Phased Xbar/S: the factors a phase is one combination of
+    insufficient_phases: list[str] | None  # Phased Xbar/S: phases charted without limits
     run_rules_applicable: bool
     bins: int  # Histogram only
     focused_stratum: str  # After AnalysisResult.focus()
@@ -134,13 +143,15 @@ class ChartPayload(TypedDict):
     """The dict structure stored in ``AnalysisResult.charts[name]``.
 
     Required: ``data``, ``statistics``, ``metadata``.
-    Optional: ``strata`` (present on stratified charts).
+    Optional: ``strata`` (present on stratified charts); ``phases`` (phased
+    Xbar/S: one row per phase, read through ``AnalysisResult.phase_statistics``).
     """
 
     data: pd.DataFrame
     statistics: ChartStatisticsValue
     metadata: ChartMetadata
     strata: NotRequired[list[str]]
+    phases: NotRequired[pd.DataFrame]
 
 
 # Type alias for the full charts dict on AnalysisResult.
diff --git a/tests/fixtures/data/wheeler_cavity_thickness.csv b/tests/fixtures/data/wheeler_cavity_thickness.csv
new file mode 100644
index 0000000..01d62fc
--- /dev/null
+++ b/tests/fixtures/data/wheeler_cavity_thickness.csv
@@ -0,0 +1,401 @@
+"Hour","Cavity","Cycle","Thickness","CAVITY LABEL","CAVITY MEAN","CAVITY STDEV","CAVITY ID","RSG LABEL-HOUR","XBAR CHART CL","S CHART CL","XBAR CHART LBL","XBAR CHART UBL","S CHART LBL","S CHART UBL","RSG SAMPLE SIZE"
+1,"A",1,15,"A",14.65,1.4705279674745848,"A",1,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+1,"A",2,16,"B",8.97,1.4948554613643728,"A",2,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+1,"A",3,17,"C",7.47,1.6805390209912117,"A",3,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+1,"A",4,16,"D",7.86,1.4694560741275022,"A",4,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+1,"A",5,18,,,,"A",5,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+2,"A",1,13,,,,"A",6,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+2,"A",2,18,,,,"A",7,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+2,"A",3,15,,,,"A",8,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+2,"A",4,15,,,,"A",9,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+2,"A",5,15,,,,"A",10,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+3,"A",1,14,,,,"A",11,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+3,"A",2,14,,,,"A",12,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+3,"A",3,18,,,,"A",13,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+3,"A",4,14,,,,"A",14,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+3,"A",5,14,,,,"A",15,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+4,"A",1,14,,,,"A",16,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+4,"A",2,15,,,,"A",17,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+4,"A",3,15,,,,"A",18,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+4,"A",4,15,,,,"A",19,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+4,"A",5,14,,,,"A",20,14.65,1.3822751182152948,12.677079701147338,16.622920298852662,0,2.8875697758124175,5
+5,"A",1,12,,,,"B",1,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+5,"A",2,13,,,,"B",2,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+5,"A",3,13,,,,"B",3,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+5,"A",4,12,,,,"B",4,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+5,"A",5,13,,,,"B",5,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+6,"A",1,14,,,,"B",6,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+6,"A",2,15,,,,"B",7,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+6,"A",3,17,,,,"B",8,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+6,"A",4,14,,,,"B",9,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+6,"A",5,13,,,,"B",10,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+7,"A",1,12,,,,"B",11,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+7,"A",2,12,,,,"B",12,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+7,"A",3,13,,,,"B",13,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+7,"A",4,13,,,,"B",14,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+7,"A",5,11,,,,"B",15,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+8,"A",1,13,,,,"B",16,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+8,"A",2,15,,,,"B",17,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+8,"A",3,16,,,,"B",18,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+8,"A",4,14,,,,"B",19,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+8,"A",5,13,,,,"B",20,8.97,1.4051426122284405,6.964440943111471,10.97555905688853,0,2.9353399220667944,5
+9,"A",1,15,,,,"C",1,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+9,"A",2,16,,,,"C",2,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+9,"A",3,17,,,,"C",3,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+9,"A",4,14,,,,"C",4,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+9,"A",5,13,,,,"C",5,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+10,"A",1,13,,,,"C",6,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+10,"A",2,13,,,,"C",7,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+10,"A",3,14,,,,"C",8,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+10,"A",4,13,,,,"C",9,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+10,"A",5,13,,,,"C",10,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+11,"A",1,13,,,,"C",11,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+11,"A",2,16,,,,"C",12,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+11,"A",3,13,,,,"C",13,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+11,"A",4,16,,,,"C",14,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+11,"A",5,18,,,,"C",15,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+12,"A",1,16,,,,"C",16,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+12,"A",2,18,,,,"C",17,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+12,"A",3,18,,,,"C",18,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+12,"A",4,16,,,,"C",19,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+12,"A",5,13,,,,"C",20,7.47,1.5796824849889777,5.215320306133623,9.724679693866378,0,3.2999533442545568,5
+13,"A",1,13,,,,"D",1,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+13,"A",2,18,,,,"D",2,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+13,"A",3,18,,,,"D",3,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+13,"A",4,14,,,,"D",4,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+13,"A",5,13,,,,"D",5,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+14,"A",1,13,,,,"D",6,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+14,"A",2,13,,,,"D",7,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+14,"A",3,14,,,,"D",8,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+14,"A",4,14,,,,"D",9,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+14,"A",5,14,,,,"D",10,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+15,"A",1,13,,,,"D",11,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+15,"A",2,15,,,,"D",12,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+15,"A",3,18,,,,"D",13,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+15,"A",4,15,,,,"D",14,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+15,"A",5,14,,,,"D",15,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+16,"A",1,15,,,,"D",16,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+16,"A",2,14,,,,"D",17,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+16,"A",3,15,,,,"D",18,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+16,"A",4,14,,,,"D",19,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+16,"A",5,12,,,,"D",20,7.86,1.3812675539011,5.888517796980562,9.831482203019439,0,2.8854649761075564,5
+17,"A",1,14,,,,,,,,,,,,
+17,"A",2,16,,,,,,,,,,,,
+17,"A",3,14,,,,,,,,,,,,
+17,"A",4,13,,,,,,,,,,,,
+17,"A",5,14,,,,,,,,,,,,
+18,"A",1,12,,,,,,,,,,,,
+18,"A",2,15,,,,,,,,,,,,
+18,"A",3,15,,,,,,,,,,,,
+18,"A",4,15,,,,,,,,,,,,
+18,"A",5,14,,,,,,,,,,,,
+19,"A",1,17,,,,,,,,,,,,
+19,"A",2,15,,,,,,,,,,,,
+19,"A",3,18,,,,,,,,,,,,
+19,"A",4,18,,,,,,,,,,,,
+19,"A",5,17,,,,,,,,,,,,
+20,"A",1,18,,,,,,,,,,,,
+20,"A",2,16,,,,,,,,,,,,
+20,"A",3,17,,,,,,,,,,,,
+20,"A",4,14,,,,,,,,,,,,
+20,"A",5,13,,,,,,,,,,,,
+1,"B",1,10,,,,,,,,,,,,
+1,"B",2,13,,,,,,,,,,,,
+1,"B",3,11,,,,,,,,,,,,
+1,"B",4,10,,,,,,,,,,,,
+1,"B",5,10,,,,,,,,,,,,
+2,"B",1,9,,,,,,,,,,,,
+2,"B",2,10,,,,,,,,,,,,
+2,"B",3,11,,,,,,,,,,,,
+2,"B",4,8,,,,,,,,,,,,
+2,"B",5,9,,,,,,,,,,,,
+3,"B",1,9,,,,,,,,,,,,
+3,"B",2,10,,,,,,,,,,,,
+3,"B",3,10,,,,,,,,,,,,
+3,"B",4,12,,,,,,,,,,,,
+3,"B",5,9,,,,,,,,,,,,
+4,"B",1,10,,,,,,,,,,,,
+4,"B",2,10,,,,,,,,,,,,
+4,"B",3,11,,,,,,,,,,,,
+4,"B",4,11,,,,,,,,,,,,
+4,"B",5,10,,,,,,,,,,,,
+5,"B",1,8,,,,,,,,,,,,
+5,"B",2,7,,,,,,,,,,,,
+5,"B",3,8,,,,,,,,,,,,
+5,"B",4,7,,,,,,,,,,,,
+5,"B",5,7,,,,,,,,,,,,
+6,"B",1,7,,,,,,,,,,,,
+6,"B",2,11,,,,,,,,,,,,
+6,"B",3,12,,,,,,,,,,,,
+6,"B",4,8,,,,,,,,,,,,
+6,"B",5,6,,,,,,,,,,,,
+7,"B",1,6,,,,,,,,,,,,
+7,"B",2,6,,,,,,,,,,,,
+7,"B",3,6,,,,,,,,,,,,
+7,"B",4,6,,,,,,,,,,,,
+7,"B",5,7,,,,,,,,,,,,
+8,"B",1,9,,,,,,,,,,,,
+8,"B",2,8,,,,,,,,,,,,
+8,"B",3,7,,,,,,,,,,,,
+8,"B",4,7,,,,,,,,,,,,
+8,"B",5,8,,,,,,,,,,,,
+9,"B",1,11,,,,,,,,,,,,
+9,"B",2,13,,,,,,,,,,,,
+9,"B",3,11,,,,,,,,,,,,
+9,"B",4,13,,,,,,,,,,,,
+9,"B",5,8,,,,,,,,,,,,
+10,"B",1,7,,,,,,,,,,,,
+10,"B",2,8,,,,,,,,,,,,
+10,"B",3,7,,,,,,,,,,,,
+10,"B",4,6,,,,,,,,,,,,
+10,"B",5,7,,,,,,,,,,,,
+11,"B",1,10,,,,,,,,,,,,
+11,"B",2,8,,,,,,,,,,,,
+11,"B",3,10,,,,,,,,,,,,
+11,"B",4,10,,,,,,,,,,,,
+11,"B",5,11,,,,,,,,,,,,
+12,"B",1,13,,,,,,,,,,,,
+12,"B",2,8,,,,,,,,,,,,
+12,"B",3,8,,,,,,,,,,,,
+12,"B",4,10,,,,,,,,,,,,
+12,"B",5,7,,,,,,,,,,,,
+13,"B",1,11,,,,,,,,,,,,
+13,"B",2,10,,,,,,,,,,,,
+13,"B",3,9,,,,,,,,,,,,
+13,"B",4,10,,,,,,,,,,,,
+13,"B",5,10,,,,,,,,,,,,
+14,"B",1,7,,,,,,,,,,,,
+14,"B",2,8,,,,,,,,,,,,
+14,"B",3,9,,,,,,,,,,,,
+14,"B",4,8,,,,,,,,,,,,
+14,"B",5,8,,,,,,,,,,,,
+15,"B",1,7,,,,,,,,,,,,
+15,"B",2,12,,,,,,,,,,,,
+15,"B",3,10,,,,,,,,,,,,
+15,"B",4,10,,,,,,,,,,,,
+15,"B",5,7,,,,,,,,,,,,
+16,"B",1,6,,,,,,,,,,,,
+16,"B",2,7,,,,,,,,,,,,
+16,"B",3,8,,,,,,,,,,,,
+16,"B",4,6,,,,,,,,,,,,
+16,"B",5,10,,,,,,,,,,,,
+17,"B",1,9,,,,,,,,,,,,
+17,"B",2,10,,,,,,,,,,,,
+17,"B",3,13,,,,,,,,,,,,
+17,"B",4,7,,,,,,,,,,,,
+17,"B",5,6,,,,,,,,,,,,
+18,"B",1,9,,,,,,,,,,,,
+18,"B",2,7,,,,,,,,,,,,
+18,"B",3,8,,,,,,,,,,,,
+18,"B",4,10,,,,,,,,,,,,
+18,"B",5,9,,,,,,,,,,,,
+19,"B",1,10,,,,,,,,,,,,
+19,"B",2,10,,,,,,,,,,,,
+19,"B",3,12,,,,,,,,,,,,
+19,"B",4,12,,,,,,,,,,,,
+19,"B",5,6,,,,,,,,,,,,
+20,"B",1,8,,,,,,,,,,,,
+20,"B",2,10,,,,,,,,,,,,
+20,"B",3,11,,,,,,,,,,,,
+20,"B",4,10,,,,,,,,,,,,
+20,"B",5,7,,,,,,,,,,,,
+1,"C",1,7,,,,,,,,,,,,
+1,"C",2,8,,,,,,,,,,,,
+1,"C",3,10,,,,,,,,,,,,
+1,"C",4,7,,,,,,,,,,,,
+1,"C",5,10,,,,,,,,,,,,
+2,"C",1,7,,,,,,,,,,,,
+2,"C",2,11,,,,,,,,,,,,
+2,"C",3,10,,,,,,,,,,,,
+2,"C",4,10,,,,,,,,,,,,
+2,"C",5,9,,,,,,,,,,,,
+3,"C",1,8,,,,,,,,,,,,
+3,"C",2,10,,,,,,,,,,,,
+3,"C",3,9,,,,,,,,,,,,
+3,"C",4,10,,,,,,,,,,,,
+3,"C",5,8,,,,,,,,,,,,
+4,"C",1,7,,,,,,,,,,,,
+4,"C",2,9,,,,,,,,,,,,
+4,"C",3,12,,,,,,,,,,,,
+4,"C",4,10,,,,,,,,,,,,
+4,"C",5,10,,,,,,,,,,,,
+5,"C",1,5,,,,,,,,,,,,
+5,"C",2,6,,,,,,,,,,,,
+5,"C",3,8,,,,,,,,,,,,
+5,"C",4,5,,,,,,,,,,,,
+5,"C",5,4,,,,,,,,,,,,
+6,"C",1,6,,,,,,,,,,,,
+6,"C",2,4,,,,,,,,,,,,
+6,"C",3,7,,,,,,,,,,,,
+6,"C",4,6,,,,,,,,,,,,
+6,"C",5,5,,,,,,,,,,,,
+7,"C",1,4,,,,,,,,,,,,
+7,"C",2,5,,,,,,,,,,,,
+7,"C",3,4,,,,,,,,,,,,
+7,"C",4,4,,,,,,,,,,,,
+7,"C",5,6,,,,,,,,,,,,
+8,"C",1,10,,,,,,,,,,,,
+8,"C",2,8,,,,,,,,,,,,
+8,"C",3,6,,,,,,,,,,,,
+8,"C",4,7,,,,,,,,,,,,
+8,"C",5,5,,,,,,,,,,,,
+9,"C",1,9,,,,,,,,,,,,
+9,"C",2,9,,,,,,,,,,,,
+9,"C",3,6,,,,,,,,,,,,
+9,"C",4,8,,,,,,,,,,,,
+9,"C",5,10,,,,,,,,,,,,
+10,"C",1,4,,,,,,,,,,,,
+10,"C",2,5,,,,,,,,,,,,
+10,"C",3,6,,,,,,,,,,,,
+10,"C",4,5,,,,,,,,,,,,
+10,"C",5,5,,,,,,,,,,,,
+11,"C",1,8,,,,,,,,,,,,
+11,"C",2,7,,,,,,,,,,,,
+11,"C",3,8,,,,,,,,,,,,
+11,"C",4,7,,,,,,,,,,,,
+11,"C",5,11,,,,,,,,,,,,
+12,"C",1,13,,,,,,,,,,,,
+12,"C",2,10,,,,,,,,,,,,
+12,"C",3,10,,,,,,,,,,,,
+12,"C",4,6,,,,,,,,,,,,
+12,"C",5,7,,,,,,,,,,,,
+13,"C",1,14,,,,,,,,,,,,
+13,"C",2,9,,,,,,,,,,,,
+13,"C",3,10,,,,,,,,,,,,
+13,"C",4,9,,,,,,,,,,,,
+13,"C",5,7,,,,,,,,,,,,
+14,"C",1,6,,,,,,,,,,,,
+14,"C",2,6,,,,,,,,,,,,
+14,"C",3,7,,,,,,,,,,,,
+14,"C",4,6,,,,,,,,,,,,
+14,"C",5,6,,,,,,,,,,,,
+15,"C",1,7,,,,,,,,,,,,
+15,"C",2,8,,,,,,,,,,,,
+15,"C",3,12,,,,,,,,,,,,
+15,"C",4,12,,,,,,,,,,,,
+15,"C",5,6,,,,,,,,,,,,
+16,"C",1,7,,,,,,,,,,,,
+16,"C",2,10,,,,,,,,,,,,
+16,"C",3,8,,,,,,,,,,,,
+16,"C",4,6,,,,,,,,,,,,
+16,"C",5,6,,,,,,,,,,,,
+17,"C",1,6,,,,,,,,,,,,
+17,"C",2,5,,,,,,,,,,,,
+17,"C",3,5,,,,,,,,,,,,
+17,"C",4,7,,,,,,,,,,,,
+17,"C",5,5,,,,,,,,,,,,
+18,"C",1,5,,,,,,,,,,,,
+18,"C",2,7,,,,,,,,,,,,
+18,"C",3,6,,,,,,,,,,,,
+18,"C",4,6,,,,,,,,,,,,
+18,"C",5,10,,,,,,,,,,,,
+19,"C",1,5,,,,,,,,,,,,
+19,"C",2,9,,,,,,,,,,,,
+19,"C",3,7,,,,,,,,,,,,
+19,"C",4,8,,,,,,,,,,,,
+19,"C",5,7,,,,,,,,,,,,
+20,"C",1,7,,,,,,,,,,,,
+20,"C",2,7,,,,,,,,,,,,
+20,"C",3,7,,,,,,,,,,,,
+20,"C",4,10,,,,,,,,,,,,
+20,"C",5,5,,,,,,,,,,,,
+1,"D",1,8,,,,,,,,,,,,
+1,"D",2,9,,,,,,,,,,,,
+1,"D",3,10,,,,,,,,,,,,
+1,"D",4,10,,,,,,,,,,,,
+1,"D",5,10,,,,,,,,,,,,
+2,"D",1,10,,,,,,,,,,,,
+2,"D",2,13,,,,,,,,,,,,
+2,"D",3,13,,,,,,,,,,,,
+2,"D",4,10,,,,,,,,,,,,
+2,"D",5,9,,,,,,,,,,,,
+3,"D",1,8,,,,,,,,,,,,
+3,"D",2,9,,,,,,,,,,,,
+3,"D",3,10,,,,,,,,,,,,
+3,"D",4,8,,,,,,,,,,,,
+3,"D",5,9,,,,,,,,,,,,
+4,"D",1,11,,,,,,,,,,,,
+4,"D",2,12,,,,,,,,,,,,
+4,"D",3,11,,,,,,,,,,,,
+4,"D",4,10,,,,,,,,,,,,
+4,"D",5,13,,,,,,,,,,,,
+5,"D",1,4,,,,,,,,,,,,
+5,"D",2,4,,,,,,,,,,,,
+5,"D",3,5,,,,,,,,,,,,
+5,"D",4,3,,,,,,,,,,,,
+5,"D",5,4,,,,,,,,,,,,
+6,"D",1,4,,,,,,,,,,,,
+6,"D",2,7,,,,,,,,,,,,
+6,"D",3,6,,,,,,,,,,,,
+6,"D",4,5,,,,,,,,,,,,
+6,"D",5,4,,,,,,,,,,,,
+7,"D",1,4,,,,,,,,,,,,
+7,"D",2,4,,,,,,,,,,,,
+7,"D",3,5,,,,,,,,,,,,
+7,"D",4,4,,,,,,,,,,,,
+7,"D",5,5,,,,,,,,,,,,
+8,"D",1,6,,,,,,,,,,,,
+8,"D",2,6,,,,,,,,,,,,
+8,"D",3,8,,,,,,,,,,,,
+8,"D",4,6,,,,,,,,,,,,
+8,"D",5,5,,,,,,,,,,,,
+9,"D",1,9,,,,,,,,,,,,
+9,"D",2,10,,,,,,,,,,,,
+9,"D",3,11,,,,,,,,,,,,
+9,"D",4,9,,,,,,,,,,,,
+9,"D",5,5,,,,,,,,,,,,
+10,"D",1,5,,,,,,,,,,,,
+10,"D",2,5,,,,,,,,,,,,
+10,"D",3,6,,,,,,,,,,,,
+10,"D",4,6,,,,,,,,,,,,
+10,"D",5,5,,,,,,,,,,,,
+11,"D",1,11,,,,,,,,,,,,
+11,"D",2,10,,,,,,,,,,,,
+11,"D",3,9,,,,,,,,,,,,
+11,"D",4,10,,,,,,,,,,,,
+11,"D",5,10,,,,,,,,,,,,
+12,"D",1,9,,,,,,,,,,,,
+12,"D",2,6,,,,,,,,,,,,
+12,"D",3,8,,,,,,,,,,,,
+12,"D",4,7,,,,,,,,,,,,
+12,"D",5,7,,,,,,,,,,,,
+13,"D",1,9,,,,,,,,,,,,
+13,"D",2,12,,,,,,,,,,,,
+13,"D",3,11,,,,,,,,,,,,
+13,"D",4,11,,,,,,,,,,,,
+13,"D",5,6,,,,,,,,,,,,
+14,"D",1,6,,,,,,,,,,,,
+14,"D",2,7,,,,,,,,,,,,
+14,"D",3,6,,,,,,,,,,,,
+14,"D",4,6,,,,,,,,,,,,
+14,"D",5,6,,,,,,,,,,,,
+15,"D",1,11,,,,,,,,,,,,
+15,"D",2,12,,,,,,,,,,,,
+15,"D",3,12,,,,,,,,,,,,
+15,"D",4,12,,,,,,,,,,,,
+15,"D",5,6,,,,,,,,,,,,
+16,"D",1,6,,,,,,,,,,,,
+16,"D",2,8,,,,,,,,,,,,
+16,"D",3,7,,,,,,,,,,,,
+16,"D",4,9,,,,,,,,,,,,
+16,"D",5,7,,,,,,,,,,,,
+17,"D",1,7,,,,,,,,,,,,
+17,"D",2,7,,,,,,,,,,,,
+17,"D",3,4,,,,,,,,,,,,
+17,"D",4,7,,,,,,,,,,,,
+17,"D",5,5,,,,,,,,,,,,
+18,"D",1,6,,,,,,,,,,,,
+18,"D",2,8,,,,,,,,,,,,
+18,"D",3,7,,,,,,,,,,,,
+18,"D",4,12,,,,,,,,,,,,
+18,"D",5,12,,,,,,,,,,,,
+19,"D",1,7,,,,,,,,,,,,
+19,"D",2,7,,,,,,,,,,,,
+19,"D",3,5,,,,,,,,,,,,
+19,"D",4,12,,,,,,,,,,,,
+19,"D",5,8,,,,,,,,,,,,
+20,"D",1,10,,,,,,,,,,,,
+20,"D",2,7,,,,,,,,,,,,
+20,"D",3,9,,,,,,,,,,,,
+20,"D",4,8,,,,,,,,,,,,
+20,"D",5,8,,,,,,,,,,,,
diff --git a/tests/fixtures/golden_masters/phased_xbar_s/S_data.parquet b/tests/fixtures/golden_masters/phased_xbar_s/S_data.parquet
new file mode 100644
index 0000000000000000000000000000000000000000..eb544fa2df114a5ede7ed81d8e1460630cf71153
GIT binary patch
literal 6089
zcmb_g4QyN06@GRd&uQZ|jrv@N1(sRG3L>QbNz#UOx_x%yzd0ds8~+VR-?N?A>3c86
zKXH;$)M>B{4VcE*U|S_tq*bMw5JOehO%*E)q;?g9QhwAS*!UTUpH)Z<CN_=joco>~
zJI{f#^(ndU-gE9f-#Ons_uMCsa}g(NV2_*F3;Hk{V3`hvIbhhfj$s&|bILj8b3x-A
zT+c8rY~0xNV>5uwU>`KD!DL-ya!vIQZf9I5W0Rmf%9z)0II8?J{pKX}@c(Gb`c9C-
zJa}Ef^!KSoI*$SfT|;M?*Q;x&1?`&7qg`%v(ls^6GHjR~?}Z@RStD!g8f1FbK@3&3
z6$jOTHf(?Z?f2-R4d{&Pc67H5^Id@ldc&`8AKyH>D|18qkaYc9`*!!-1v7x313ESu
zX4@z|V~*B(cKXZ^soBi9`YvNy%o{Fa+|b#mpn^+%tgFGL)|yLJ>?#gXcQ(Zpch&&;
zn}fFCv&W>(!DG@o#CNcLVF(I$$7)!^wCc9wTzm}5ivg(3Y|k>3)l}73P1U&;rl*ON
zf_Drn;QDp<2s@^8?YCQ65ci~kz1vXDi=XRH#fOfbnt$|n<k$!Qes1CI*T_F_nLSbZ
z;bHR4A9=qgzV#CM1f+K7DRQ6O@vW18I7`NacAeWqe)<!gPW%lSxbe5aHvJiLc{2Ug
zKVSD0d9wTALiPvGlLJraT1}4;UrVC*H2Lyv?@k=Q;$`yE$2xw;IRfUbhUhB<?HwF>
z#T+o4Bf499tcDlJ??3Fh{`-$VY2I=8JONbGFT0*4fmfF{KIDIlyazUJd-ql1+5OjN
zY%iT95BX8?A{o*>ZMy$YM0Z#3o%g;%Ufyg3&(D$8U$-42U)rMU=zX3HLm>Zp@H}~7
zM0b!mWqx@0SLE*)to^(i<r%Uv)2yZ{A68S95#(Ac{k(a%o4v7`c8;zKHzR{*NvrAB
z(*&9O+s!YMtAK(Jd!8dt|6}Ltzxdg6s;_6!$A;9inyM?TCdJTcf<9`cpP7v9?8Y#w
zr^#sQQqvFfuY<MJ=U-K0HL1P*h^LiOyw+|d|8xONyZw6hBVCy7rWpHm3tMmJU?1Nk
zXAX`4r|*$n5b6`%Z&Y)|Y(4n~8Nc&hNV>lenxXf<Lr%{>djDlxJ|J!StG{D8|0Zdt
zx%)oZbpDS()!Sz4t?v`c-Suw~%#n|8d51iY*=u@_D4F^ulDLiCULRD|s?}7@Qw&Ng
zr9pLxN>GM%?O)M*Cr9SzXYM=l;F+gRd~)GGw5<HEOfpTGWVB3n9vEP(#$qa)k$AS9
zVcV|MTU%xdbLIS%Oc?i0k>s;{F;kToV2(~^*tPgBDsN*dBbPD-Ww1ufi)ycoJzQWG
zGK+I^dP>Ymxl)m_1UkDq@31pnY+L65V{GqKQCYnuVEOi1%SnU%$7}4S&#KX#y^stH
ztH<nWhu3z?ea22>r=1Ca>kZav1l0QDddq&j{X~oXlmX4JR&TYYzuVGa+MqhFld0JY
zI?MYy``ur#pE7oRuG*>rx~e{awN;cNMLA-$>^IuqX|o?U0gASM;J$06u<sX#qHZIF
zI;^)G(c2$cZ$JJ;fYO!$gwjec8VZwy;+)V~o&e%D*k5Z09PROMn#%hx9tTp8;=HiN
z^2{3hiH-I{9e|@fk%*&}ekM*e4H^@?F0w-w`$ufE$W}k3t8bR&h084gyUFs%CfHG{
zG0)3szG!{4b+0War!xyvsX4JM$wk|>-QU=2+XFJ_*$TWold<)7+hktKfIu#bB8-Y9
zzEFapgZAh0&~*>eeg^x_eK){BF*rr7XH>h5W~Or)F+J5NQ~a7LE#?u=R+e+AxwPU!
zS;{Zg6m!CCCRIY$q)drV^Ccdp@3m!Qn4XhoaP8zS-!<E8`|6x&^uk7~6`E*pWLf-w
z`GQpK>jjt8_-wR}XwhooON%CJk*H_|N=_~f_TyB!f!&(8j$BPZMJ)sst$>&o=lG@|
zb&5}eQSJvV*0TH(144PxqN&153<VWN3#AG#F;v4Xss+`kSHRaY;vy3TD*#_A6>>tk
zL{I--I1IK_A;XvE3iPOoIibK87V)I!7x_YAu7Y08X9~sKoQzV}V4tJUVZ(K*Gxf?o
zcc0s~kFmCKk;rwiYn2ZX>m40i@Y%-sIgXove|)~qas8ZUhVyaU4CkGpGLNT0UqRtW
zO+P*XeH_C>tx$cC<0>4igEtx(9)&)RIWF7}b3I#wl2}gqqy2H$*t`&kmIQe;cU{f{
z+=dE#Y+z4tc5Wsphq8i{7K6FUOi&ug3b9F?Q%=VgieWF8OSvQ2lso}4uA+|TQ{LcU
zaCp?Ejzwkc8;b8l8!DDKlozCknDIxK5GUp?<^;)CO74WU{sE`puM~scki^HL#iVbj
zoQOF^oF~Z9qTmHPvIqPEZ0C?86ZcSEP>z=Q*w}p1Kbb>pm|IE(#OZWEtWdk&(0oG5
ziwXCHqK`!#iC83?_K#ppz(=&4a1HH&*ku~K7@dsd9Q{Gr+35Q+Pc^@SQhHHvM=J_X
z$$T=luqW=S^2Ozn?$CTX&QaWy-v_+=mg4TnKni>X?sDnaKpwHdPmCM*byhOby4a+l
z94|#@(_X4CCF7AS*l;9jd@T%xy%@tXX3C$<rd$daNl6q1e`LDGu{TuV!3Pzul7RhA
z#K-Enate|-E4U)ny>|+6&p^G;f-7`WGB)N=Vs?sYzYlCgvbB8*IEukQ-M&IxtMMb`
zXk5+w$P4m>YfOo|2`^$fJ)HPFJ}j-~LuzqVA0)v)lwHPw{9M6BQd*c##A|+D)czz~
z(Q?`!U0ltd27R-9my$)n<#UY9s`;$V_Z5C$)E*U6{tC_0COK2nHOt`@95mUh75+A!
zKU!0!p+-RcSw2^Cv?3^buc-xa-mAK1d9xz+OSETFbU<(oz}a5N3o+55=8V=iDI{iM
zLuII4MLy=trvsy_o_(k{pSFkf`QJ3Jl)CD3qywRRLXJ9;aazZ>B9~Ba7n3ocLvRPx
z8n*y*LQu2!(7FTl5&BS@l-zeE<Dn%U{HFD9Hs}h^a^+EP2x<$|`dZyeiYa+4FSz=t
zjnQ2HQh0VWHw||W+*=$MP-`<>LteaZ{1v=*DtI?tjpI?A<V|2#z(%>bXuB}Jb;}mK
zz9xe7CWFfh2eidz7`uKAx2P)ORf-*M5}V;|TjuH5!!6)gk<$B*4!L1`rQ@C=7Zm6v
zj8|o79J&;Ll!C$~49_w?L;?5<6@l}XRXi^@vZSsb8Kyx{64a>ZdJNSwn}mH%#v!7y
zLJcpo0?$w9*WqC|Fp9jw!i2&v&YAJ<QSm)omM%(RTNfYg3P?5F{}3#!PbvF>BK*97
z$5q3hq@XMDJrv(F2X}*r^A~YG@CmCGe=6K4$|3Hm(;iO*B>)Cboj(*%y$4Z6m(=OM
t##p=aeBq{YrgTjzpTCBFWOb__QQgYtRp%Y>TTRnH<_-8IdmMhm|1b2AxuyUB

literal 0
HcmV?d00001

diff --git a/tests/fixtures/golden_masters/phased_xbar_s/S_metadata.json b/tests/fixtures/golden_masters/phased_xbar_s/S_metadata.json
new file mode 100644
index 0000000..06793b0
--- /dev/null
+++ b/tests/fixtures/golden_masters/phased_xbar_s/S_metadata.json
@@ -0,0 +1,53 @@
+{
+  "center_col": "center",
+  "chart_type": "S",
+  "group_col": "group",
+  "insufficient_phases": null,
+  "lane_boundaries": [
+    {
+      "label": "F1_1_F2_2",
+      "position": 6,
+      "variables": [
+        "phase"
+      ]
+    },
+    {
+      "label": "F1_2_F2_1",
+      "position": 12,
+      "variables": [
+        "phase"
+      ]
+    },
+    {
+      "label": "F1_2_F2_2",
+      "position": 18,
+      "variables": [
+        "phase"
+      ]
+    },
+    {
+      "label": "F1_3_F2_1",
+      "position": 24,
+      "variables": [
+        "phase"
+      ]
+    },
+    {
+      "label": "F1_3_F2_2",
+      "position": 30,
+      "variables": [
+        "phase"
+      ]
+    }
+  ],
+  "n_mode": "actual",
+  "n_sigma": 3.0,
+  "phase_by": [
+    "factor 1",
+    "factor 2"
+  ],
+  "phased": true,
+  "run_rules_applicable": false,
+  "single_point_phases": 0,
+  "value_col": "s"
+}
diff --git a/tests/fixtures/golden_masters/phased_xbar_s/S_statistics.json b/tests/fixtures/golden_masters/phased_xbar_s/S_statistics.json
new file mode 100644
index 0000000..6a89b33
--- /dev/null
+++ b/tests/fixtures/golden_masters/phased_xbar_s/S_statistics.json
@@ -0,0 +1,7 @@
+{
+  "N": null,
+  "center": null,
+  "limits_vary": true,
+  "lpl": null,
+  "upl": null
+}
diff --git a/tests/fixtures/golden_masters/phased_xbar_s/Xbar_data.parquet b/tests/fixtures/golden_masters/phased_xbar_s/Xbar_data.parquet
new file mode 100644
index 0000000000000000000000000000000000000000..35f76be2a265fd5f3a2c88d49d95b3c6d5781257
GIT binary patch
literal 6294
zcmb_h3v64}8NPO$%xMyucD3hJ84IRZof;B9lC*)s&b1TYN0T&h6TeoL`r3}|^xhlS
zk2q=R)(yJR@zMqZv{mA<hZsVLNnO=V73vsc?7;-pUOXnoD#XyXs?#(fGzR<sb8W}=
zO(1B!s?R;=|Ihiq@Bhz#{&SU1vH_RHWI1EDTr`eZcnh<cVNRNMwJ{9icg?zH{cdPn
zLmL^!jg1GJ0c-}b8R~<^Jru2J%<kENq1}udWtJ!?Pca=EH=WXc%s@vJdiXolys-<U
zFb`f=G6Q{OBVDHegs!2pIyRQqPz(AsU8lM|NYXt!WMM30mdRdlqSMl1Y3Uwfwzq)~
z6}1fq6^Ax$0tX%U8lmM4Ep4QyeI(Y+KhisPVfW;=@#*AE;qNK?-#mI?`&}>t;&akq
z365FXX?P}_t+mnV>i|zXIv97~R;IaQ(^kd<oh=$FQ0lX|>nOEWDcP{A5u$WH5!UE5
zLGZT;cH#VKWt(tXX+wOcrEd(J!rif%nwVC@Zk&scNqaE>wav191<GcwXl&-nTr0DE
z6)6SZJfeXc+wPUfT&DAcwHa~GnJkZ+DtYmc@lyEir~kcu|APV<{M?U)cB4pM`pw(>
zE;Nhe@$P8i`hJl-|7&<;fxKvc<(<!eRV25Z`{l!<&-%$@|2#h~IeD_M86-d8$e|lH
z{NdZ@II{ef!)5N}$cxR+9-h2@c*_%<fp0xv&hf-BXWqf_t-BJh^CX}C+COgI%abQ_
z(CHP(Pl3hVe;IA<ycBEah|3J{(|+>kvwL=)nDUdW0ek-4QL^&{d*F>uj<g$he8u{f
zKs@-wMe>FDvxR$>L~{LWd+)seN|AgLxW5MEaoO1|dr0vfn<;n~M?Q7@Z&PO`_*Prm
z#FIa5Yx&|guI9-#@bbOF5p?9?W`Q)Ooy}Y^V>6d&5)8av`caecfJLh0AOs6z^|nA5
zgW<L(If5S9ZsoN&XZ>1S_lsl>9|w;t>oU=1rjdskHZ$EDSc5F}(htqeV@As{s2x2t
zJC7N<e-4|#J}Ho9^CjOy0y*%;3-(ujB#@s#BAmX3BR}hbO+3VDTX{1_E`0Y}&-6fM
z;2!Ml;0f+!VuzpLp1pL1pG-Gj^&LMs{jdN2cIh2K%O|Ty{yuf)%HwS!na1TWaOBSW
z)y=p3nL|H!S$5YR1Db)b0**o1BKoP9Qa`(8?HHKaI<GSNW=0nmQx81${a-zQ_PvYm
z(?S5weyPUUOE;o-@7ySoUo1ZL;MTn&>Bl@76v=)(DX(%QJbdc8#fu!7M1~hQ^1CI-
z_&@T5o}d5lt<Jyx(661Q^L{NB=loh?J|$|W<zbOv%H96BNbd99F28w%C!fLNC(~VR
z@Rp`DthAXa5wMNk>ZR-R)@kXjmL2B(4#Vzpv0)5rR-KB*W`4i?AOXFU5?HA^)3+U;
zK5=|nI1W{(>*OF~YstscNhN0KWGwC17;Vj|+(I#X4Kvma$eJXZj^&dTi3#QyRK~Ia
z-liofmeaZ|@uXTvmIqC8R@TH#MeO044k@|3peAPJjFKtj87tq_-F1h9>9(|Y4KgjA
zU1ea4(aKxDyuo_T<hZZNVg9hYrK=ZeCu1=(&~(CUxAlRRu9hwb!vp&|>(mKq{b{51
zn9*^z+3}nS&98M&=uO|&TxZ%;CT>uv*^372C4=MZA9Fm{(*2R{1Rc<)YaQ57Mxn8$
zK##RpzuDrrwZn173_$v7hiRafLK>_OsJxRD>QSTh38UkMjgIqI0F=IDA(URaZYazX
ziu3y>>&s1!hd%B&wh3_bB?EEv()Hs&7FKXxZF0P8b-cU<aP)^2arDv;#i?XMeU4K^
z$cl>w$9p!%X&WHwPbVVkr5}iBeWb_AJIvOTTVY3SE!miwh~;h1v>vf%)I@S=HohPi
z6*X_auIIBy><2*xJ$o*urjqvF9=jS-k|0ovvJ9hqA(kt^&`J9XS?GF(Xg`U4*U_6`
zpc$N{){|ws^=9TWNjWiFFVlF<7M8P!XD_Om_(DRXP*k$ZRmFldpNtofnvyKU60t%I
zrXR5<RhXVpQ@D0!+JEgX`_URQ^(^ePW}%5XA}iwmiwlP8XfIGI`)8-6M2}V#UtTm-
zk3>ZcC>gadG=Njpy69HLHNn+6gk@ZhQjs<Q&B+U~RiSDW9|Wdd6neN7`DF%#BBMtu
z3okPi)EPZgS$LVD>aJ5gs7~DgU(JhkrVG}ASfP;1NW}uJ0Y~6S*yFimtgw)y2UgBV
zxma!)&uex$mdh=akkw)`m(MJyD0L6@Is2S;TvsMjtL*didF)3STRR&F><wL~eU#W9
z*xZlLt`x_z>=gX)*~78}tT)B_SvJM`QdH*k*6B+q9IfgnOVGzLJk$!+3oKh=VI92D
z$jA)zam=#gQ!v-NQ&8k$lnV}o-4lxvA1p}fcxG?L3$Yb*vC!Z_VSXVcsA5`D60(pf
zr37U#Ern)qPB9T$%8&WjOxzPl$JHs2v8D2OG42zFgpqM~c`R#V|8RI8+9+em!&ylQ
z$Vo1E7;!?Ld`43Ih3Gz5%MH3Du9O#iq7n-Q^HKkBF%ojgI8Rc8dC3QMR4?!XY}c?e
z8TL}0pazSv(8OYto5>(H%q_%uc`m`rC2H3vE=H8B9PvzP`cTjr2?f##ZWMijcm#_P
z_wYgRU8TOu!I?nDIUuO6dhRQ6s`3()#IocGmb5rUv(eDf!LYj$FE$hPh>M9ZOL60z
zAL8ym9QFhT<G>YSmq~;Mvxp5m(Qk;atCS4Z_@)eJVoGp6;iLLOG#p5S4QHepuchHJ
zANsI@8RydJxLb=wRFP$g3(Qqx>=R2d;84aZMPR=p;fdN<xg<rNm)wEM-n*o*cd(YT
z<Q5M_LlaKTXP2Dd{9q%HuI>}>%nN+YzJgzC#Yf6ezgFXsmDDNsgy#1uagmFOk;q5k
zFc4Av%WLsSa$-8ZHWx+WhSMu}kf#kXi7HErk#Lpoy7npJ4i*z!a7|8i`qlDYiRLA@
z-#IZ~&TD<%H}GE9UghIliRSAnxl`4xmdgz>SY^)|_|~68T3hCzRzMwEIbUk9BxyOI
zs|j%KE4tP4r@{Bjv}a0iP;w8#8DGjuA=w#UrunPqrigMfG+cxlmXC#8*#tknW(|P4
z^Fe!9dk$93E3Llzoe5sdM%17)8m9GpC&m)$?s7EbcS;_iT<ex#4&=Iekk%illhB76
zrTO8GhQ-4%;7#k`yx<<2XN%)L5o!$7{A&G*%5ilfE4c@#jq%LD;j#Jg%$#<U;hn}-
z=0X?lAsFK|H8KU)(8w+vFJmL({p8JKS6hnh2;0Bd-;WA5jT(B1Y4;Z$vJ?1<H0}dl
zeZG0TR8j4v_9o|dXyefldP!<h?Y_iW>^M$!LB6vo!0bmGDK>>yGb-36Xjyg!8+w0Y
z2MNL>@xXgeu8b$J`zFilr^YEX4MO#5hxO<xHz7-EK#oJ-RmJl<8}Q<2zXp%;rS%Ec
zfgOlvM6-tzQ@(@M^;*0d*Hcpx;)8QQ9@fWM9}ejQpBwPFtNdqH=Z{j}y8w5EjPsXq
zKEwxBYkW#Hn7GSyJxuX>12_}F;HmM40+#PZRM91x{)QVETTeEYJ5)>-u8n82*U}HJ
dp7IA(kM^n6bqD+nXVrhqyYOfEN%$X${{i9k5_SLp

literal 0
HcmV?d00001

diff --git a/tests/fixtures/golden_masters/phased_xbar_s/Xbar_metadata.json b/tests/fixtures/golden_masters/phased_xbar_s/Xbar_metadata.json
new file mode 100644
index 0000000..45e976a
--- /dev/null
+++ b/tests/fixtures/golden_masters/phased_xbar_s/Xbar_metadata.json
@@ -0,0 +1,58 @@
+{
+  "center_col": "center",
+  "chart_type": "Xbar",
+  "group_col": "group",
+  "insufficient_phases": null,
+  "lane_boundaries": [
+    {
+      "label": "F1_1_F2_2",
+      "position": 6,
+      "variables": [
+        "factor 1",
+        "factor 2"
+      ]
+    },
+    {
+      "label": "F1_2_F2_1",
+      "position": 12,
+      "variables": [
+        "factor 1",
+        "factor 2"
+      ]
+    },
+    {
+      "label": "F1_2_F2_2",
+      "position": 18,
+      "variables": [
+        "factor 1",
+        "factor 2"
+      ]
+    },
+    {
+      "label": "F1_3_F2_1",
+      "position": 24,
+      "variables": [
+        "factor 1",
+        "factor 2"
+      ]
+    },
+    {
+      "label": "F1_3_F2_2",
+      "position": 30,
+      "variables": [
+        "factor 1",
+        "factor 2"
+      ]
+    }
+  ],
+  "n_mode": "actual",
+  "n_sigma": 3.0,
+  "phase_by": [
+    "factor 1",
+    "factor 2"
+  ],
+  "phased": true,
+  "run_rules_applicable": false,
+  "single_point_phases": 0,
+  "value_col": "xbar"
+}
diff --git a/tests/fixtures/golden_masters/phased_xbar_s/Xbar_statistics.json b/tests/fixtures/golden_masters/phased_xbar_s/Xbar_statistics.json
new file mode 100644
index 0000000..6a89b33
--- /dev/null
+++ b/tests/fixtures/golden_masters/phased_xbar_s/Xbar_statistics.json
@@ -0,0 +1,7 @@
+{
+  "N": null,
+  "center": null,
+  "limits_vary": true,
+  "lpl": null,
+  "upl": null
+}
diff --git a/tests/test_analysis_golden_master.py b/tests/test_analysis_golden_master.py
index 607fc48..923c567 100644
--- a/tests/test_analysis_golden_master.py
+++ b/tests/test_analysis_golden_master.py
@@ -274,6 +274,11 @@ def _make_residual_chart():
     return synthetic.make_design(1, K1=3, K2=2, T=4, n_min=3, n_max=3, seed=42)
 
 
+def _make_phased_xbar_s():
+    """Scenario: phased Xbar/S — a centre line and limits per factor combination (#144)."""
+    return synthetic.make_design(1, K1=3, K2=2, T=6, n_min=2, n_max=4, seed=42)
+
+
 def _make_single_obs_strata():
     """Scenario 7: Stratum with n=1 — edge case for mR computation."""
     # SDS 3 has mixed replication. Use specific params to get some n=1 cells.
@@ -356,6 +361,14 @@ def _execute_single_obs_strata():
     return study.execute(chart='X', by=[], companion=True)
 
 
+def _execute_phased_xbar_s():
+    """Execute scenario: phased Xbar/S companion pair (#144)."""
+    df = _make_phased_xbar_s()
+    pb = ProcessBehavior(df)
+    study = pb.formulate(response='y', factors=['factor 1', 'factor 2'], time='time')
+    return study.execute(chart='Xbar', phased=True, companion=True)
+
+
 def _execute_pathological_ordering():
     """Execute scenario 8: reverse-sorted input, X companion."""
     df = _make_pathological_ordering()
@@ -404,6 +417,10 @@ def test_single_obs_strata(self):
         result = _execute_single_obs_strata()
         _assert_result_matches_snapshot('single_obs_strata', result)
 
+    def test_phased_xbar_s(self):
+        result = _execute_phased_xbar_s()
+        _assert_result_matches_snapshot('phased_xbar_s', result)
+
     def test_pathological_ordering(self):
         result = _execute_pathological_ordering()
         _assert_result_matches_snapshot('pathological_ordering', result)
diff --git a/tests/test_phased_limits.py b/tests/test_phased_limits.py
index d4c5cf6..82b435d 100644
--- a/tests/test_phased_limits.py
+++ b/tests/test_phased_limits.py
@@ -182,10 +182,12 @@ def test_phased_lane_boundaries_preserved(self, sds1_study):
 class TestPhasedValidation:
     """Test that invalid phased configurations raise appropriate errors."""
 
-    def test_phased_requires_xmr_or_r(self, sds1_study):
-        """phased=True with Xbar should raise ValidationError."""
-        with pytest.raises(ValidationError, match='only valid for X or mR'):
-            sds1_study.execute(chart='Xbar', phased=True)
+    def test_phased_requires_a_limits_chart(self, sds1_study):
+        """phased=True is refused on Histogram and accepted on Xbar/S (#144) as on X/mR."""
+        with pytest.raises(ValidationError, match='only valid for'):
+            sds1_study.execute(chart='Histogram', phased=True)
+        result = sds1_study.execute(chart='Xbar', phased=True)
+        assert result.charts['Xbar']['metadata']['phased'] is True
 
     def test_phased_stratified_allowed_with_collapsed_factors(self, sds1_study):
         """phased=True with by=['factor 1'] should work (factor 2 is collapsed)."""
diff --git a/tests/test_phased_xbar_s.py b/tests/test_phased_xbar_s.py
new file mode 100644
index 0000000..c4b3ebc
--- /dev/null
+++ b/tests/test_phased_xbar_s.py
@@ -0,0 +1,536 @@
+"""Phased (staged) Xbar/S charts: one centre line and one set of limits per phase (issue #144).
+
+A phase is one combination of the grouping factors, charted the way the per-condition charts
+of ``by=[time]`` chart one stratum, and the phases sit side by side on one axis. Dr. Bishop
+settled the two design questions on 2026-10-18: each stage has its own S-bar and sigma-hat
+("treat each stage separately as if it were a single analysis"), and with other factors pooled
+into the subgroup the centre is "the average of the rational subgroup averages for that stage".
+He also sent Wheeler's four-cavity thickness data with Minitab's per-cavity limits; that file is
+``tests/fixtures/data/wheeler_cavity_thickness.csv`` (Hour 1-20, Cavity A-D, Cycle 1-5,
+Thickness, plus VAS's per-cavity CL/LBL/UBL columns on its first four rows).
+"""
+
+from __future__ import annotations
+
+import time as _time
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import pytest
+
+import processbehavior as pb
+from processbehavior.calibration import Calibration
+from processbehavior.datasets.synthetic import make_design
+from processbehavior.exceptions import CalibrationNotSupportedError, ValidationError
+from processbehavior.spc_constants import b3, b4, c4
+
+REFERENCE_CSV = Path(__file__).resolve().parent.parent / 'validation' / 'PBTESTDATABASE_T100.csv'
+CAVITY_CSV = Path(__file__).resolve().parent / 'fixtures' / 'data' / 'wheeler_cavity_thickness.csv'
+FACTORS = ['FACTOR 1', 'FACTOR 2']
+TIME = 'PRODUCTION TIME'
+XBAR_COLS = ['n', 'xbar', 'center', 'lpl', 'upl', 'beyond_limits']
+S_COLS = ['n', 's', 'center', 'lpl', 'upl', 'beyond_limits']
+VALUE_COLS = {'Xbar': XBAR_COLS, 'S': S_COLS}
+
+
+@pytest.fixture(scope='module')
+def reference_df():
+    return pd.read_csv(REFERENCE_CSV, na_values=['*'])
+
+
+@pytest.fixture(scope='module')
+def reference_studies(reference_df):
+    return {
+        response: pb.formulate(reference_df, response=response, factors=FACTORS, time=TIME)
+        for response in ('PM SDS 1', 'PM SDS 3')
+    }
+
+
+@pytest.fixture(scope='module')
+def cavity_df():
+    return pd.read_csv(CAVITY_CSV)
+
+
+@pytest.fixture(scope='module')
+def cavity_study(cavity_df):
+    return pb.formulate(cavity_df, response='Thickness', time='Hour', factors=['Cavity'])
+
+
+@pytest.fixture
+def sds1_study():
+    df = make_design(1, K1=3, K2=2, T=6, n_min=2, n_max=4, seed=42)
+    return pb.formulate(df, response='y', time='time', factors=['factor 1', 'factor 2'])
+
+
+def _staged(values_by_stage: dict[str, list[list[float]]]) -> pd.DataFrame:
+    """A staged study: ``{stage: [readings at time 1, readings at time 2, ...]}``."""
+    rows = []
+    for stage, subgroups in values_by_stage.items():
+        for t, readings in enumerate(subgroups, start=1):
+            rows += [{'stage': stage, 'time': t, 'y': float(v)} for v in readings]
+    return pd.DataFrame(rows)
+
+
+def _focus_concat(result, chart: str) -> pd.DataFrame:
+    """The per-condition charts of a stratified result, one after another, in strata order."""
+    return pd.concat([result.focus(k).get_chart(chart) for k in result.strata], ignore_index=True)
+
+
+def _frame(result, chart: str) -> pd.DataFrame:
+    return result.get_chart(chart)[VALUE_COLS[chart]].reset_index(drop=True)
+
+
+# =============================================================================
+# Oracles: the phased chart is the per-condition charts laid side by side
+# =============================================================================
+
+
+class TestOracles:
+    @pytest.mark.parametrize('chart', ['Xbar', 'S'])
+    @pytest.mark.parametrize('response', ['PM SDS 1', 'PM SDS 3'])
+    def test_phased_lanes_equal_stratified_focus(self, reference_studies, response, chart):
+        study = reference_studies[response]
+        phased = study.execute(chart=chart, phased=True)
+        stratified = study.execute(chart=chart, by=[TIME])
+        expected = _focus_concat(stratified, chart)
+
+        pd.testing.assert_frame_equal(_frame(phased, chart), expected[VALUE_COLS[chart]].reset_index(drop=True))
+
+        data = phased.get_chart(chart)
+        phases = [k for k in stratified.strata for _ in range(len(stratified.focus(k).get_chart(chart)))]
+        assert list(data['phase']) == phases
+        assert list(data['group']) == [f'{k}_{t}' for k, t in zip(phases, expected[TIME], strict=True)]
+
+    @pytest.mark.parametrize('chart', ['Xbar', 'S'])
+    @pytest.mark.parametrize('response', ['PM SDS 1', 'PM SDS 3'])
+    def test_phase_statistics_equal_stratified_statistics(self, reference_studies, response, chart):
+        study = reference_studies[response]
+        phases = study.execute(chart=chart, phased=True).phase_statistics(chart).set_index('phase')
+        stratified = study.execute(chart=chart, by=[TIME])
+        assert list(phases.index) == stratified.strata
+        for stratum in stratified.strata:
+            stats = stratified.statistics(chart, stratum=stratum)
+            row = phases.loc[stratum]
+            for key in ('center', 'N', 'lpl', 'upl'):
+                if stats[key] is None:
+                    assert np.isnan(row[key]), (stratum, key, row[key])
+                else:
+                    assert row[key] == pytest.approx(stats[key]), (stratum, key)
+
+    @pytest.mark.parametrize('chart', ['Xbar', 'S'])
+    def test_pooled_phase_equals_reformulated_study(self, reference_df, chart):
+        """by=['FACTOR 1', time] with FACTOR 2 pooled: each phase is the chart of a study without FACTOR 2.
+
+        Checked on the balanced reference data and on an unbalanced copy, where the mean of the
+        pooled subgroup means (Bishop's rule) differs from the mean of the fine cell means.
+        """
+        df = reference_df
+        thinned = df.drop(
+            df[(df['FACTOR 1'] == 1) & (df['FACTOR 2'] == 1) & (df[TIME] <= 50)].groupby(TIME).head(2).index
+        )
+        for frame in (df, thinned):
+            two_factor = pb.formulate(frame, response='PM SDS 1', factors=FACTORS, time=TIME)
+            phased = two_factor.execute(chart=chart, by=['FACTOR 1', TIME], phased=True)
+            one_factor = pb.formulate(frame, response='PM SDS 1', factors=['FACTOR 1'], time=TIME)
+            expected = _focus_concat(one_factor.execute(chart=chart, by=[TIME]), chart)
+            pd.testing.assert_frame_equal(_frame(phased, chart), expected[VALUE_COLS[chart]].reset_index(drop=True))
+
+        # The non-phased pooled chart keeps its study-wide centre: the mean of the fine cell means.
+        two_factor = pb.formulate(thinned, response='PM SDS 1', factors=FACTORS, time=TIME)
+        pooled = two_factor.execute(chart='Xbar', by=['FACTOR 1', TIME]).get_statistics('Xbar')['center']
+        fine = thinned.groupby(FACTORS + [TIME])['PM SDS 1'].mean().mean()
+        plotted = thinned.groupby(['FACTOR 1', TIME])['PM SDS 1'].mean().mean()
+        assert pooled == pytest.approx(round(fine, 3))
+        assert abs(fine - plotted) > 1e-3  # the two rules really differ here
+
+    def test_by_none_equals_by_empty_and_the_full_by_in_any_order(self, reference_studies):
+        study = reference_studies['PM SDS 3']
+        base = study.execute(chart='Xbar', phased=True, companion=True)
+        for by in ([], FACTORS + [TIME], ['FACTOR 2', 'FACTOR 1', TIME]):
+            other = study.execute(chart='Xbar', by=by, phased=True, companion=True)
+            for chart in ('Xbar', 'S'):
+                pd.testing.assert_frame_equal(_frame(base, chart), _frame(other, chart))
+                assert list(base.get_chart(chart)['phase']) == list(other.get_chart(chart)['phase'])
+                pd.testing.assert_frame_equal(base.phase_statistics(chart), other.phase_statistics(chart))
+                assert base.charts[chart]['metadata']['phase_by'] == other.charts[chart]['metadata']['phase_by']
+
+    def test_wheeler_cavity_oracle(self, cavity_df, cavity_study):
+        """Tom's file carries VAS's per-cavity limits; the phased chart reproduces them to the digit."""
+        result = cavity_study.execute(chart='Xbar', phased=True, companion=True)
+        chart_cols = [
+            'XBAR CHART CL',
+            'XBAR CHART LBL',
+            'XBAR CHART UBL',
+            'S CHART CL',
+            'S CHART LBL',
+            'S CHART UBL',
+            'RSG SAMPLE SIZE',
+        ]
+        charts = cavity_df.dropna(subset=['CAVITY ID'])
+        assert (charts.groupby('CAVITY ID')[chart_cols].nunique() == 1).all().all()  # one value per cavity
+        vas = charts.groupby('CAVITY ID')[chart_cols].first()
+        xbar = result.phase_statistics('Xbar').set_index('phase')
+        s = result.phase_statistics('S').set_index('phase')
+        assert list(xbar.index) == ['A', 'B', 'C', 'D']
+        for cavity in 'ABCD':
+            assert xbar.loc[cavity, 'center'] == pytest.approx(vas.loc[cavity, 'XBAR CHART CL'], abs=5e-4)
+            assert xbar.loc[cavity, 'lpl'] == pytest.approx(vas.loc[cavity, 'XBAR CHART LBL'], abs=5e-4)
+            assert xbar.loc[cavity, 'upl'] == pytest.approx(vas.loc[cavity, 'XBAR CHART UBL'], abs=5e-4)
+            assert s.loc[cavity, 'center'] == pytest.approx(vas.loc[cavity, 'S CHART CL'], abs=5e-4)
+            assert s.loc[cavity, 'lpl'] == pytest.approx(vas.loc[cavity, 'S CHART LBL'], abs=5e-4)
+            assert s.loc[cavity, 'upl'] == pytest.approx(vas.loc[cavity, 'S CHART UBL'], abs=5e-4)
+            assert xbar.loc[cavity, 'N'] == vas.loc[cavity, 'RSG SAMPLE SIZE']
+        assert xbar.loc['A', 'center'] == 14.65 and xbar.loc['D', 'upl'] == 9.831  # the Minitab screenshot
+
+        stratified = cavity_study.execute(chart='Xbar', by=['Hour'], companion=True)
+        for chart in ('Xbar', 'S'):
+            expected = _focus_concat(stratified, chart)[VALUE_COLS[chart]].reset_index(drop=True)
+            pd.testing.assert_frame_equal(_frame(result, chart), expected)
+
+
+# =============================================================================
+# Hand-computed
+# =============================================================================
+
+
+class TestHandComputed:
+    def test_two_stage_example(self):
+        stages = {
+            'A': [[10, 12, 14], [11, 13, 12], [9, 10, 14], [12, 12, 15]],
+            'B': [[20, 22, 21], [19, 23, 22], [21, 21, 24], [20, 20, 20]],
+        }
+        study = pb.formulate(_staged(stages), response='y', time='time', factors=['stage'], precision=9)
+        result = study.execute(chart='Xbar', phased=True, companion=True)
+        xbar, s = result.get_chart('Xbar'), result.get_chart('S')
+        phases = result.phase_statistics('Xbar').set_index('phase')
+
+        for stage, subgroups in stages.items():
+            means = np.array([np.mean(g) for g in subgroups])
+            stds = np.array([np.std(g, ddof=1) for g in subgroups])
+            centre, s_bar = means.mean(), stds.mean()
+            sigma_hat = s_bar / c4(3)
+            half = 3 * sigma_hat / np.sqrt(3)
+            rows = xbar[xbar['phase'] == stage]
+            assert rows['xbar'].to_numpy() == pytest.approx(means, abs=1e-9)
+            assert rows['center'].to_numpy() == pytest.approx(centre, abs=1e-9)
+            assert rows['lpl'].to_numpy() == pytest.approx(centre - half, abs=1e-9)
+            assert rows['upl'].to_numpy() == pytest.approx(centre + half, abs=1e-9)
+            srows = s[s['phase'] == stage]
+            assert srows['s'].to_numpy() == pytest.approx(stds, abs=1e-9)
+            assert srows['center'].to_numpy() == pytest.approx(s_bar, abs=1e-9)
+            assert srows['lpl'].to_numpy() == pytest.approx(b3(3) * s_bar, abs=1e-9)  # 0 at n = 3
+            assert srows['upl'].to_numpy() == pytest.approx(b4(3) * s_bar, abs=1e-9)
+            row = phases.loc[stage]
+            assert row['subgroups'] == 4 and row['N'] == 3
+            assert (row['center'], row['s_bar'], row['sigma_hat']) == pytest.approx(
+                (centre, s_bar, sigma_hat), abs=1e-9
+            )
+            assert (row['lpl'], row['upl']) == pytest.approx((centre - half, centre + half), abs=1e-9)
+
+    def test_unequal_n_uses_eq_11_5_within_the_phase(self):
+        stages = {
+            'A': [[10, 12, 14], [11, 13, 12, 16], [9, 10, 14, 11], [12, 12, 15]],
+            'B': [[20, 22, 21], [19, 23, 22], [21, 21, 24], [20, 20, 25]],
+        }
+        study = pb.formulate(_staged(stages), response='y', time='time', factors=['stage'], precision=9)
+        result = study.execute(chart='Xbar', phased=True)
+        xbar = result.get_chart('Xbar')
+        phases = result.phase_statistics('Xbar').set_index('phase')
+
+        a = stages['A']
+        stds = np.array([np.std(g, ddof=1) for g in a])
+        sizes = np.array([len(g) for g in a])
+        sigma_hat = np.mean(stds / np.array([c4(n) for n in sizes]))  # Eq 11-5 within the phase
+        rows = xbar[xbar['phase'] == 'A']
+        assert (rows['upl'] - rows['center']).to_numpy() == pytest.approx(3 * sigma_hat / np.sqrt(sizes), abs=1e-9)
+        assert phases.loc['A', 'sigma_hat'] == pytest.approx(sigma_hat, abs=1e-9)
+        assert all(np.isnan(phases.loc['A', k]) for k in ('N', 'lpl', 'upl'))  # the limits step within A
+        assert phases.loc['B', 'N'] == 3 and not np.isnan(phases.loc['B', 'upl'])  # B has one pair
+        assert result.get_statistics('Xbar') == {
+            'N': None,
+            'center': None,
+            'lpl': None,
+            'upl': None,
+            'limits_vary': True,
+        }
+
+    def test_average_n_mode_is_applied_per_phase(self):
+        stages = {
+            'A': [[10, 12, 14], [11, 13, 12, 16], [9, 10, 14, 11], [12, 12, 15]],
+            'B': [[20, 22], [19, 23, 22, 25], [21, 21, 24, 26], [20, 20]],
+        }
+        study = pb.formulate(_staged(stages), response='y', time='time', factors=['stage'], precision=9)
+        result = study.execute(chart='Xbar', phased=True, companion=True, n_mode='average')
+        xbar = result.get_chart('Xbar')
+        phases = result.phase_statistics('Xbar').set_index('phase')
+        for stage, subgroups in stages.items():
+            n_avg = np.mean([len(g) for g in subgroups])
+            s_bar = np.mean([np.std(g, ddof=1) for g in subgroups])
+            half = 3 * (s_bar / c4(int(n_avg))) / np.sqrt(n_avg)  # Eq 11-20, as the unphased chart applies it
+            rows = xbar[xbar['phase'] == stage]
+            assert phases.loc[stage, 'N'] == pytest.approx(n_avg)
+            assert rows['upl'].nunique() == 1
+            assert (rows['upl'] - rows['center']).iloc[0] == pytest.approx(half, abs=1e-9)
+        assert 'n_avg' not in result.charts['Xbar']['metadata']
+        assert 'n_avg' not in result.charts['S']['metadata']
+
+    def test_n_sigma_two_halves_the_width_ratio(self, sds1_study):
+        three = sds1_study.execute(chart='Xbar', phased=True).get_chart('Xbar')
+        two = sds1_study.execute(chart='Xbar', phased=True, n_sigma=2).get_chart('Xbar')
+        assert ((two['upl'] - two['center']) / (three['upl'] - three['center'])).to_numpy() == pytest.approx(
+            2 / 3, abs=1e-2
+        )
+
+
+# =============================================================================
+# Edge cases
+# =============================================================================
+
+
+class TestEdgeCases:
+    def test_one_reading_subgroups_sit_on_xbar_and_not_on_s(self, reference_studies):
+        result = reference_studies['PM SDS 3'].execute(chart='Xbar', phased=True, companion=True)
+        xbar, s = result.get_chart('Xbar'), result.get_chart('S')
+        singles = xbar[xbar['n'] == 1]
+        assert len(singles) > 0
+        assert (s['n'] >= 2).all()
+        assert len(s) == len(xbar) - len(singles)
+        phases = result.phase_statistics('Xbar').set_index('phase')
+        for _, row in singles.iterrows():
+            # N_r = 1: center ± n_sigma · sigma_hat (10-12 manual Eqs 11-16, 11-17)
+            assert row['upl'] - row['center'] == pytest.approx(3 * phases.loc[row['phase'], 'sigma_hat'], abs=3e-3)
+
+    def test_single_subgroup_phase_has_finite_limits(self):
+        stages = {'A': [[10, 12, 14], [11, 13, 12], [9, 10, 14]], 'B': [[20, 22, 21]]}
+        result = pb.formulate(_staged(stages), response='y', time='time', factors=['stage']).execute(
+            chart='Xbar', phased=True, companion=True
+        )
+        b = result.get_chart('Xbar').query("phase == 'B'")
+        assert len(b) == 1 and np.isfinite(b[['lpl', 'upl']].to_numpy()).all()
+        assert result.charts['Xbar']['metadata']['single_point_phases'] == 1
+        assert result.phase_statistics('Xbar').set_index('phase').loc['B', 'subgroups'] == 1
+
+    def test_insufficient_phase_keeps_its_points_without_limits(self):
+        stages = {'A': [[10, 12, 14], [11, 13, 12], [9, 10, 14]], 'B': [[20], [19], [21]]}
+        study = pb.formulate(_staged(stages), response='y', time='time', factors=['stage'])
+        result = study.execute(chart='Xbar', phased=True, companion=True)
+        xbar, s = result.get_chart('Xbar'), result.get_chart('S')
+        b = xbar[xbar['phase'] == 'B']
+        assert len(b) == 3
+        assert b[['center', 'lpl', 'upl']].isna().all().all()
+        assert (b['beyond_limits'] == 0).all()
+        assert 'B' not in set(s['phase'])
+        for chart in ('Xbar', 'S'):
+            assert result.charts[chart]['metadata']['insufficient_phases'] == ['B']
+        row = result.phase_statistics('Xbar').set_index('phase').loc['B']
+        assert row['subgroups'] == 3 and np.isnan(row['center']) and np.isnan(row['upl'])
+        assert 'B' not in set(result.phase_statistics('S')['phase'])
+        flagged = result.detect_signals('Xbar').flagged_observations
+        assert flagged.isdisjoint(set(b.index))
+
+    def test_no_replicated_phase_raises(self):
+        stages = {'A': [[10], [11], [9]], 'B': [[20], [19], [21]]}
+        study = pb.formulate(_staged(stages), response='y', time='time', factors=['stage'])
+        with pytest.raises(ValidationError):
+            study.execute(chart='Xbar', phased=True)
+
+
+# =============================================================================
+# Contracts
+# =============================================================================
+
+
+class TestContracts:
+    def test_solo_charts_equal_the_companion_pair(self, reference_studies):
+        study = reference_studies['PM SDS 3']
+        pair = study.execute(chart='Xbar', phased=True, companion=True)
+        for chart in ('Xbar', 'S'):
+            solo = study.execute(chart=chart, phased=True)
+            pd.testing.assert_frame_equal(solo.get_chart(chart), pair.get_chart(chart))
+            assert solo.get_statistics(chart) == pair.get_statistics(chart)
+            assert solo.charts[chart]['metadata'] == pair.charts[chart]['metadata']
+            pd.testing.assert_frame_equal(solo.phase_statistics(chart), pair.phase_statistics(chart))
+
+    def test_statistics_and_metadata_shape(self, sds1_study):
+        result = sds1_study.execute(chart='Xbar', phased=True, companion=True, n_mode='average')
+        for chart, value_col in (('Xbar', 'xbar'), ('S', 's')):
+            assert result.get_statistics(chart) == {
+                'N': None,
+                'center': None,
+                'lpl': None,
+                'upl': None,
+                'limits_vary': True,
+            }
+            meta = result.charts[chart]['metadata']
+            assert set(meta) == {
+                'chart_type',
+                'value_col',
+                'center_col',
+                'n_sigma',
+                'n_mode',
+                'group_col',
+                'phased',
+                'run_rules_applicable',
+                'phase_by',
+                'lane_boundaries',
+                'single_point_phases',
+                'insufficient_phases',
+            }
+            assert meta['value_col'] == value_col and meta['group_col'] == 'group'
+            assert meta['phased'] is True and meta['run_rules_applicable'] is False
+            assert meta['phase_by'] == ['factor 1', 'factor 2']
+            data = result.get_chart(chart)
+            assert list(data.columns) == ['group', 'phase', 'n', value_col, 'center', 'lpl', 'upl', 'beyond_limits']
+            assert list(result.phase_statistics(chart).columns) == [
+                'phase',
+                'subgroups',
+                'N',
+                'center',
+                's_bar',
+                'sigma_hat',
+                'lpl',
+                'upl',
+            ]
+
+    def test_lane_positions_are_the_first_rows_of_phases(self, reference_studies):
+        result = reference_studies['PM SDS 3'].execute(chart='Xbar', phased=True, companion=True)
+        for chart in ('Xbar', 'S'):
+            data = result.get_chart(chart).reset_index(drop=True)
+            starts = [i for i in range(1, len(data)) if data.loc[i, 'phase'] != data.loc[i - 1, 'phase']]
+            boundaries = result.charts[chart]['metadata']['lane_boundaries']
+            assert [b['position'] for b in boundaries] == starts
+            assert [b['label'] for b in boundaries] == [data.loc[i, 'phase'] for i in starts]
+        xbar_positions = [b['position'] for b in result.charts['Xbar']['metadata']['lane_boundaries']]
+        s_positions = [b['position'] for b in result.charts['S']['metadata']['lane_boundaries']]
+        assert xbar_positions != s_positions  # one-reading subgroups shift the S chart's dividers
+
+    def test_chart_table_has_phase_after_subgroup(self, sds1_study):
+        result = sds1_study.execute(chart='Xbar', phased=True, companion=True)
+        for chart in ('Xbar', 'S'):
+            table = result.chart_table(chart)
+            assert list(table.columns) == ['subgroup', 'phase', 'n', 'value', 'center', 'lpl', 'upl', 'signal']
+
+    def test_phase_statistics_refuses_unphased_and_xmr(self, sds1_study):
+        with pytest.raises(ValidationError, match='not phased'):
+            sds1_study.execute(chart='Xbar').phase_statistics('Xbar')
+        with pytest.raises(ValidationError, match='phased Xbar and S'):
+            sds1_study.execute(chart='X', by=[], phased=True).phase_statistics('X')
+
+    def test_detect_signals_rule_1_against_each_phase(self, reference_studies):
+        result = reference_studies['PM SDS 1'].execute(chart='Xbar', phased=True, companion=True)
+        for chart in ('Xbar', 'S'):
+            signals = result.detect_signals(chart)
+            data = result.get_chart(chart)
+            assert signals.rules_evaluated == ['rule_1']
+            assert not signals.is_partial
+            assert signals.flagged_observations == set(data.index[data['beyond_limits'] != 0])
+
+    def test_stage_order_follows_the_labels(self):
+        readings = [[10, 12, 14], [11, 13, 12]]
+        by_name = _staged({'before': readings, 'after': [[20, 22, 21], [19, 23, 22]]})
+        result = pb.formulate(by_name, response='y', time='time', factors=['stage']).execute(chart='Xbar', phased=True)
+        assert list(result.phase_statistics('Xbar')['phase']) == ['after', 'before']  # strings sort
+
+        sortable = by_name.assign(stage=by_name['stage'].map({'before': '1 before', 'after': '2 after'}))
+        result = pb.formulate(sortable, response='y', time='time', factors=['stage']).execute(chart='Xbar', phased=True)
+        assert list(result.phase_statistics('Xbar')['phase']) == ['1 before', '2 after']
+
+        ordered = by_name.assign(stage=pd.Categorical(by_name['stage'], categories=['before', 'after'], ordered=True))
+        result = pb.formulate(ordered, response='y', time='time', factors=['stage']).execute(chart='Xbar', phased=True)
+        assert list(result.phase_statistics('Xbar')['phase']) == ['before', 'after']
+
+    def test_unphased_output_is_unchanged(self, reference_studies):
+        study = reference_studies['PM SDS 3']
+        plain = study.execute(chart='Xbar', companion=True)
+        phased = study.execute(chart='Xbar', phased=True, companion=True)
+        for chart, value_col in (('Xbar', 'xbar'), ('S', 's')):
+            a = plain.get_chart(chart).reset_index(drop=True)
+            b = phased.get_chart(chart).reset_index(drop=True)
+            assert list(a['group']) == list(b['group'])
+            pd.testing.assert_series_equal(a[value_col], b[value_col])
+            pd.testing.assert_series_equal(a['n'], b['n'])
+            assert 'phase' not in a.columns
+            assert 'phased' not in plain.charts[chart]['metadata']
+
+    @pytest.mark.slow
+    def test_a_thousand_phases_is_quick(self):
+        rng = np.random.default_rng(0)
+        rows = [
+            {'stage': f'{k:04d}', 'time': t, 'y': rng.normal(k % 7, 1)}
+            for k in range(1000)
+            for t in range(3)
+            for _ in range(2)
+        ]
+        study = pb.formulate(pd.DataFrame(rows), response='y', time='time', factors=['stage'])
+        started = _time.perf_counter()
+        result = study.execute(chart='Xbar', phased=True, companion=True)
+        assert _time.perf_counter() - started < 5
+        assert len(result.phase_statistics('Xbar')) == 1000
+
+
+# =============================================================================
+# Validation
+# =============================================================================
+
+
+class TestValidation:
+    @pytest.mark.parametrize(
+        ('kwargs', 'fragment'),
+        [
+            (dict(chart='Histogram'), 'only valid for'),
+            (dict(chart='Xbar', by=[TIME]), 'already gives each factor combination its own chart'),
+            (dict(chart='S', by=[TIME]), 'use by=None with phased=True'),
+            (dict(chart='Xbar', by=['FACTOR 1']), 'requires the time variable in by='),
+            (dict(chart='Xbar', by=FACTORS), 'requires the time variable in by='),
+            (dict(chart='Xbar', by=[TIME, 'FACTOR 1']), 'time variable last in by='),
+            (dict(chart='Xbar', value='R5'), 'charts the response'),
+            (dict(chart='S', value='R1', recentered=True), 'charts the response'),
+        ],
+    )
+    def test_refused_requests(self, reference_studies, kwargs, fragment):
+        with pytest.raises(ValidationError, match=fragment):
+            reference_studies['PM SDS 1'].execute(phased=True, **kwargs)
+
+    def test_refused_without_factors(self, reference_df):
+        study = pb.formulate(reference_df, response='PM SDS 1', time=TIME)
+        with pytest.raises(ValidationError, match='requires factors to define phases'):
+            study.execute(chart='Xbar', phased=True)
+
+    @pytest.mark.parametrize('by', [None, [], FACTORS + [TIME], ['FACTOR 2', 'FACTOR 1', TIME], ['FACTOR 1', TIME]])
+    def test_accepted_shapes(self, reference_studies, by):
+        result = reference_studies['PM SDS 1'].execute(chart='Xbar', by=by, phased=True, companion=True)
+        assert result.charts['Xbar']['metadata']['phased'] is True
+        expected_phase_by = ['FACTOR 1'] if by == ['FACTOR 1', TIME] else FACTORS
+        assert result.charts['Xbar']['metadata']['phase_by'] == expected_phase_by
+
+    def test_calibration_is_refused_with_the_chart_in_the_context(self, reference_studies):
+        study = reference_studies['PM SDS 1']
+        cal = Calibration(label='ref', mean=10.0, sigma=1.0)
+        for chart in ('Xbar', 'S'):
+            with pytest.raises(CalibrationNotSupportedError) as exc:
+                study.execute(chart=chart, phased=True, calibration=cal)
+            assert exc.value.context == f'phased {chart}'
+            assert study.supports_calibration(chart=chart, phased=True) is False
+
+
+# =============================================================================
+# Plotting
+# =============================================================================
+
+
+@pytest.mark.plotting
+class TestPlotting:
+    def test_stepped_lines_dividers_and_note(self, cavity_study):
+        result = cavity_study.execute(chart='Xbar', phased=True, companion=True)
+        for chart in ('Xbar', 'S'):
+            fig = result.plot(chart=chart)
+            stepped = [t.name for t in fig.data if t.name and 'varies' in t.name]
+            assert any('CL' in name for name in stepped)
+            assert any('UPL' in name for name in stepped)
+            assert any('LPL' in name for name in stepped)
+            texts = [a.text for a in fig.layout.annotations]
+            assert 'Process limits computed per phase' in texts
+            assert {'B', 'C', 'D'} <= set(texts)  # divider labels name the phases
+            assert len(fig.layout.shapes) >= 3
+            assert not any(t.name and 'Zone' in t.name for t in fig.data)
diff --git a/tests/test_result_tabular.py b/tests/test_result_tabular.py
index 6274451..cfb6211 100644
--- a/tests/test_result_tabular.py
+++ b/tests/test_result_tabular.py
@@ -273,7 +273,7 @@ def test_signal_symbol_map_is_complete(self):
         assert set(_SIGNAL_SYMBOLS.keys()) == {-1, 0, 1}
 
     def test_non_value_cols_includes_typical_metadata(self):
-        for col in ("rsg", "center", "lpl", "upl", "beyond_limits", "n", "obs_id"):
+        for col in ("rsg", "center", "lpl", "upl", "beyond_limits", "n", "phase", "obs_id"):
             assert col in _NON_VALUE_COLS
 
     def test_preferred_statistic_cols_order(self):
diff --git a/tests/test_supports_calibration.py b/tests/test_supports_calibration.py
index bb58606..3ab8ad7 100644
--- a/tests/test_supports_calibration.py
+++ b/tests/test_supports_calibration.py
@@ -47,7 +47,7 @@ def _requests():
         for by in (None, [], ['FACTOR 1'], FACTORS, [TIME], FACTORS + [TIME]):
             for value in (None, 'R5', 'R2'):
                 for recentered in ((False, True) if value else (False,)):
-                    for phased in ((False, True) if chart in ('X', 'mR') else (False,)):
+                    for phased in (False, True):
                         kwargs = dict(
                             chart=chart, by=by, recentered=recentered, phased=phased,
                         )
@@ -136,6 +136,19 @@ def test_phased_xmr_is_refused_only_when_collapsing(study):
     assert study.supports_calibration(chart='X', by=[], phased=True) is False
 
 
+def test_phased_xbar_s_is_refused(study):
+    """Phased Xbar/S gives each phase its own limits, so one calibration cannot apply (#144)."""
+    assert study.supports_calibration(chart='Xbar', phased=True) is False
+    assert study.supports_calibration(chart='S', phased=True) is False
+    assert study.supports_calibration(chart='Xbar', by=['FACTOR 1', TIME], phased=True) is False
+    with pytest.raises(CalibrationNotSupportedError) as exc:
+        study.execute(chart='Xbar', phased=True, calibration=CAL)
+    assert exc.value.context == 'phased Xbar'
+    with pytest.raises(CalibrationNotSupportedError) as exc:
+        study.execute(chart='S', phased=True, calibration=CAL)
+    assert exc.value.context == 'phased S'
+
+
 def test_companion_does_not_change_the_answer(study):
     """Xbar+S together must agree with Xbar alone — they share the grouping."""
     for by in ([], [TIME], FACTORS):

From 398aedd0a246618a1859a6c75c35b0a1cbbbe990 Mon Sep 17 00:00:00 2001
From: Chris Nicholas <4948774+cnicholas@users.noreply.github.com>
Date: Sun, 18 Oct 2026 13:50:34 -0400
Subject: [PATCH 3/4] fix(data-preparation): X/mR lane order honours an ordered
 Categorical factor (#144)

The canonical row order in build_keys sorted the raw factor labels, so a stage
column with categories before < after gave Xbar/S subgroups in that order (the
aggregation groups on the factor column) but X/mR lanes, phases and strata in
the order after, before. The canonical sort now uses an ordered Categorical's
category codes. Every other input sorts exactly as before: the parallel key is
built only when an ordered Categorical is present (xbar_s_sweep: 0 differing
lines; golden masters untouched; gate 1547 / 0 / 4; full suite 2638 passed).
---
 CHANGELOG.md                                 |  4 ++
 processbehavior/data_preparation.py          | 15 +++-
 tests/test_ordered_categorical_lane_order.py | 73 ++++++++++++++++++++
 3 files changed, 90 insertions(+), 2 deletions(-)
 create mode 100644 tests/test_ordered_categorical_lane_order.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4f64ca9..e8863f7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -58,6 +58,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - **`execute()` no longer documents `phased=True` as requiring `by=[]`.** Stratified phasing of X/mR
   (`by=['factor 1'], phased=True`) has been accepted since 75bcdd1; the docstring never caught up. (#144)
 - **The companion mR of a phased X chart reports `single_point_phases`.** Only the solo mR did. (#144)
+- **X/mR lane order honours an ordered pandas Categorical factor.** A stage column with categories
+  `before < after` already ordered the Xbar/S subgroups that way, but the canonical row order sorted
+  the raw labels, so X/mR lanes, phases and `strata` came out `after, before`. The canonical sort now
+  uses an ordered Categorical's category codes; plain string and numeric factors sort as before. (#144)
 
 ## [0.4.0] - 2026-10-14
 
diff --git a/processbehavior/data_preparation.py b/processbehavior/data_preparation.py
index 0e3db11..f9e4bfa 100644
--- a/processbehavior/data_preparation.py
+++ b/processbehavior/data_preparation.py
@@ -529,8 +529,19 @@ def build_keys(self, df: pd.DataFrame, spec: FormulationSpec) -> pd.DataFrame:
         else:
             out['cell_key'] = out['rsg_key']
 
-        # 3. Sort by canonical order: (cell_key, obs_id)
-        out = out.sort_values(['cell_key', 'obs_id'], kind='stable')
+        # 3. Sort by canonical order: (cell_key, obs_id). An ordered Categorical factor sorts by
+        #    its category order rather than its labels, so a stage labelled before < after keeps
+        #    that order on every chart and in every stratum list, not only in the Xbar/S
+        #    aggregation (issue #144). Every other column sorts as the tuple of its values,
+        #    exactly as before: the parallel key exists only when an ordered Categorical is present.
+        key_cols = (k_vars + [t]) if (k_vars and t) else ([t] if t else list(k_vars))
+        ordered_cols = {c for c in key_cols if isinstance(out[c].dtype, pd.CategoricalDtype) and out[c].cat.ordered}
+        if ordered_cols:
+            parts = [out[c].cat.codes.astype(object) if c in ordered_cols else out[c].astype(object) for c in key_cols]
+            out['_order_key'] = list(map(tuple, pd.concat(parts, axis=1).to_numpy()))
+            out = out.sort_values(['_order_key', 'obs_id'], kind='stable').drop(columns=['_order_key'])
+        else:
+            out = out.sort_values(['cell_key', 'obs_id'], kind='stable')
 
         # 4. Assign sort_key AFTER sorting (canonical ordering 0..N-1)
         out['sort_key'] = np.arange(len(out), dtype=np.int64)
diff --git a/tests/test_ordered_categorical_lane_order.py b/tests/test_ordered_categorical_lane_order.py
new file mode 100644
index 0000000..259cd1e
--- /dev/null
+++ b/tests/test_ordered_categorical_lane_order.py
@@ -0,0 +1,73 @@
+"""An ordered pandas Categorical factor orders the lanes on every chart (issue #144, PR 4).
+
+The Xbar/S aggregation groups on the factor column and so followed an ordered Categorical's
+category order; the canonical row order in ``DataPreparation.build_keys`` sorted the raw labels,
+so X/mR lanes, phases and ``strata`` disagreed (``after, before`` for a stage ``before < after``).
+Plain string and numeric factors keep their previous order, which the golden masters and the
+Bishop gate pin.
+"""
+
+from __future__ import annotations
+
+import pandas as pd
+
+import processbehavior as pb
+
+
+def _staged_frame(labels=('before', 'after')) -> pd.DataFrame:
+    rows = []
+    for base, stage in zip((10, 20), labels, strict=True):
+        for t in range(1, 4):
+            for v in (0, 1, 2):
+                rows.append({'stage': stage, 'time': t, 'y': float(base + t + v)})
+    return pd.DataFrame(rows)
+
+
+def _lane_order(study, **kwargs) -> list[str]:
+    data = study.execute(chart='X', by=[], phased=True, **kwargs).get_chart('X')
+    return list(dict.fromkeys(data['phase']))
+
+
+class TestOrderedCategorical:
+    def test_every_view_follows_the_category_order(self):
+        df = _staged_frame()
+        ordered = df.assign(stage=pd.Categorical(df['stage'], categories=['before', 'after'], ordered=True))
+        study = pb.formulate(ordered, response='y', time='time', factors=['stage'])
+
+        assert _lane_order(study) == ['before', 'after']
+        assert study.execute(chart='X', by=['stage']).strata == ['before', 'after']
+        xbar = study.execute(chart='Xbar', phased=True)
+        assert list(xbar.phase_statistics('Xbar')['phase']) == ['before', 'after']
+        boundaries = study.execute(chart='X', by=[]).charts['X']['metadata']['lane_boundaries']
+        assert [(b['position'], b['label']) for b in boundaries] == [(9, 'after')]
+
+        frame = study._ads.analysis_dataset.sort_values('sort_key')
+        assert list(frame['stage'].astype(str).unique()) == ['before', 'after']
+        assert list(frame['sort_key']) == list(range(len(frame)))
+
+    def test_plain_strings_still_sort_lexically(self):
+        study = pb.formulate(_staged_frame(), response='y', time='time', factors=['stage'])
+        assert _lane_order(study) == ['after', 'before']
+        assert study.execute(chart='X', by=['stage']).strata == ['after', 'before']
+        assert list(study.execute(chart='Xbar', phased=True).phase_statistics('Xbar')['phase']) == ['after', 'before']
+
+    def test_sortable_labels_need_no_categorical(self):
+        df = _staged_frame(labels=('1 before', '2 after'))
+        study = pb.formulate(df, response='y', time='time', factors=['stage'])
+        assert _lane_order(study) == ['1 before', '2 after']
+
+    def test_two_factors_with_one_ordered(self):
+        df = _staged_frame()
+        two = pd.concat([df.assign(machine=1), df.assign(machine=2, y=df['y'] + 0.5)], ignore_index=True)
+        two['stage'] = pd.Categorical(two['stage'], categories=['before', 'after'], ordered=True)
+        study = pb.formulate(two, response='y', time='time', factors=['stage', 'machine'])
+        assert _lane_order(study) == ['before_1', 'before_2', 'after_1', 'after_2']
+        assert study.execute(chart='X', by=['stage']).strata == ['before', 'after']
+
+    def test_rows_within_a_cell_keep_source_order(self):
+        df = _staged_frame()
+        ordered = df.assign(stage=pd.Categorical(df['stage'], categories=['before', 'after'], ordered=True))
+        study = pb.formulate(ordered, response='y', time='time', factors=['stage'])
+        frame = study._ads.analysis_dataset.sort_values('sort_key')
+        for _, cell in frame.groupby(['stage', 'time'], observed=True, sort=False):
+            assert list(cell['obs_id']) == sorted(cell['obs_id'])

From 3454269b606cab3417514621d972087d3e68fb92 Mon Sep 17 00:00:00 2001
From: Chris Nicholas <4948774+cnicholas@users.noreply.github.com>
Date: Sun, 18 Oct 2026 13:47:15 -0400
Subject: [PATCH 4/4] feat(charts): phase column and phase_statistics on phased
 X/mR (#144)

Phased X and mR charts gain what phased Xbar/S got: a phase column (the lane
label) on the chart data and in chart_table, and result.phase_statistics(chart)
with one row per phase (phase, points, center, mR_bar, lpl, upl). A stratified
phased chart keeps one table per stratum, read with stratum= or through
focus(). The companion mR of a phased X chart now reports single_point_phases,
which only the solo mR did.

Output additions only: the per-phase arithmetic, row order, limits and signals
are unchanged, and phased=False charts are untouched (xbar_s_sweep: 0 differing
lines; gate 1547 / 0 / 4; full suite 2633 passed).
---
 CHANGELOG.md                       |  16 +++--
 docs/reference/api.md              |  24 +++----
 processbehavior/analysis.py        |  69 ++++++++++++++++---
 processbehavior/analysis_result.py |  55 ++++++++++++----
 processbehavior/focused_result.py  |   6 ++
 tests/test_phased_limits.py        | 102 ++++++++++++++++++++++++++---
 tests/test_phased_xbar_s.py        |   6 +-
 7 files changed, 227 insertions(+), 51 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index edd254d..4f64ca9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -25,8 +25,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   Minitab's staged Xbar-S chart. Residual values, calibrations and `by=[time]` are refused with
   messages that say what to use instead; a phase with no replicated subgroup keeps its points without
   limits and is named in `metadata['insufficient_phases']`. (#144)
-- **`result.phase_statistics(chart)`.** One row per phase of a phased Xbar/S chart: `phase`,
-  `subgroups`, `N`, `center`, `s_bar`, `sigma_hat`, `lpl`, `upl`. (#144)
+- **`result.phase_statistics(chart, stratum=None)`.** One row per phase of a phased chart. Xbar/S:
+  `phase`, `subgroups`, `N`, `center`, `s_bar`, `sigma_hat`, `lpl`, `upl`. X/mR: `phase`, `points`,
+  `center`, `mR_bar`, `lpl`, `upl`, per stratum on a stratified chart (`stratum=`, or through
+  `focus()`). (#144)
 
 ### Changed
 - **`formulate(response=..., factors=[])` with no time is the same study.** It used to treat the
@@ -36,12 +38,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - **Xbar and S chart data carry `n`.** `get_chart('Xbar')` / `get_chart('S')`, and the Excel chart
   sheets, gain an `n` column after the subgroup label. Unstratified Xbar/S metadata names the label
   column as `group_col` (#143).
-- **Phased Xbar/S chart data and `chart_table` carry a `phase` column** after the subgroup label, and
-  the metadata carries `phase_by`, `lane_boundaries`, `single_point_phases` and `insufficient_phases`.
-  `n_avg` is not reported on a phased chart, whose `n_mode='average'` is applied per phase; the
-  per-phase mean size is `phase_statistics(chart)['N']`. (#144)
 - **Lane dividers on an Xbar or S chart are labelled with the phase.** These charts have no time
   column and so no cell band under the axis, which is where X charts name their lanes. (#144)
+- **Phased chart data and `chart_table` carry a `phase` column** (the lane label) after the subgroup
+  label, on X/mR as well as Xbar/S. Phased Xbar/S metadata carries `phase_by`, `lane_boundaries`,
+  `single_point_phases` and `insufficient_phases`; `n_avg` is not reported on a phased chart, whose
+  `n_mode='average'` is applied per phase, so the per-phase mean size is `phase_statistics(chart)['N']`.
+  (#144)
 
 ### Fixed
 - **Xbar/S `chart_table` keeps the subgroup label and `n` on every `by=`.** With a `by=` other than
@@ -54,6 +57,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   centres and statistics are unchanged (#143).
 - **`execute()` no longer documents `phased=True` as requiring `by=[]`.** Stratified phasing of X/mR
   (`by=['factor 1'], phased=True`) has been accepted since 75bcdd1; the docstring never caught up. (#144)
+- **The companion mR of a phased X chart reports `single_point_phases`.** Only the solo mR did. (#144)
 
 ## [0.4.0] - 2026-10-14
 
diff --git a/docs/reference/api.md b/docs/reference/api.md
index 06ababb..24b6b5b 100644
--- a/docs/reference/api.md
+++ b/docs/reference/api.md
@@ -330,17 +330,19 @@ directly, and is not a top-level export.)
 ### Phases
 
 ```python
-result.phase_statistics(chart)   # -> DataFrame, one row per phase (phased Xbar/S)
-```
-
-A chart executed with `phased=True` gives each phase — one combination of the grouping
-factors, a contiguous run of subgroups in time order — its own center line and limits.
-On Xbar/S the chart data gains a `phase` column (the lane label) next to the subgroup
-label, and `chart_table` shows it. `phase_statistics(chart)` returns the numbers behind
-the stepped lines: `phase`, `subgroups`, `N` (the subgroup size the phase's limits assume;
-`NaN` when sizes vary within it), `center`, `s_bar`, `sigma_hat`, `lpl`, `upl` (`NaN` when
-the limits step within the phase). It raises `ValidationError` on a chart that is not
-phased, and on phased X/mR, whose per-phase values are on each row of `get_chart()`.
+result.phase_statistics(chart, stratum=None)   # -> DataFrame, one row per phase
+```
+
+A chart executed with `phased=True` gives each phase — one factor combination, a
+contiguous run of points in time order — its own center line and limits. The chart data
+gains a `phase` column (the lane label) next to the subgroup label, and `chart_table`
+shows it. `phase_statistics(chart)` returns the numbers behind the stepped lines. On
+Xbar/S: `phase`, `subgroups`, `N` (the subgroup size the phase's limits assume; `NaN`
+when sizes vary within it), `center`, `s_bar`, `sigma_hat`, `lpl`, `upl` (`NaN` when the
+limits step within the phase). On X/mR: `phase`, `points`, `center`, `mR_bar`, `lpl`,
+`upl`. A stratified phased X/mR chart (`by=['factor 1'], phased=True`) has phases per
+stratum: pass `stratum=`, or call it on `focus(stratum)`. It raises `ValidationError` on a
+chart that is not phased.
 
 Phased metadata carries `phased=True`, `run_rules_applicable=False`, `phase_by` (the
 factors a phase is one combination of), `lane_boundaries` (the divider positions),
diff --git a/processbehavior/analysis.py b/processbehavior/analysis.py
index c8a22b8..cb68c48 100644
--- a/processbehavior/analysis.py
+++ b/processbehavior/analysis.py
@@ -332,6 +332,25 @@ def _raise_no_replication(chart: str, sds, insufficient: list, total: int) -> No
     )
 
 
+def _mr_phase_table(phase_stats: pd.DataFrame, out: pd.DataFrame, keys: list[str], round_to: int) -> pd.DataFrame:
+    """One row per phase of a phased X or mR chart: label, points, centre, mR-bar and limits.
+
+    ``phase_stats`` is the per-phase aggregation keyed by ``keys`` (``['_phase_id']``, or
+    ``[stratify_col, '_phase_id']``); ``out`` carries the ``phase`` label on every row.
+    """
+    labels = out.groupby(keys, sort=False)['phase'].first().reset_index()
+    table = phase_stats.merge(labels, on=keys, how='left').rename(columns={'_n': 'points', 'mR': 'mR_bar'})
+    return table[keys + ['phase', 'points', 'center', 'mR_bar', 'lpl', 'upl']].round(round_to)
+
+
+def _split_phase_table(table: pd.DataFrame, stratify_col: str) -> dict:
+    """The per-stratum phase tables of a stratified phased chart, keyed by stratum."""
+    return {
+        stratum: group.drop(columns=[stratify_col, '_phase_id']).reset_index(drop=True)
+        for stratum, group in table.groupby(stratify_col, sort=False)
+    }
+
+
 def _build_stratified_mr_payload(
     *,
     mr_spec: _MRChartSpec,
@@ -345,6 +364,7 @@ def _build_stratified_mr_payload(
     phased: bool = False,
     insufficient_strata: list | None = None,
     single_point_phases: int | None = None,
+    phases: dict | None = None,
 ) -> dict:
     """Build the {chart_type: {...}} payload for a stratified MR-family chart.
 
@@ -379,6 +399,7 @@ def _build_stratified_mr_payload(
             'statistics': chart_statistics,
             'metadata': metadata,
             'strata': published_strata,
+            **({'phases': phases} if phased and phases is not None else {}),
         },
     }
 
@@ -392,6 +413,7 @@ def _build_ungrouped_mr_payload(
     lane_boundaries,
     phased: bool = False,
     single_point_phases: int | None = None,
+    phases: pd.DataFrame | None = None,
 ) -> dict:
     """Build the {chart_type: {...}} payload for an ungrouped MR-family chart."""
     metadata: dict = {
@@ -410,6 +432,7 @@ def _build_ungrouped_mr_payload(
             'data': chart_out,
             'statistics': chart_statistics,
             'metadata': metadata,
+            **({'phases': phases} if phased and phases is not None else {}),
         },
     }
 
@@ -2239,10 +2262,13 @@ def _calculate_mr_chart_from_precomputed(
             phase_stats = precomputed['phase_stats']
 
             # Recompute R-specific limits per phase within each stratum
-            r_phase_stats = phase_stats[[stratify_col, '_phase_id', 'mR']].copy()
+            r_phase_stats = phase_stats[[stratify_col, '_phase_id', '_n', 'mR']].copy()
             r_lims = _mr_family_limits(mr_spec, mean=0, mR=r_phase_stats['mR'])
             r_phase_stats = _join_limits(r_phase_stats, r_lims)
             r_phase_stats['center'] = r_phase_stats['mR']
+            r_phases = _split_phase_table(
+                _mr_phase_table(r_phase_stats, out, [stratify_col, '_phase_id'], spec.round_to), stratify_col
+            )
 
             # Replace X limits with mR limits
             out = out.drop(columns=['center', 'lpl', 'upl', 'beyond_limits'], errors='ignore')
@@ -2261,7 +2287,7 @@ def _calculate_mr_chart_from_precomputed(
             # Format output
             chart_out = self._build_output_columns(
                 df=out,
-                value_cols=extra_cols + [plot_col, 'center', 'lpl', 'upl', 'beyond_limits'],
+                value_cols=extra_cols + ['phase', plot_col, 'center', 'lpl', 'upl', 'beyond_limits'],
             )
 
             chart_statistics = {
@@ -2286,6 +2312,8 @@ def _calculate_mr_chart_from_precomputed(
                 strata=strata,
                 lane_boundaries=xmr_lane_boundaries,
                 phased=True,
+                single_point_phases=precomputed.get('single_point_phases'),
+                phases=r_phases,
             )
 
         # --- Non-phased precomputed path ---
@@ -2396,12 +2424,14 @@ def _calculate_mr_chart_stratified(  # noqa: C901
             else:
                 out['_phase_key'] = list(map(tuple, out[collapsed_factors].to_numpy()))
 
-            # 2. Assign within-stratum phase IDs via transitions of _phase_key
+            # 2. Assign within-stratum phase IDs via transitions of _phase_key, and the
+            #    phase's display label (the lane label) on every row
             out['_phase_id'] = (
                 out.groupby(stratify_col, sort=False)['_phase_key']
                 .transform(lambda s: (s != s.shift()).cumsum())
                 .astype('int64')
             )
+            out['phase'] = _lane_labels(out, collapsed_factors)
 
             # 3. Per-phase mR (reset at phase boundaries within each stratum)
             out['mr'] = out.groupby([stratify_col, '_phase_id'], sort=False)[mr_source_col].diff().abs()
@@ -2455,12 +2485,15 @@ def _calculate_mr_chart_stratified(  # noqa: C901
 
             # 11. Track single-point phases
             n_single_point = int((phase_stats['_n'] < 2).sum())
+            phases = _split_phase_table(
+                _mr_phase_table(phase_stats, out, [stratify_col, '_phase_id'], spec.round_to), stratify_col
+            )
 
             # Format output
             extra_cols = [c for c in stratify_by if c != spec.rsg_var_name and c != spec.time_var]
             chart_out = self._build_output_columns(
                 df=out,
-                value_cols=extra_cols + [plot_col, 'center', 'lpl', 'upl', 'beyond_limits'],
+                value_cols=extra_cols + ['phase', plot_col, 'center', 'lpl', 'upl', 'beyond_limits'],
             )
 
             result.update(_build_stratified_mr_payload(
@@ -2474,12 +2507,14 @@ def _calculate_mr_chart_stratified(  # noqa: C901
                 lane_boundaries=lane_boundaries,
                 phased=True,
                 single_point_phases=n_single_point,
+                phases=phases,
             ))
 
             if _return_intermediates:
                 result['_intermediates'] = {
                     'out': out,
                     'phase_stats': phase_stats,
+                    'single_point_phases': n_single_point,
                     'stratify_col': stratify_col,
                     'stratify_by': stratify_by,
                     'strata': strata,
@@ -2614,7 +2649,7 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
                 xmr_lane_boundaries = _precomputed['lane_boundaries']
 
                 # Recompute R-specific limits per phase
-                r_phase_stats = phase_stats[['_phase_id', 'mR']].copy()
+                r_phase_stats = phase_stats[['_phase_id', '_n', 'mR']].copy()
                 r_lims = r_phase_stats.apply(
                     lambda row: calculate_limits(
                         mean=0,
@@ -2628,6 +2663,11 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
                 )
                 r_phase_stats = _join_limits(r_phase_stats, r_lims)
                 r_phase_stats['center'] = r_phase_stats['mR']
+                r_phases = (
+                    _mr_phase_table(r_phase_stats, out, ['_phase_id'], spec.round_to)
+                    .drop(columns=['_phase_id'])
+                    .reset_index(drop=True)
+                )
 
                 # Replace X limits with mR limits
                 out = out.drop(columns=['center', 'lpl', 'upl'], errors='ignore')
@@ -2646,7 +2686,7 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
 
                 chart_out = self._build_output_columns(
                     df=out,
-                    value_cols=[plot_col, 'center', 'lpl', 'upl', 'beyond_limits'],
+                    value_cols=['phase', plot_col, 'center', 'lpl', 'upl', 'beyond_limits'],
                 )
 
                 chart_statistics = {
@@ -2665,6 +2705,8 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
                     plot_col=plot_col,
                     lane_boundaries=xmr_lane_boundaries,
                     phased=True,
+                    single_point_phases=_precomputed.get('single_point_phases'),
+                    phases=r_phases,
                 ))
                 return result
 
@@ -2738,8 +2780,10 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
         # === Phased limits path ===
         if phased and collapsed_factors:
             self._reject_calibration('phased X/mR')
-            # 1. Assign phase_id: contiguous runs of the same rsg_key
+            # 1. Assign phase_id: contiguous runs of the same rsg_key, and the phase's display
+            #    label (the lane label) on every row
             out['_phase_id'] = (out['rsg_key'] != out['rsg_key'].shift()).cumsum().astype('int64')
+            out['phase'] = _lane_labels(out, collapsed_factors)
 
             # 2. Per-phase moving range (reset at boundaries)
             out['mr'] = out.groupby('_phase_id', sort=False)[mr_source_col].diff().abs()
@@ -2806,7 +2850,7 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
             # 9. Format output
             chart_out = self._build_output_columns(
                 df=out,
-                value_cols=[plot_col, 'center', 'lpl', 'upl', 'beyond_limits'],
+                value_cols=['phase', plot_col, 'center', 'lpl', 'upl', 'beyond_limits'],
             )
 
             # 10. Statistics: limits_vary=True signals stepped rendering
@@ -2818,8 +2862,13 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
                 'limits_vary': True,
             }
 
-            # 11. Track single-point phases for visibility
+            # 11. Track single-point phases for visibility, and the per-phase table
             n_single_point = int((phase_stats['_n'] < 2).sum())
+            phases = (
+                _mr_phase_table(phase_stats, out, ['_phase_id'], spec.round_to)
+                .drop(columns=['_phase_id'])
+                .reset_index(drop=True)
+            )
 
             result.update(_build_ungrouped_mr_payload(
                 mr_spec=mr_spec,
@@ -2829,6 +2878,7 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
                 lane_boundaries=lane_boundaries,
                 phased=True,
                 single_point_phases=n_single_point,
+                phases=phases,
             ))
 
             if _return_intermediates:
@@ -2839,6 +2889,7 @@ def _calculate_mr_chart_ungrouped(  # noqa: C901
                     'is_stratified': False,
                     'phased': True,
                     'phase_stats': phase_stats,
+                    'single_point_phases': n_single_point,
                 }
 
             return result
diff --git a/processbehavior/analysis_result.py b/processbehavior/analysis_result.py
index 3a3f2b8..13cd389 100644
--- a/processbehavior/analysis_result.py
+++ b/processbehavior/analysis_result.py
@@ -656,21 +656,34 @@ def statistics(self, chart: str, stratum: str | None = None) -> dict:
             )
         return dict(stats[stratum])
 
-    def phase_statistics(self, chart: str) -> pd.DataFrame:
+    def phase_statistics(self, chart: str, stratum: str | None = None) -> pd.DataFrame:
         """One row per phase of a phased chart: the numbers behind its stepped lines.
 
-        Columns: ``phase`` (the lane label, as in ``chart_table``), ``subgroups`` (points
-        charted in it), ``N`` (the subgroup size its limits assume: the common size, the mean
-        size under ``n_mode='average'``, or NaN when sizes vary within the phase), ``center``,
-        ``s_bar``, ``sigma_hat``, ``lpl`` and ``upl`` (NaN when the limits step within the
-        phase; ``get_chart`` carries them per row). A phase with no replicated subgroup keeps
-        its ``subgroups`` and is NaN elsewhere; ``metadata['insufficient_phases']`` names it.
+        Xbar/S columns: ``phase`` (the lane label, as in ``chart_table``), ``subgroups``
+        (points charted in it), ``N`` (the subgroup size its limits assume: the common size,
+        the mean size under ``n_mode='average'``, or NaN when sizes vary within the phase),
+        ``center``, ``s_bar``, ``sigma_hat``, ``lpl`` and ``upl`` (NaN when the limits step
+        within the phase; ``get_chart`` carries them per row). A phase with no replicated
+        subgroup keeps its ``subgroups`` and is NaN elsewhere; ``metadata['insufficient_phases']``
+        names it.
+
+        X/mR columns: ``phase``, ``points``, ``center``, ``mR_bar``, ``lpl``, ``upl``. On the
+        mR chart ``center`` is ``mR_bar``.
+
+        Parameters
+        ----------
+        chart : str
+            Chart name, e.g. ``'Xbar'`` or ``'X'``.
+        stratum : str, optional
+            Required on a stratified phased chart (X/mR with ``by=[factor]``), where each
+            stratum has its own phases; rejected otherwise. ``focus(stratum)`` gives the same
+            table without the argument.
 
         Raises
         ------
         ValidationError
-            On a chart that is not phased, and on phased X/mR, whose per-phase values are on
-            each row of ``get_chart`` (a later release adds this table there; issue #144).
+            On a chart that is not phased, or when a stratum is needed and not given, given
+            and not needed, or unknown.
         ChartNotAvailableError
             If the chart is not in this result.
         """
@@ -680,16 +693,30 @@ def phase_statistics(self, chart: str) -> pd.DataFrame:
                 f"Chart '{name}' not found. Available charts: {self.all_charts}", chart=name, available=self.all_charts
             )
         payload = self.charts[name]
-        if not payload.get('metadata', {}).get('phased'):
+        if not payload.get('metadata', {}).get('phased') or payload.get('phases') is None:
             raise ValidationError(
                 f'Chart {chart!r} is not phased, so it has no per-phase statistics. '
                 'execute(..., phased=True) gives each phase its own centre line and limits.'
             )
-        phases = payload.get('phases')
-        if phases is None:
+        phases = payload['phases']
+        if isinstance(phases, dict):
+            available = list(phases)
+            if stratum is None:
+                raise ValidationError(
+                    f'Chart {chart!r} is stratified, so a stratum is required to return one set of '
+                    f'phases.\nAvailable strata: {available}\n'
+                    f'Use result.phase_statistics({chart!r}, stratum={available[0]!r}) '
+                    f'or result.focus({available[0]!r}).phase_statistics({chart!r}).'
+                )
+            if stratum not in phases:
+                raise ValidationError(
+                    f'Stratum {stratum!r} not found in chart {chart!r}.\nAvailable strata: {available}'
+                )
+            return phases[stratum].copy()
+        if stratum is not None:
             raise ValidationError(
-                'phase_statistics is available on phased Xbar and S charts; a phased X or mR chart '
-                f"carries each phase's centre and limits on its rows (get_chart({chart!r}))."
+                f'This result is not stratified, so stratum={stratum!r} does not apply.\n'
+                f'Use result.phase_statistics({chart!r}) with no stratum.'
             )
         return phases.copy()
 
diff --git a/processbehavior/focused_result.py b/processbehavior/focused_result.py
index 1bc4982..a6823e4 100644
--- a/processbehavior/focused_result.py
+++ b/processbehavior/focused_result.py
@@ -99,6 +99,11 @@ def _build_focused_chart_info(
     focused_stats = _extract_stratum_statistics(chart_name, chart_info, stratum)
     focused_lb = _extract_stratum_lane_boundaries(parent_metadata, stratum)
 
+    # A phased stratified chart keeps one per-phase table per stratum; the focused chart
+    # carries its own, in the unstratified shape phase_statistics() reads.
+    phases = chart_info.get("phases")
+    focused_phases = {"phases": phases[stratum]} if isinstance(phases, dict) and stratum in phases else {}
+
     return {
         "data": focused_data,
         "statistics": focused_stats,
@@ -108,6 +113,7 @@ def _build_focused_chart_info(
             "focused_stratum": stratum,
             "lane_boundaries": focused_lb,
         },
+        **focused_phases,
     }
 
 
diff --git a/tests/test_phased_limits.py b/tests/test_phased_limits.py
index 82b435d..685ee6d 100644
--- a/tests/test_phased_limits.py
+++ b/tests/test_phased_limits.py
@@ -470,16 +470,22 @@ def test_explicit_rules_run_against_each_phase_centre(self):
 
         # Phase 1 sits above its centre 0, phase 2 below its centre 10: one run each.
         # Phase 2's last two points fall in its lower zone A (7 to 8, sigma 1).
-        data = pd.DataFrame({
-            'x': [1.0] * 8 + [9.0] * 6 + [7.5, 7.5],
-            'center': [0.0] * 8 + [10.0] * 8,
-            'upl': [3.0] * 8 + [13.0] * 8,
-            'lpl': [-3.0] * 8 + [7.0] * 8,
-        })
+        data = pd.DataFrame(
+            {
+                'x': [1.0] * 8 + [9.0] * 6 + [7.5, 7.5],
+                'center': [0.0] * 8 + [10.0] * 8,
+                'upl': [3.0] * 8 + [13.0] * 8,
+                'lpl': [-3.0] * 8 + [7.0] * 8,
+            }
+        )
         stats = {'N': None, 'center': None, 'lpl': None, 'upl': None, 'limits_vary': True}
         signals = SignalDetector().detect(
-            data, stats, SignalConfig(enabled_rules=['rule_2', 'rule_4'], min_observations=1),
-            value_col='x', chart_type='X', run_rules_applicable=False,
+            data,
+            stats,
+            SignalConfig(enabled_rules=['rule_2', 'rule_4'], min_observations=1),
+            value_col='x',
+            chart_type='X',
+            run_rules_applicable=False,
         )
         assert signals.rules_evaluated == ['rule_2', 'rule_4']
         assert set(signals.by_rule['rule_4']['obs_id']) == {7, 15}
@@ -488,3 +494,83 @@ def test_explicit_rules_run_against_each_phase_centre(self):
     def test_unphased_x_keeps_every_rule(self, sds1_study):
         signals = sds1_study.execute(chart='X', by=[]).detect_signals(chart='X')
         assert signals.rules_evaluated == [f'rule_{i}' for i in range(1, 9)]
+
+
+# =============================================================================
+# Phase column and per-phase statistics (#144, parity with phased Xbar/S)
+# =============================================================================
+
+
+class TestPhaseColumnAndStatistics:
+    """Phased X/mR name each row's phase and expose one row per phase, as phased Xbar/S do."""
+
+    def test_phase_column_is_the_lane_label(self, sds1_study):
+        result = sds1_study.execute(chart='X', by=[], phased=True, companion=True)
+        for chart in ('X', 'mR'):
+            data = result.get_chart(chart)
+            assert list(data.columns) == [
+                'time',
+                'rsg',
+                'phase',
+                data.columns[3],
+                'center',
+                'lpl',
+                'upl',
+                'beyond_limits',
+                'obs_id',
+            ]
+            assert list(data['phase']) == list(data['rsg'].astype(str))
+            # one centre per phase, and every row of a phase shares it
+            assert data.groupby('phase', sort=False)['center'].nunique().eq(1).all()
+            assert list(result.chart_table(chart).columns[:3]) == ['Obs', 'subgroup', 'phase']
+
+    def test_phase_statistics_x_and_mr(self, sds1_study):
+        result = sds1_study.execute(chart='X', by=[], phased=True, companion=True)
+        x, mr = result.phase_statistics('X'), result.phase_statistics('mR')
+        assert list(x.columns) == ['phase', 'points', 'center', 'mR_bar', 'lpl', 'upl']
+        assert list(x.columns) == list(mr.columns)
+        data = result.get_chart('X')
+        assert list(x['phase']) == list(dict.fromkeys(data['phase']))  # display order, one row each
+        assert x['points'].sum() == len(data)
+        for _, row in x.iterrows():
+            rows = data[data['phase'] == row['phase']]
+            assert rows['center'].iloc[0] == row['center']
+            assert (rows['lpl'].iloc[0], rows['upl'].iloc[0]) == (row['lpl'], row['upl'])
+        assert (mr['center'] == mr['mR_bar']).all()
+        assert (mr['lpl'] == 0).all()
+        pd.testing.assert_series_equal(mr['mR_bar'], x['mR_bar'])
+
+    def test_companion_mr_carries_single_point_phases(self, sds1_study):
+        result = sds1_study.execute(chart='X', by=[], phased=True, companion=True)
+        x_meta, mr_meta = result.charts['X']['metadata'], result.charts['mR']['metadata']
+        assert 'single_point_phases' in mr_meta
+        assert mr_meta['single_point_phases'] == x_meta['single_point_phases']
+        assert x_meta['single_point_phases'] == int((result.phase_statistics('X')['points'] < 2).sum())
+
+    def test_stratified_phases_by_stratum_and_focus(self, sds1_study):
+        result = sds1_study.execute(chart='X', by=['factor 1'], phased=True, companion=True)
+        for chart in ('X', 'mR'):
+            assert 'phase' in result.get_chart(chart).columns
+            with pytest.raises(ValidationError, match='stratum is required'):
+                result.phase_statistics(chart)
+            with pytest.raises(ValidationError, match='not found'):
+                result.phase_statistics(chart, stratum='no such stratum')
+            for stratum in result.strata:
+                table = result.phase_statistics(chart, stratum=stratum)
+                assert list(table.columns) == ['phase', 'points', 'center', 'mR_bar', 'lpl', 'upl']
+                focused = result.focus(stratum)
+                pd.testing.assert_frame_equal(focused.phase_statistics(chart), table)
+                data = focused.get_chart(chart)
+                assert list(table['phase']) == list(dict.fromkeys(data['phase']))
+                assert table['points'].sum() == len(data)
+            assert 'single_point_phases' in result.charts[chart]['metadata']
+
+    def test_unphased_charts_have_no_phase(self, sds1_study):
+        result = sds1_study.execute(chart='X', by=[], companion=True)
+        for chart in ('X', 'mR'):
+            assert 'phase' not in result.get_chart(chart).columns
+            with pytest.raises(ValidationError, match='not phased'):
+                result.phase_statistics(chart)
+        phased = sds1_study.execute(chart='X', by=[], phased=True)
+        with pytest.raises(ValidationError, match='does not apply'):
+            phased.phase_statistics('X', stratum='F1_1_F2_1')
diff --git a/tests/test_phased_xbar_s.py b/tests/test_phased_xbar_s.py
index c4b3ebc..4eb51c7 100644
--- a/tests/test_phased_xbar_s.py
+++ b/tests/test_phased_xbar_s.py
@@ -411,11 +411,11 @@ def test_chart_table_has_phase_after_subgroup(self, sds1_study):
             table = result.chart_table(chart)
             assert list(table.columns) == ['subgroup', 'phase', 'n', 'value', 'center', 'lpl', 'upl', 'signal']
 
-    def test_phase_statistics_refuses_unphased_and_xmr(self, sds1_study):
+    def test_phase_statistics_refuses_unphased(self, sds1_study):
         with pytest.raises(ValidationError, match='not phased'):
             sds1_study.execute(chart='Xbar').phase_statistics('Xbar')
-        with pytest.raises(ValidationError, match='phased Xbar and S'):
-            sds1_study.execute(chart='X', by=[], phased=True).phase_statistics('X')
+        with pytest.raises(ValidationError, match='does not apply'):
+            sds1_study.execute(chart='Xbar', phased=True).phase_statistics('Xbar', stratum='1_1')
 
     def test_detect_signals_rule_1_against_each_phase(self, reference_studies):
         result = reference_studies['PM SDS 1'].execute(chart='Xbar', phased=True, companion=True)