Audit: test validity

plsdo/cli.py

@@ -81,8 +81,9 @@ def pls_main(argv=None):
         default=1.96,
         type=float,
         help=(
-            "Plot loading bars only for features with |bootstrap ratio| > "
-            "THRESHOLD (default: 1.96). Does not affect CSV outputs."
+            "Bootstrap-ratio reliability threshold (default: 1.96). Controls "
+            "both which latent variables survive and which loading bars are "
+            "plotted. Per-component CSV outputs are written regardless."
         ),
     )
     run_common.add_argument(

plsdo/core.py

@@ -49,6 +49,7 @@ def fit(self):
         """Run PLS: cross-covariance, SVD, loadings, and subject scores."""
         self.xcorr = self.X.T @ self.Y / (self.n_subjects - 1)
         self._decompose()
+        self._fix_component_signs()
         self.u_loadings = self.u * self.s[np.newaxis, :]
         self.vt_loadings = self.s[:, np.newaxis] * self.vt
         self.x_scores = self.X @ self.u
@@ -63,6 +64,22 @@ def _decompose(self):
         """
         self.u, self.s, self.vt = np.linalg.svd(self.xcorr, full_matrices=False)
 
+    def _fix_component_signs(self):
+        """Pin each component's arbitrary global sign deterministically.
+
+        A PLS component's sign is not scientifically meaningful, but
+        ``np.linalg.svd`` can return a different one across BLAS builds, which
+        would flip loadings, scores, and bootstrap ratios machine-to-machine.
+        Flip each component so its largest-magnitude X loading is positive,
+        making all package outputs reproducible across machines. ``u`` and
+        ``vt`` share a component's sign, so both are flipped together.
+        """
+        max_idx = np.argmax(np.abs(self.u), axis=0)
+        signs = np.sign(self.u[max_idx, np.arange(self.u.shape[1])])
+        signs[signs == 0] = 1.0
+        self.u = self.u * signs
+        self.vt = self.vt * signs[:, np.newaxis]
+
     def _check_fitted(self):
         """Raise if fit() has not been called."""
         if not self._fitted:
@@ -137,14 +154,11 @@ def bootstrap(self, n_bootstraps: int = 10000) -> None:
             aligned_u_load = boot_u_load @ Q
             aligned_vt_load = Q.T @ boot_vt_load
 
-            # Sign correction
-            signs = np.sign(
-                np.sum(aligned_vt_load * self.vt_loadings, axis=1, keepdims=True)
-            )
-            signs[signs == 0] = 1.0
-
-            u_distribution.append(aligned_u_load * signs.T)
-            vt_distribution.append(aligned_vt_load * signs)
+            # No separate sign correction: orthogonal_procrustes returns an
+            # unconstrained orthogonal matrix (reflections allowed), so the
+            # alignment above already resolves each component's arbitrary sign.
+            u_distribution.append(aligned_u_load)
+            vt_distribution.append(aligned_vt_load)
 
         self.u_se = np.std(np.stack(u_distribution, axis=2), axis=2, ddof=1)
         self.vt_se = np.std(np.stack(vt_distribution, axis=2), axis=2, ddof=1)
@@ -172,7 +186,7 @@ def filter_lvs(self, bsr_threshold: float = 1.96) -> None:
         if not self._bootstrapped:
             raise RuntimeError("Call .bootstrap() before .filter_lvs().")
 
-        significant = self.p_values < 0.05
+        significant = self.significant_lvs
 
         # Check if any feature exceeds threshold on X side
         x_reliable = np.any(np.abs(self.u_bootstrap_ratios) > bsr_threshold, axis=0)

plsdo/cross_validate.py

@@ -2,14 +2,20 @@
 
 import numpy as np
 import pandas as pd
-from sklearn.cross_decomposition import PLSRegression
-from sklearn.metrics import (
-    accuracy_score,
-    balanced_accuracy_score,
-    confusion_matrix,
-)
-from sklearn.model_selection import RepeatedStratifiedKFold
-from sklearn.preprocessing import StandardScaler
+try:
+    from sklearn.cross_decomposition import PLSRegression
+    from sklearn.metrics import (
+        accuracy_score,
+        balanced_accuracy_score,
+        confusion_matrix,
+    )
+    from sklearn.model_selection import RepeatedStratifiedKFold
+    from sklearn.preprocessing import StandardScaler
+except ImportError as exc:
+    raise ImportError(
+        "scikit-learn is required for cross-validation but is not installed. "
+        "Install the optional dependency with: pip install 'plsdo[cv]'"
+    ) from exc
 
 from plsdo.io import corrected_pvalue
 
@@ -45,6 +51,10 @@ def run_cv(
         mean_accuracy, mean_balanced_accuracy, fold_results (DataFrame),
         true_labels, pred_labels, confusion_matrix
     """
+    # Deliberate X/Y flip, opposite to `plsdo discriminatory`: the continuous
+    # data ``X`` is the *predictor* and the dummy-coded groups are the
+    # *target*, so PLSRegression.predict yields predicted group scores that
+    # argmax into class labels.
     n_groups = len(np.unique(labels))
     Y_dummy = np.eye(n_groups)[labels]
 

plsdo/pipeline.py

@@ -39,13 +39,35 @@
 logger = logging.getLogger("plsdo")
 
 
+def _dependency_versions() -> dict[str, str]:
+    """Versions of the numerical dependencies, for reproducibility.
+
+    scikit-learn is optional (only the cross-validate path needs it), so it is
+    reported as "not installed" when absent rather than failing.
+    """
+    import numpy
+    import scipy
+
+    versions = {"numpy": numpy.__version__, "scipy": scipy.__version__}
+    try:
+        import sklearn
+
+        versions["scikit-learn"] = sklearn.__version__
+    except ImportError:
+        versions["scikit-learn"] = "not installed"
+    return versions
+
+
 def _write_log(output_dir: Path, params: dict, notes: list[str] | None = None) -> None:
     """Write a log.txt with run parameters and optional trailing notes."""
     log_path = output_dir / "log.txt"
     with open(log_path, "w") as f:
         f.write("PLS analysis log\n")
         f.write(f"Version: {__version__}\n")
         f.write(f"Timestamp: {datetime.now().isoformat()}\n")
+        f.write("\nLibrary versions:\n")
+        for name, version in _dependency_versions().items():
+            f.write(f"  {name}: {version}\n")
         f.write("\nParameters:\n")
         for k, v in params.items():
             f.write(f"  {k}: {v}\n")
@@ -131,9 +153,12 @@ def run_pipeline(
     all_plots : bool
         If True, generate additional diagnostic plots.
     bsr_threshold : float
-        Plot loading bars only for features with |bootstrap ratio|
-        exceeding this threshold. Default 1.96 (≈ 95% CI under the
-        standard-normal approximation). CSV outputs are unaffected.
+        Bootstrap-ratio reliability threshold. Default 1.96 (≈ 95% CI under
+        the standard-normal approximation). Controls both which latent
+        variables survive ``filter_lvs`` (a surviving LV needs at least one
+        feature with |bootstrap ratio| > threshold on each side) and which
+        loading bars are plotted. The loading and bootstrap-ratio CSVs are
+        written for every component regardless.
     verbose_feature_limit : int, optional
         Maximum number of features before verbose plots (except scree)
         are skipped. Defaults to ``VERBOSE_FEATURE_LIMIT`` (100).
@@ -211,7 +236,7 @@ def run_pipeline(
     model.fit()
     model.permutation_test(n_perms=n_perms)
     model.bootstrap(n_bootstraps=n_bootstraps)
-    model.filter_lvs()
+    model.filter_lvs(bsr_threshold=bsr_threshold)
 
     # --- Save data CSVs ---
     _save_csv(
@@ -544,6 +569,9 @@ def cross_validate_pipeline(
 
     # --- Run CV ---
     logger.info("Running %d-fold CV with %d repeats...", n_folds, n_repeats)
+    # X/Y flip: the continuous Y-matrix is the CV *predictor* and the
+    # demographic groups are the classification *target* (opposite to the
+    # discriminatory run_pipeline convention).
     cv_result = run_cv(
         Y,
         labels,

plsdo/plotting.py

@@ -366,6 +366,9 @@ def _box_strip_facet(
     )
     if rotate_xticklabels:
         for ax in g.axes.flat:
+            # Pin the tick locator before relabelling: set_xticklabels alone
+            # warns when the number of ticks is not fixed first.
+            ax.set_xticks(ax.get_xticks())
             ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha="right")
     g.add_legend()
     _finalise(g, out_path, dpi)

pyproject.toml

@@ -61,6 +61,13 @@ path = "plsdo/__init__.py"
 testpaths = ["tests"]
 filterwarnings = [
     "error::UserWarning:seaborn",
+    # seaborn 0.13.2 (the latest release, Jan 2024) calls matplotlib's bxp()
+    # with vert=, which matplotlib >= 3.10 deprecates in favour of orientation=.
+    # The fix is merged on seaborn's main branch (PR #3820) but unreleased, and
+    # seaborn has not cut a release in over two years — so there is no version
+    # to bump to and no firm timeline. Ignore this third-party warning; remove
+    # this line if/when seaborn next releases.
+    "ignore:vert.*deprecated:PendingDeprecationWarning",
 ]
 
 [tool.coverage.run]

tests/test_cli.py

@@ -68,7 +68,14 @@ def test_discriminatory_without_group_col_errors(self, data_dir, tmp_path, capsy
         captured = capsys.readouterr()
         assert "requires --group-col" in captured.err.lower()
 
-    def test_corr_alias_runs(self, data_dir, tmp_path):
+    def test_corr_alias_dispatches_correlational(self, data_dir, tmp_path, monkeypatch):
+        # The alias only needs to resolve to the correlational subcommand;
+        # the full run is exercised elsewhere. Stub the pipeline to keep this a
+        # fast dispatch check.
+        captured = {}
+        monkeypatch.setattr(
+            "plsdo.pipeline.run_pipeline", lambda **kw: captured.update(kw)
+        )
         pls_main(
             [
                 "corr",
@@ -80,18 +87,17 @@ def test_corr_alias_runs(self, data_dir, tmp_path):
                 str(data_dir / "demographics.csv"),
                 "--output",
                 str(tmp_path / "out"),
-                "--n-perms",
-                "10",
-                "--n-bootstraps",
-                "10",
-                "--subject-id",
-                "subject_id",
             ]
         )
-        assert (tmp_path / "out" / "data").exists()
+        assert captured["method"] == "correlational"
 
-    def test_discrim_alias_runs(self, data_dir, tmp_path):
-        out = tmp_path / "out_discrim"
+    def test_discrim_alias_dispatches_discriminatory(
+        self, data_dir, tmp_path, monkeypatch
+    ):
+        captured = {}
+        monkeypatch.setattr(
+            "plsdo.pipeline.run_pipeline", lambda **kw: captured.update(kw)
+        )
         pls_main(
             [
                 "discrim",
@@ -101,17 +107,11 @@ def test_discrim_alias_runs(self, data_dir, tmp_path):
                 str(data_dir / "demographics.csv"),
                 "--group-col",
                 "group",
-                "--subject-id",
-                "subject_id",
                 "--output",
-                str(out),
-                "--n-perms",
-                "10",
-                "--n-bootstraps",
-                "10",
+                str(tmp_path / "out_discrim"),
             ]
         )
-        assert (out / "data").exists()
+        assert captured["method"] == "discriminatory"
 
     def test_group_col_and_groups_mutually_exclusive(
         self,
@@ -312,8 +312,14 @@ def test_runs_successfully(self, data_dir, tmp_path):
         assert (out / "data").exists()
         assert (out / "log.txt").exists()
 
-    def test_cv_alias_runs(self, data_dir, tmp_path):
-        out = tmp_path / "cv_alias"
+    def test_cv_alias_dispatches_cross_validate(self, data_dir, tmp_path, monkeypatch):
+        # The alias only needs to resolve to the cross-validate subcommand; the
+        # full run is exercised by test_runs_successfully. Stub the pipeline.
+        called = {}
+        monkeypatch.setattr(
+            "plsdo.pipeline.cross_validate_pipeline",
+            lambda **kw: called.update(kw, dispatched=True),
+        )
         pls_main(
             [
                 "cv",
@@ -323,19 +329,12 @@ def test_cv_alias_runs(self, data_dir, tmp_path):
                 str(data_dir / "demographics.csv"),
                 "--group-col",
                 "group",
-                "--subject-id",
-                "subject_id",
                 "--output",
-                str(out),
-                "--n-folds",
-                "3",
-                "--n-repeats",
-                "2",
-                "--n-permutations",
-                "10",
+                str(tmp_path / "cv_alias"),
             ]
         )
-        assert (out / "data").exists()
+        assert called["dispatched"] is True
+        assert called["group_col"] == "group"
 
     def test_accepts_groups_yaml(self, data_dir, tmp_path):
         out = tmp_path / "cv_yaml"