From fa6d0de19f2558e668bd87f16cbb11d4a96082d5 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:53:58 +0000
Subject: [PATCH 01/15] docs: fix docstring in designs_response_surface.py

Add caveat noting that a numeric alpha is only honored for the
fractional-cube CCD path; cube='full' silently coerces it to 'orthogonal'.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/experiments/designs_response_surface.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/process_improve/experiments/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py
index bc7406c..8e5f364 100644
--- a/src/process_improve/experiments/designs_response_surface.py
+++ b/src/process_improve/experiments/designs_response_surface.py
@@ -47,6 +47,12 @@ def dispatch_ccd(  # noqa: PLR0913
         Axial distance.  Accepted string values: ``"rotatable"``,
         ``"face_centered"``, ``"orthogonal"``.  A numeric value sets
         alpha directly.  Defaults to ``"orthogonal"``.
+
+        .. note::
+           A numeric ``alpha`` is only honored when ``cube="fractional"``.
+           For ``cube="full"`` (the default) the underlying pyDOE3
+           ``ccdesign`` call does not accept an arbitrary axial distance,
+           so a numeric value is silently treated as ``"orthogonal"``.
     cube : str
         How to build the cube (factorial) portion: ``"full"`` (default) uses
         the complete 2^k factorial; ``"fractional"`` uses a resolution-V (or

From 89f6edf6f092e538eb216fd554fc65bee0a54d04 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:54:48 +0000
Subject: [PATCH 02/15] docs: fix docstring in designs.py

Propagate the same numeric-alpha caveat: a numeric alpha is only honored
for cube='fractional'; cube='full' silently coerces it to 'orthogonal'.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/experiments/designs.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/process_improve/experiments/designs.py b/src/process_improve/experiments/designs.py
index 1c97837..e16d0c5 100644
--- a/src/process_improve/experiments/designs.py
+++ b/src/process_improve/experiments/designs.py
@@ -332,6 +332,12 @@ def generate_design(  # noqa: PLR0913
     alpha : str, float, or None
         Axial distance for CCD designs: ``"rotatable"``,
         ``"face_centered"``, ``"orthogonal"``, or a numeric value.
+
+        .. note::
+           A numeric ``alpha`` is only honored when ``cube="fractional"``.
+           For ``cube="full"`` (the default) the underlying pyDOE3
+           ``ccdesign`` call does not accept an arbitrary axial distance,
+           so a numeric value is silently treated as ``"orthogonal"``.
     cube : str
         For CCD designs, how to build the cube (factorial) portion:
         ``"full"`` (default) uses the complete 2^k factorial; ``"fractional"``

From 6fe63ece7cc5fc8dcea5fc1743f7b01f91b50d29 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:55:20 +0000
Subject: [PATCH 03/15] docs: fix docstring in _pca.py

PCA.select_n_components: document return_consensus in Parameters and
add the four extra Returns keys it exposes.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/multivariate/_pca.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/src/process_improve/multivariate/_pca.py b/src/process_improve/multivariate/_pca.py
index ee9d632..a676f3a 100644
--- a/src/process_improve/multivariate/_pca.py
+++ b/src/process_improve/multivariate/_pca.py
@@ -1067,6 +1067,13 @@ def select_n_components(  # noqa: PLR0913, PLR0915, C901
             step. Ignored under ``cv_scheme="row_wise"``.
         random_state : int, optional
             Seed for the ekf element-fold permutation.
+        return_consensus : bool, default False
+            When ``True``, also cross-check the CV recommendation against
+            two cheap alternative selectors: Minka's PPCA MLE
+            (:meth:`minka_mle`) and Horn's parallel analysis
+            (:meth:`parallel_analysis`). The result Bunch then gains the
+            ``minka_n_components``, ``parallel_analysis_n_components``,
+            ``consensus``, and ``consensus_counts`` keys (see Returns).
         threshold : float, optional
             Deprecated. The original Wold PRESS-ratio cutoff. Passing it
             emits a :class:`DeprecationWarning`; the value is ignored. Use
@@ -1106,6 +1113,17 @@ def select_n_components(  # noqa: PLR0913, PLR0915, C901
             - ``cv_scheme`` - the scheme used (``"ekf"`` or ``"row_wise"``).
             - ``selection_rule`` - the rule used to pick ``n_components``.
 
+            When ``return_consensus=True``, the Bunch additionally carries:
+
+            - ``minka_n_components`` - the Minka PPCA MLE estimate (int).
+            - ``parallel_analysis_n_components`` - Horn's parallel-analysis
+              estimate (int).
+            - ``consensus`` - ``"agree"`` if the three integer estimates
+              (CV recommendation, Minka, parallel analysis) span at most
+              1, otherwise ``"disagree"``.
+            - ``consensus_counts`` - the tuple
+              ``(recommended, minka_n, parallel_analysis_n)``.
+
         References
         ----------
         Bro, R., Kjeldahl, K., Smilde, A. K., & Kiers, H. A. L. (2008).

From 8b317db60396512a31771d0df86867d97874c996 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:55:45 +0000
Subject: [PATCH 04/15] docs: fix docstring in _pls.py

PLS.select_n_components: complete the truncated selection_mode Returns
sentence (mode is None when selection_distribution is None).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/multivariate/_pls.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/process_improve/multivariate/_pls.py b/src/process_improve/multivariate/_pls.py
index 353ec74..cd4e25a 100644
--- a/src/process_improve/multivariate/_pls.py
+++ b/src/process_improve/multivariate/_pls.py
@@ -1311,7 +1311,7 @@ def select_n_components(  # noqa: C901, PLR0912, PLR0913, PLR0915
               distribution signals a confident recommendation; a flat or
               multi-modal one flags it for review.
             - ``selection_mode`` - the most-voted component count, or
-              ``None`` when ``selection_distribution`` is.
+              ``None`` when ``selection_distribution`` is ``None``.
             - ``selection_is_stable`` - ``True`` iff the modal vote share
               meets ``stability_threshold``; ``None`` when no distribution
               was computed.

From 2e32e2c2f94a94d53e70a4c1403b2e29ba21563b Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:56:45 +0000
Subject: [PATCH 05/15] docs: fix docstring in _robust_regression.py

robust_regression: add t_value to always-present keys and describe the
degenerate early-return dict contents.

multiple_linear_regression: align Returns with OLS.to_dict() output.
Add R2_regression_based, R2_residual_based, k, conf_interval_intercept;
document the degenerate-path shapes for coefficients and conf_intervals.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 .../regression/_robust_regression.py          | 19 +++++++++++++++++--
 1 file changed, 17 insertions(+), 2 deletions(-)

diff --git a/src/process_improve/regression/_robust_regression.py b/src/process_improve/regression/_robust_regression.py
index 4333d81..4b19eb5 100644
--- a/src/process_improve/regression/_robust_regression.py
+++ b/src/process_improve/regression/_robust_regression.py
@@ -130,11 +130,19 @@ def robust_regression(  # noqa: PLR0913, PLR0915
         k:                        the number of model parameters (2 if fit_intercept else 1)
         fitted_values:            the N predicted values, one per row in y
         residuals:                the N residuals
+        t_value:                  the two-sided critical t-value at ``conflevel`` used
+                                  to build the confidence and prediction intervals
         conf_intervals:           K rows x 2 columns (lower, upper) confidence intervals
         conf_interval_intercept:  (lower, upper) confidence interval for the intercept
         pi_range:                 prediction intervals above and below, over the range of data
         leverage:                 the hat-matrix diagonal (leverage) for each observation
         influence:                Cook-style influence values for each observation
+
+    On the degenerate early-return path (fewer than three usable observations
+    in ``x`` or ``y``), the returned dictionary is the initialization stub:
+    ``N`` is ``None``, ``x_ssq`` is absent, and ``leverage``, ``influence``,
+    ``pi_range``, ``t_value``, ``fitted_values``, ``residuals``, and the
+    coefficient/interval arrays hold ``np.nan``.
     """
 
     out: dict[str, Any] = {
@@ -301,16 +309,23 @@ def multiple_linear_regression(  # noqa: PLR0913
     Returns a dictionary of outputs. Keys always present::
 
         N:                        number of observations actually used to fit
-        coefficients:             a vector of K coefficients, one for each column in X
+        coefficients:             a vector of K coefficients, one for each column in X;
+                                  shape ``[np.nan]`` on the degenerate/unfitted path
         intercept:                returned if fit_intercept==True
         standard_errors:          a vector of K standard errors, one per column in X
         standard_error_intercept: standard error for the intercept
         R2:                       the R^2 value
+        R2_regression_based:      R^2 computed as ``RegSS / TSS`` (added post-fit)
+        R2_residual_based:        R^2 computed as ``1 - RSS / TSS`` (added post-fit)
+        k:                        the number of model parameters (added post-fit)
         SE:                       the model's standard error
         fitted_values:            the N predicted values, one per row in y
         residuals:                the N residuals
         t_value:                  the t-values for the standard errors
-        conf_intervals:           K rows x 2 columns (lower, upper) confidence intervals
+        conf_intervals:           K rows x 2 columns (lower, upper) confidence intervals;
+                                  shape ``[np.nan, np.nan]`` on the degenerate/unfitted path
+        conf_interval_intercept:  (lower, upper) confidence interval for the intercept
+                                  (added post-fit)
 
     Keys present only for single-feature ``X`` (and only when ``fit_intercept``
     is True and there is enough non-degenerate data)::

From bfe8977c914f2e4a8eafd8df0dc976c70341c201 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:57:32 +0000
Subject: [PATCH 06/15] docs: fix docstring in _mbpls.py

MBPLS: expand the Attributes section to list every fitted attribute
actually set in fit() (n_samples_, n_targets_, n_features_in_,
feature_names_in_, preproc_, y_preproc_, super_hotellings_t2_,
super_vip_, block_spe_, block_hotellings_t2_, block_vip_, and the
r2_x_/r2_y_ family).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/multivariate/_mbpls.py | 37 ++++++++++++++++++++++
 1 file changed, 37 insertions(+)

diff --git a/src/process_improve/multivariate/_mbpls.py b/src/process_improve/multivariate/_mbpls.py
index 7d305f1..5261a93 100644
--- a/src/process_improve/multivariate/_mbpls.py
+++ b/src/process_improve/multivariate/_mbpls.py
@@ -83,6 +83,19 @@ class MBPLS(_HotellingsT2LimitMixin, RegressorMixin, BaseEstimator):
         Ordered list of X-block names (the keys of the input dict).
     block_widths_ : dict[str, int]
         Number of variables in each X-block.
+    n_samples_ : int
+        Number of rows fitted.
+    n_targets_ : int
+        Number of Y columns.
+    n_features_in_ : int
+        Total number of X variables summed across blocks.
+    feature_names_in_ : np.ndarray
+        Concatenated column names, one per feature, in block order.
+    preproc_ : dict[str, MCUVScaler]
+        Per-block preprocessors used to mean-centre and unit-variance
+        scale each X-block.
+    y_preproc_ : MCUVScaler
+        Preprocessor used on Y.
     super_scores_ : pd.DataFrame, shape (n_samples, n_components)
         Super-block (consensus) X-scores ``T``.
     super_y_scores_ : pd.DataFrame, shape (n_samples, n_components)
@@ -91,6 +104,11 @@ class MBPLS(_HotellingsT2LimitMixin, RegressorMixin, BaseEstimator):
         Super-block weights ``w_super``; rows indexed by block name.
     super_y_loadings_ : pd.DataFrame, shape (n_targets, n_components)
         Y-block loadings ``c``.
+    super_hotellings_t2_ : pd.DataFrame, shape (n_samples, n_components)
+        Cumulative Hotelling's T^2 on the super-scores per component.
+    super_vip_ : pd.Series
+        Variable-importance in projection for each X-block, indexed by
+        block name.
     block_scores_ : dict[str, pd.DataFrame]
         Per-block X-scores ``t_b``, each shape ``(n_samples, n_components)``.
     block_weights_ : dict[str, pd.DataFrame]
@@ -99,12 +117,31 @@ class MBPLS(_HotellingsT2LimitMixin, RegressorMixin, BaseEstimator):
     block_loadings_ : dict[str, pd.DataFrame]
         Per-block X-loadings ``p_b`` (used for deflation), each shape
         ``(K_b, n_components)``.
+    block_spe_ : dict[str, pd.DataFrame]
+        Per-block squared prediction error per sample and component.
+    block_hotellings_t2_ : dict[str, pd.DataFrame]
+        Per-block cumulative Hotelling's T^2 per sample and component.
+    block_vip_ : dict[str, pd.Series]
+        Per-block variable-importance in projection, indexed by variable
+        name inside each block.
     predictions_ : pd.DataFrame, shape (n_samples, n_targets)
         In-sample Y predictions on the *original* scale.
     explained_variance_ : np.ndarray, shape (n_components,)
         Variance of the super-score per component (ddof=1).
     scaling_factor_for_super_scores_ : pd.Series
         ``sqrt(explained_variance_)`` per component.
+    r2_x_per_block_cumulative_ : pd.DataFrame, shape (n_blocks, n_components)
+        Cumulative R^2X per block and component.
+    r2_x_per_block_per_component_ : pd.DataFrame, shape (n_blocks, n_components)
+        Incremental R^2X per block and component.
+    r2_x_per_variable_ : dict[str, pd.DataFrame]
+        Cumulative R^2X per variable within each block.
+    r2_y_cumulative_ : pd.Series, shape (n_components,)
+        Cumulative R^2Y per component.
+    r2_y_per_component_ : pd.Series, shape (n_components,)
+        Incremental R^2Y per component.
+    r2_y_per_variable_ : pd.DataFrame, shape (n_targets, n_components)
+        Cumulative R^2Y per Y-variable and component.
     fitting_info_ : dict
         Per-component iteration count and timing.
     has_missing_data_ : bool

From f6b4cdedd4b57a6c3054763e218891f734fdac18 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:58:04 +0000
Subject: [PATCH 07/15] docs: fix docstring in _mbpca.py

MBPCA: expand the Attributes section to list every fitted attribute
individually with its shape/type, replacing the '(as MBPLS)' shorthand.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/multivariate/_mbpca.py | 58 ++++++++++++++++++----
 1 file changed, 48 insertions(+), 10 deletions(-)

diff --git a/src/process_improve/multivariate/_mbpca.py b/src/process_improve/multivariate/_mbpca.py
index 4c17f39..4b7f5ac 100644
--- a/src/process_improve/multivariate/_mbpca.py
+++ b/src/process_improve/multivariate/_mbpca.py
@@ -82,16 +82,54 @@ class MBPCA(_HotellingsT2LimitMixin, TransformerMixin, BaseEstimator):
 
     Attributes (after fitting)
     --------------------------
-    block_names_, block_widths_                       (as MBPLS)
-    super_scores_                                     DataFrame (N x A)
-    super_loadings_                                   DataFrame (B x A)
-    block_scores_, block_loadings_                    dict[str, DataFrame]
-    r2_x_per_block_cumulative_, r2_x_per_block_per_component_
-    r2_x_per_variable_                                dict[str, DataFrame]
-    block_vip_
-    block_spe_, block_hotellings_t2_, super_hotellings_t2_
-    explained_variance_, scaling_factor_for_super_scores_
-    fitting_info_, has_missing_data_, algorithm_
+    block_names_ : list[str]
+        Ordered list of X-block names (the keys of the input dict).
+    block_widths_ : dict[str, int]
+        Number of variables in each X-block.
+    n_samples_ : int
+        Number of rows fitted.
+    n_features_in_ : int
+        Total number of X variables summed across blocks.
+    feature_names_in_ : np.ndarray
+        Concatenated column names, one per feature, in block order.
+    preproc_ : dict[str, MCUVScaler]
+        Per-block preprocessors used to mean-centre and unit-variance
+        scale each X-block.
+    super_scores_ : pd.DataFrame, shape (n_samples, n_components)
+        Super-block (consensus) scores ``T``.
+    super_loadings_ : pd.DataFrame, shape (n_blocks, n_components)
+        Super-block loadings ``p_super``; rows indexed by block name.
+    super_hotellings_t2_ : pd.DataFrame, shape (n_samples, n_components)
+        Cumulative Hotelling's T^2 on the super-scores per component.
+    block_scores_ : dict[str, pd.DataFrame]
+        Per-block scores ``t_b``, each shape ``(n_samples, n_components)``.
+    block_loadings_ : dict[str, pd.DataFrame]
+        Per-block loadings ``p_b``, each shape ``(K_b, n_components)``.
+    block_spe_ : dict[str, pd.DataFrame]
+        Per-block squared prediction error per sample and component.
+    block_hotellings_t2_ : dict[str, pd.DataFrame]
+        Per-block cumulative Hotelling's T^2 per sample and component.
+    block_vip_ : dict[str, pd.Series]
+        Per-block variable-importance in projection, indexed by variable
+        name inside each block.
+    r2_x_per_block_cumulative_ : pd.DataFrame, shape (n_blocks, n_components)
+        Cumulative R^2X per block and component.
+    r2_x_per_block_per_component_ : pd.DataFrame, shape (n_blocks, n_components)
+        Incremental R^2X per block and component.
+    r2_x_per_variable_ : dict[str, pd.DataFrame]
+        Cumulative R^2X per variable within each block.
+    explained_variance_ : np.ndarray, shape (n_components,)
+        Variance of the super-score per component (ddof=1).
+    scaling_factor_for_super_scores_ : pd.Series
+        ``sqrt(explained_variance_)`` per component.
+    fitting_info_ : dict
+        Per-component iteration count and timing.
+    has_missing_data_ : bool
+        Whether any X-block had NaN values.
+    algorithm_ : str
+        The resolved algorithm actually used for the fit. With
+        ``algorithm="auto"``, this is ``"dense"`` for complete data
+        and ``"nipals"`` for NaN-containing data.
 
     Notes
     -----

From dcd1acfd34d16f578fa657f18c94b11b7d975e50 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:58:44 +0000
Subject: [PATCH 08/15] docs: fix docstring in _tpls.py

TPLS: expand the Attributes section from two fields to the full set of
attributes set in fit() (score/loading families, spe/spe_limit nested
dicts, mask arrays, name maps, etc.) and add a note that TPLS does not
follow the sklearn trailing-underscore convention for its fitted
attributes.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/multivariate/_tpls.py | 84 +++++++++++++++++++++--
 1 file changed, 80 insertions(+), 4 deletions(-)

diff --git a/src/process_improve/multivariate/_tpls.py b/src/process_improve/multivariate/_tpls.py
index 93955bd..9b3a820 100644
--- a/src/process_improve/multivariate/_tpls.py
+++ b/src/process_improve/multivariate/_tpls.py
@@ -238,11 +238,87 @@ class TPLS(RegressorMixin, BaseEstimator):
     Attributes
     ----------
     n_samples : int
-        The number of samples (rows) in the training data
-
+        Number of samples (rows) in the training data.
     n_substances : int
-        The number of substances (columns) in the training data, i.e. the number of materials in the F matrix.
-
+        Number of materials (columns) in the F matrix, summed across groups.
+    n_conditions : int
+        Number of process-condition columns summed across Z blocks.
+    n_outputs : int
+        Number of quality-indicator columns summed across Y blocks.
+    is_fitted_ : bool
+        Set to True once ``fit()`` completes.
+    tolerance_ : float
+        Convergence tolerance used by the inner NIPALS loop
+        (``sqrt(np.finfo(float).eps)``).
+    fitting_statistics : dict
+        Per-component ``iterations``, ``convergance_tolerance`` and
+        ``milliseconds`` lists.
+    preproc_ : dict
+        Nested per-block per-group preprocessors, indexed as
+        ``preproc_[block][group]`` where ``block`` is ``"D"``, ``"F"``,
+        ``"Y"`` or ``"Z"``.
+    sums_of_squares_ : list[dict]
+        Per-component per-block sums of squares.
+    r2_frac : list[dict]
+        Per-component per-block cumulative R^2 fractions.
+    feature_importance : dict
+        Per-block per-group variable-importance (populated for ``"D"``
+        and ``"F"``).
+    d_mats, f_mats, z_mats, y_mats : dict[str, np.ndarray]
+        Deflated block matrices (per group for D/F, per Z-block / Y-block
+        for Z/Y).
+    not_na_d, not_na_f, not_na_z, not_na_y : dict[str, np.ndarray]
+        Boolean observed-value masks matching the shapes above.
+    observation_names : pd.Index
+        Row index shared across the F/Z/Y blocks.
+    property_names, material_names : dict[str, list[str]]
+        Column and row names of each D group.
+    condition_names, quality_names : dict[str, list[str]]
+        Column names of each Z-block and Y-block respectively.
+    t_scores_super : pd.DataFrame, shape (n_samples, n_components)
+        Super-block scores ``T``.
+    r_loadings_f : dict[str, pd.DataFrame]
+        F-block weights per group.
+    w_loadings_z : dict[str, pd.DataFrame]
+        Z-block weights per Z-block.
+    w_loadings_super : pd.DataFrame
+        Super-block weights, rows ``["Z", "F"]`` (or just ``["F"]`` when
+        there are no process conditions).
+    s_loadings_d, v_loadings_d : dict[str, pd.DataFrame]
+        D-block scores and loadings per group.
+    p_loadings_f : dict[str, pd.DataFrame]
+        F-block loadings per group (used for deflation).
+    p_loadings_z : dict[str, pd.DataFrame]
+        Z-block loadings per Z-block.
+    q_loadings_y : dict[str, pd.DataFrame]
+        Y-block loadings per Y-block.
+    hat_ : dict[str, pd.DataFrame]
+        Per-Y-block predictions on the preprocessed (centred / scaled) scale.
+    hat : dict[str, pd.DataFrame]
+        Per-Y-block in-sample predictions on the *original* scale.
+    spe : dict[str, dict[str, pd.DataFrame or pd.Series]]
+        Nested ``spe[block][group]`` squared prediction error tables.
+    spe_limit : dict[str, dict[str, Callable]]
+        Nested ``spe_limit[block][group]`` callables. Each is a
+        :func:`functools.partial` over :func:`spe_calculation` bound to
+        the block's SPE array; call it with a confidence level to obtain
+        the SPE limit. This is a deliberate divergence from the flat
+        ``spe_limit`` method on PCA / PLS.
+    hotellings_t2 : pd.DataFrame, shape (n_samples, n_components)
+        Cumulative Hotelling's T^2 per super-component.
+    scaling_factor_for_scores : pd.Series
+        Per-component scaling factor used by ellipse / T^2 helpers.
+
+    .. note::
+       TPLS deliberately does **not** follow the sklearn trailing-
+       underscore convention for its fitted attributes. Names such as
+       ``t_scores_super``, ``spe``, ``hotellings_t2``, and the
+       ``*_loadings_*`` family are written without the trailing ``_``
+       to keep the chemometrics symbol names readable. Attributes
+       that are set in ``__init__`` and refined during ``fit`` (for
+       example ``is_fitted_``, ``preproc_``, ``tolerance_``,
+       ``required_blocks_``, ``required_inputs_``) do carry the
+       underscore.
 
     Example
     -------

From 80bbac158dc9d73629f19456817839a7b119d126 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:59:09 +0000
Subject: [PATCH 09/15] docs: fix docstring in _opls.py

OPLS.spe_: clarify that the returned DataFrame's per-column values are
broadcasts of the final-component SPE, not a per-component progression
like the PCA/PLS spe_.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/multivariate/_opls.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/process_improve/multivariate/_opls.py b/src/process_improve/multivariate/_opls.py
index bfcdccc..aa84bc1 100644
--- a/src/process_improve/multivariate/_opls.py
+++ b/src/process_improve/multivariate/_opls.py
@@ -99,7 +99,13 @@ class OPLS(_LatentVariableModel, RegressorMixin, TransformerMixin, BaseEstimator
         The orthogonal-signal-corrected X (``X`` with the orthogonal variation
         removed), on the scaled fitting scale.
     spe_ : pd.DataFrame
-        Per-row SPE after reconstructing X from all components.
+        Per-row SPE after reconstructing X from all components. Unlike
+        the PCA / PLS ``spe_`` (which stores a per-component
+        progression, one column per component), the OPLS ``spe_`` only
+        holds the final-component SPE and broadcasts that single
+        column across every ``t_predictive`` / ``t_orthogonal_i``
+        column, so all columns are identical. This keeps the shape
+        aligned with :attr:`scores_` for the inherited SPE plots.
     hotellings_t2_ : pd.DataFrame
         Cumulative Hotelling's T2 over the combined score space.
 

From b68b468ba4f87de84fe2a25b6c39f4877f93bdd3 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 07:59:30 +0000
Subject: [PATCH 10/15] docs: fix docstring in _limits.py

spe_calculation: document the SEC-21/#270 fallback: when variance or
centre of the squared SPE is at or below epsqrt, the Jackson-Mudholkar
chi-square approximation degenerates and the limit falls back to
sqrt(center_spe).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/multivariate/_limits.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/process_improve/multivariate/_limits.py b/src/process_improve/multivariate/_limits.py
index 6af9a11..3303627 100644
--- a/src/process_improve/multivariate/_limits.py
+++ b/src/process_improve/multivariate/_limits.py
@@ -91,6 +91,16 @@ def spe_calculation(spe_values: np.ndarray, conf_level: float = 0.95) -> float:
     float
         The limit, above which we judge observations in the model to have a different correlation
         structure than those values which were used to build the model.
+
+    Notes
+    -----
+    When either the variance or the centre of the squared SPE values is
+    at or below ``epsqrt`` (a perfect-fit training set where ``A == K``,
+    or an all-equal SPE column), the Jackson-Mudholkar chi-square
+    approximation degenerates. In that case the limit falls back to
+    ``sqrt(center_spe)``: there is no spread to bound, so any value
+    above the centre is by construction out of family. See SEC-21
+    (#270), sub-item 3.
     """
     if not 0.0 < conf_level < 1.0:
         raise ValueError(f"conf_level must lie in (0, 1); got {conf_level}.")

From f81f8f726532dd44052fd57732661b4be275d231 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 08:00:34 +0000
Subject: [PATCH 11/15] docs: fix docstring in metrics.py

t_value / t_value_cdf: replace undefined 'v' in the doctest examples
with 'v=10' and correct the -inf / +inf renderings (scipy returns
-inf / +inf at p in {0, 1}, not NaN, and Python prints them lowercase).

ttest_paired already documents that its 'Standard deviation' key holds
the standard error of the mean difference, so no change is needed
there.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/univariate/metrics.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/process_improve/univariate/metrics.py b/src/process_improve/univariate/metrics.py
index 8033374..5fbdebc 100644
--- a/src/process_improve/univariate/metrics.py
+++ b/src/process_improve/univariate/metrics.py
@@ -27,18 +27,18 @@ def t_value(p: float, v: float) -> float:
     Since the cumulative distribution passes symmetrically through the x-axis at 0.0 for any
     number of degrees of freedom
 
-    >>> t_value(0.5, v)
+    >>> t_value(0.5, v=10)
     0.0
 
     Zero fractional area under the curve is always at :math:`-\infty`:
 
-    >>> t_value(0.0, v)
-    -Inf
+    >>> t_value(0.0, v=10)
+    -inf
 
     100% fractional area is always at :math:`+\infty`:
 
-    >>> t_value(1.0, v)
-    +Inf
+    >>> t_value(1.0, v=10)
+    inf
 
     See also
     --------
@@ -57,17 +57,17 @@ def t_value_cdf(z: float, v: float) -> float:
     The cumulative distribution is symmetric through the x-axis at 0.0 for any number of degrees
     of freedom, so half of the area lies below zero:
 
-    >>> t_value_cdf(0.0, v)
+    >>> t_value_cdf(0.0, v=10)
     0.5
 
     Zero fractional area under the curve is at :math:`-\infty`:
 
-    >>> t_value_cdf(-np.inf, v)
+    >>> t_value_cdf(-np.inf, v=10)
     0.0
 
     100% fractional area is at :math:`+\infty`:
 
-    >>> t_value_cdf(np.inf, v)
+    >>> t_value_cdf(np.inf, v=10)
     1.0
 
     See also

From 7b1777c4a8ba6f860f78afdaad83ee093dcafab9 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 08:01:01 +0000
Subject: [PATCH 12/15] docs: fix docstring in _elbow_peak.py

find_elbow_point: add an explicit Returns section that covers both
the -1 (all-missing) and NaN (no consensus intersection) return paths.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/bivariate/_elbow_peak.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/process_improve/bivariate/_elbow_peak.py b/src/process_improve/bivariate/_elbow_peak.py
index 85cf25f..3d7ca14 100644
--- a/src/process_improve/bivariate/_elbow_peak.py
+++ b/src/process_improve/bivariate/_elbow_peak.py
@@ -16,7 +16,7 @@ def find_elbow_point(x: np.ndarray, y: np.ndarray, max_iter: int = 41) -> int |
     Find the elbow point when plotting numeric entries in `x` vs numeric values in list `y`.
 
     Return the index into the vectors `x` and `y` [the vectors must have the same length], where
-    the elbow point occurs. Returns -1 if every value in `x` or `y` is missing.
+    the elbow point occurs.
 
     Using a robust linear fit, sorts the samples in X (independent variable)
     and takes the first 5 samples from the left, and the last 5 from the right,
@@ -32,6 +32,14 @@ def find_elbow_point(x: np.ndarray, y: np.ndarray, max_iter: int = 41) -> int |
     Will probably not work well on few data points. If so, try fitting a spline
     to the raw data and then repeat with the interpolated data.
 
+    Returns
+    -------
+    int or float
+        The 0-based index of the elbow point in the (sorted) vectors.
+        Returns ``-1`` if every value in ``x`` or ``y`` is missing.
+        Returns ``np.nan`` when the intersection sweep produced only
+        NaNs (for example, when every candidate line pair was
+        near-parallel), so no consensus intersection could be formed.
     """
     start = 5
     # assert divmod(max_iter, 2)[1]  # must be odd number; to ensure we calculate the median later

From bc18ce3c740dd45c446dbbb16ec03f9d0edfdc8b Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 08:01:26 +0000
Subject: [PATCH 13/15] docs: fix docstring in analysis.py

analyze_experiment: enumerate every key that the always-present
model_summary dict actually carries (formula, n_obs, n_terms,
model_rank, rank_deficient, df_model, df_residual, mse_residual),
in addition to the R^2 and adequate-precision entries that were
already documented.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/experiments/analysis.py | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/src/process_improve/experiments/analysis.py b/src/process_improve/experiments/analysis.py
index afcc7ad..60b5b84 100644
--- a/src/process_improve/experiments/analysis.py
+++ b/src/process_improve/experiments/analysis.py
@@ -188,7 +188,22 @@ def analyze_experiment(  # noqa: PLR0912, PLR0913, PLR0915, C901
     -------
     dict[str, Any]
         Results keyed by analysis type.  Always includes ``"model_summary"``
-        with R², adj-R², pred-R², and adequate precision.
+        with the keys:
+
+        - ``formula`` - the resolved patsy formula that was fitted.
+        - ``r_squared`` - R^2 of the fit.
+        - ``r_squared_adj`` - adjusted R^2.
+        - ``r_squared_pred`` - prediction R^2 (leave-one-out style).
+        - ``adequate_precision`` - signal-to-noise ratio (>= 4 is
+          considered adequate).
+        - ``n_obs`` - number of observations used to fit.
+        - ``n_terms`` - number of columns in the model matrix.
+        - ``model_rank`` - numerical rank of the model matrix; less
+          than ``n_terms`` implies aliasing / rank deficiency.
+        - ``rank_deficient`` - ``True`` if ``model_rank < n_terms``.
+        - ``df_model`` - model degrees of freedom.
+        - ``df_residual`` - residual degrees of freedom.
+        - ``mse_residual`` - mean squared error of the residuals.
 
     Examples
     --------

From bd4d2a20330849cf1b7ec371b19c3bb0f1b5dd32 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 08:01:45 +0000
Subject: [PATCH 14/15] docs: fix docstring in control_charts.py

ControlChart.__init__: mark 'cusum' as an unimplemented future variant;
only 'hw' and 'xbar.no.subgroup' are dispatched by calculate_limits.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 src/process_improve/monitoring/control_charts.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/process_improve/monitoring/control_charts.py b/src/process_improve/monitoring/control_charts.py
index 78d902a..9d63a68 100644
--- a/src/process_improve/monitoring/control_charts.py
+++ b/src/process_improve/monitoring/control_charts.py
@@ -59,7 +59,10 @@ def __init__(self, style: str = "robust", variant: str = "HW") -> None:
                 'xbar.no.subgroup' [Shewhart chart, with no subgroups]. In other words, each
                 observation is independently plotted on the control chart.
 
-                'cusum' (CUmulative SUM) chart, which uses all the history of the chart.
+                A pure 'cusum' (CUmulative SUM) chart is a planned future variant but
+                is not currently implemented; the Holt-Winters ('hw') default already
+                blends CUSUM-style infinite history with Shewhart-style
+                instantaneous behaviour via its lambda parameters.
         """
         self.style = style.strip()
         self.variant = variant.strip().lower()

From 3b6f868dd2bc06486ffa418b3785db1261a67e1a Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 14 Aug 2026 08:03:16 +0000
Subject: [PATCH 15/15] chore: bump version to 1.66.2 and log docstring audit

Docs-only PATCH bump. Sync CITATION.cff to the same version and today's
date; add a Documentation entry under [Unreleased] in CHANGELOG.md
summarising the 17 docstring corrections in this PR.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0132CNPVroNNG679mikPzqwc
---
 CHANGELOG.md   | 25 +++++++++++++++++++++++++
 CITATION.cff   |  4 ++--
 pyproject.toml |  2 +-
 3 files changed, 28 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c8f0bc1..6adb378 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -92,6 +92,31 @@ those changes.
 
 - The `trade_off_table` MCP tool is untouched: same tool name, same `runs` and
   `factors` schema keys, same output. Only its internal call site moved.
+
+### Documentation
+
+- Docstring audit correcting drift between NumPy-style docstrings and the
+  actual runtime behaviour, no runtime changes:
+  `dispatch_ccd` and `generate_design` now warn that a numeric `alpha`
+  is only honored for the fractional-cube CCD path;
+  `PCA.select_n_components` documents `return_consensus` and the four
+  extra keys it exposes; `PLS.select_n_components` completes the
+  truncated `selection_mode` sentence; `MBPLS`, `MBPCA`, and `TPLS`
+  now list every fitted attribute set in `fit()` (TPLS also notes that
+  it deliberately does not follow the sklearn trailing-underscore
+  convention); `OPLS.spe_` clarifies that per-column values are
+  broadcasts of the final-component SPE; `spe_calculation` documents
+  the SEC-21 / #270 low-variance fallback; `robust_regression` and
+  `multiple_linear_regression` describe their full dict outputs
+  (including the degenerate/unfitted-path shapes and the
+  `R2_regression_based` / `R2_residual_based` / `k` /
+  `conf_interval_intercept` keys); `t_value` and `t_value_cdf` use a
+  concrete `v=10` in doctest examples with correct `-inf` / `inf`
+  renderings; `find_elbow_point` documents the secondary NaN return
+  path; `analyze_experiment` enumerates every key on the always-
+  present `model_summary` dict; `ControlChart.__init__` marks
+  `'cusum'` as an unimplemented future variant.
+
 ## [1.66.1] - 2026-08-09
 
 ### Changed
diff --git a/CITATION.cff b/CITATION.cff
index e6e4e1d..d91c294 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -12,8 +12,8 @@ authors:
 repository-code: "https://github.com/kgdunn/process-improve"
 url: "https://kgdunn.github.io/process-improve/"
 license: MIT
-version: 1.66.1
-date-released: "2026-08-09"
+version: 1.66.2
+date-released: "2026-08-14"
 keywords:
   - chemometrics
   - multivariate analysis
diff --git a/pyproject.toml b/pyproject.toml
index f73ad03..bce2b63 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "process-improve"
-version = "1.66.1"
+version = "1.66.2"
 description = 'Designed Experiments; Latent Variables (PCA, PLS, multivariate methods with missing data); Process Monitoring; Batch data analysis.'
 readme = "README.md"
 license = "MIT"