docs: correct docstring/implementation drift across multivariate, regression, experiments, and monitoring

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

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

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"

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

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
     --------

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"``

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

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()

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}.")

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
     -----

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

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.
 

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).

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.