PLS-DA — Partial Least Squares Discriminant Analysis

Task	Classification
Method key (analysis_type)	plsda
Prediction on new data	Yes — see Prediction
Library	Own class PLSDA in analysis.py, built on sklearn.cross_decomposition.PLSRegression, sklearn.decomposition.PCA and scipy.stats.chi2

When to use

PLS-DA is the standard chemometric classifier:

• many correlated features (colour channels, time points), possibly more features than samples;
• you want a model you can look at: scores, loadings, feature weights;
• with the soft rule, you want the model to say "this sample belongs to none of the known classes" instead of forcing a class.

How it works

The implementation follows Pomerantsev & Rodionova (2018) and the PyChemAuth package.

1. Encode the classes. The labels of $K$ classes become a one-hot matrix $Y$ ($n \times K$). $X$ is centred (and scaled, if Scale Features is on); $Y$ is centred.
2. PLS2 regression. PLSRegression with $a$ components finds latent variables that explain the covariance between $X$ and $Y$, and predicts $\hat{Y} = X B$.
3. PCA of the predictions. The rows of $\hat{Y}$ (back on the one-hot scale) are projected by PCA onto $K - 1$ components, giving scores $T$. The class centres $c_k$ are the projections of the ideal one-hot vectors (the rows of the identity matrix).
4. Distances. For each sample $t$ and class $k$, the squared Mahalanobis distance $$d_k(t) = (t - c_k)^\top S^{-1} (t - c_k)$$ is computed, with a covariance matrix $S$ that depends on the rule.
5. Decision rule:

Rule	style	Covariance $S$	Assignment
Hard	hard	One matrix for all classes: the variance of $T$ along each PCA component (a diagonal matrix)	The class with the smallest $d_k$. Every sample gets exactly one class.
Soft	soft	One matrix per class: the scatter of the class samples around their centre	All classes with $d_k < \chi^2_{1-\alpha}(K-1)$. If there are several, the closest one is reported; if there are none, the sample is not assigned (UNKNOWN).

The difference between the two rules is explained in Hard vs soft PLS-DA.

Probabilities (predict_proba) are derived from the distances: $p_k = e^{-d_k/2}$, normalised to sum to 1 for the hard rule, and divided by $\sqrt{\det(2\pi S_k)}$ (capped at 1) for the soft rule. Soft probabilities therefore do not sum to 1.

Parameters

Defaults as set by the analysis dialog.

Parameter	UI label	Default	Notes
n_components	Number of Components	$\min(p, 7)$, $p$ = number of feature columns	Number of PLS components $a$. Reduced to $\min(n_\text{train} - 1, p)$ if larger. Used only when the search is off.
style	Classification Style	hard	hard or soft
search_optimal_n_components	Search optimal n_components	on	See Choosing the number of components
search_optimal_n_components.metric	Optimization Metric	accuracy	Accuracy, F1 Macro, F1 Weighted, Precision Macro, Recall Macro; for two classes also Sensitivity, Specificity, Balanced Accuracy, R2
search_optimal_n_components.min_n	Min Components	2
search_optimal_n_components.max_n	Max Components	$\min(p, 10)$
alpha	Significance Level (α)	0.05	Soft rule only: the share of samples of a class that may fall outside its acceptance area
gamma	Outlier Threshold (γ)	0.01	Stored with the model, not used by the current implementation
scale_x	Scale Features	on	Divide each feature by its standard deviation after centring
feature_selection	Enable Feature Selection	off	See Feature selection.
average_replicates	Average technical replicates	off	See Averaging replicates.

Choosing the number of components

With Search optimal n_components on (the default), the analysis runs the full cross-validation for every number of components from Min Components to Max Components and keeps the one with the best value of the Optimization Metric. The number that was chosen and the score of every candidate are saved in optimal_n_components_search.

warning

The reported CV metrics are the ones of the best candidate, measured on the same folds that were used to choose it. They are therefore slightly optimistic, more so the wider the search range. To check the result, run a permutation test: it uses the chosen number of components.

F1 Weighted is offered in the dialog, but the cross-validation does not compute it yet: with this metric the analysis fails. Use F1 Macro instead. Balanced accuracy is the mean of CV sensitivity and CV specificity.

Preprocessing

• $X$ is centred, and scaled to unit variance with Scale Features on. Scaling is fitted on the training rows of each fold only.
• Before the model: replicate averaging and feature selection, if they are on. See Order of steps.

Results and metrics

The common classification metrics are described in Metrics. With the soft rule, a sample that is not assigned (UNKNOWN) counts as an error in accuracy. For binary problems AUC uses the probability of the positive class, and two extra metrics are reported:

• R² (training) and Q² (cv_r2, cross-validated): the class labels are coded as 0 and 1 and compared with the predicted class codes (an unassigned sample counts as 0.5): $$Q^2 = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2} .$$

PLS-DA-specific fields:

Field	Meaning
transformed_data	Scores $T$ of the samples in the $(K-1)$-dimensional PCA space of the predictions (final model)
loadings	PLS weights x_weights_ (one row per feature, one column per PLS component), or correlation loadings if loading_type = 'correlation'
feature_importance	Sum of the absolute PLS weights of each feature over all components
prediction_probabilities	See Probabilities above
model_info	Parameters actually used, including the chosen n_components
optimal_n_components_search	best_n_components, metric and the list of (components, score) pairs, if the search was on

Visualizations

• Score plot and loading plot
• Confusion matrix and ROC curve
• Permutation plot, with the extra metric R² / Q²
• Time series and heatmap of the input data

For two classes the PCA space of the predictions has $K - 1 = 1$ dimension: the score plot has one meaningful axis.

Source code

The model is the class PLSDA in chrometrica/analysis/analysis.py. Fitting:

analysis.py · PLSDA.fit() · lines 666–761

def fit(self, X, y):
    """
    Fit PLS-DA model to training data following pychemauth algorithm.
    """
    X = np.array(X, dtype=np.float64)
    y = np.array(y)

    # Step 1: Preprocess data (following pychemauth exactly)
    self._X_train_scaled = self._stable_scale(X, fit=True)
    self._y_train_encoded = self._encode_y(y, fit=True)

    # Step 2: Perform PLS2 regression
    k = len(self.encoder.categories_[0])  # number of classes

    # Bound n_components properly
    upper_bound = min(self._X_train_scaled.shape[0] - 1, self._X_train_scaled.shape[1])
    self.n_components = min(self.n_components, upper_bound)

    self.pls = PLSRegression(
        n_components=self.n_components,
        max_iter=10000,
        tol=1e-9,
        scale=False  # Already scaled
    )
    self.pls.fit(self._X_train_scaled, self._y_train_encoded)

    # Step 3: Get predictions and perform PCA (CRITICAL: follow pychemauth exactly)
    # Get PLS predictions and inverse transform Y scaling to get back to one-hot scale
    y_hat_train = self.y_scaler.inverse_transform(
        self.pls.predict(self._X_train_scaled)
    )

    # Perform PCA on y_hat_train with k-1 components
    n_pca_components = k - 1
    self.pca = PCA(n_components=n_pca_components, random_state=self.random_state)
    T_train = self.pca.fit_transform(y_hat_train)

    # Step 4: Compute class centers in PCA space (CRITICAL: pychemauth method)
    # Class centers are projections of identity matrix (one-hot vectors)
    identity_matrix = np.eye(k)
    self.class_centers = self.pca.transform(identity_matrix)

    # Step 5: Store class masks for training data
    self.class_masks = {}
    for i in range(k):
        class_label = self.encoder.categories_[0][i]
        self.class_masks[i] = (y == class_label)

    # Step 6: Compute covariance matrices based on classification style
    if self.style == 'hard':
        # For hard classification: L = diag(explained_variance)
        # This matches pychemauth: L = np.cov(T_train.T) which equals diag(explained_variance)
        self.cov_matrix = np.eye(len(self.pca.explained_variance_)) * self.pca.explained_variance_

        # Handle 1D case (binary classification)
        if self.cov_matrix.ndim == 0:
            self.cov_matrix = np.array([[self.cov_matrix]])

    else:  # soft classification
        # Compute within-class scatter matrices (following pychemauth exactly)
        self.within_class_cov = {}

        for i in range(k):
            mask = self.class_masks[i]
            if np.sum(mask) > 0:
                # Get class samples in T space and center around class center
                t_class = T_train[mask] - self.class_centers[i]

                # Compute scatter matrix as in pychemauth
                S_i = np.zeros((T_train.shape[1], T_train.shape[1]), dtype=np.float64)
                for j in range(t_class.shape[0]):
                    t_vec = t_class[j, :].reshape(-1, 1)
                    S_i += np.dot(t_vec, t_vec.T)
                S_i /= t_class.shape[0]  # Divide by number of samples

                # Check if positive definite and regularize if needed
                try:
                    np.linalg.cholesky(S_i)
                    self.within_class_cov[i] = S_i
                except np.linalg.LinAlgError:
                    # Add small regularization to diagonal
                    eigenvals = np.linalg.eigvals(S_i)
                    min_eigenval = np.min(eigenvals)
                    if min_eigenval <= 0:
                        reg_param = abs(min_eigenval) + 1e-6
                        S_i += np.eye(S_i.shape[0]) * reg_param
                    self.within_class_cov[i] = S_i
            else:
                # Fallback for empty classes
                self.within_class_cov[i] = np.eye(T_train.shape[1]) * np.mean(self.pca.explained_variance_)

    self.is_fitted = True
    self.n_classes_ = k
    self.classes_ = self.encoder.categories_[0]

    return self

Distances and the decision rule:

analysis.py · PLSDA.predict() · lines 827–857

def predict(self, X):
    """
    Predict class labels for samples in X.
    """
    if not self.is_fitted:
        raise ValueError("Model must be fitted before prediction")

    distances = self.mahalanobis(X)

    if self.style == 'hard':
        # Assign to closest class
        class_indices = np.argmin(distances, axis=1)
        predictions = [self.classes_[idx] for idx in class_indices]
    else:
        # Soft classification with critical distance
        d_crit = stats.chi2.ppf(1 - self.alpha, self.n_classes_ - 1)

        predictions = []
        for row in distances:
            # Sort classes by distance and find those within critical distance
            class_distance_pairs = list(zip(self.classes_, row))
            class_distance_pairs.sort(key=lambda x: x[1])  # Sort by distance

            candidate_classes = [cls for cls, d in class_distance_pairs if d < d_crit]

            if candidate_classes:
                predictions.append(candidate_classes)
            else:
                predictions.append(["NOT_ASSIGNED"])

    return predictions

Full source of the PLSDA class

analysis.py · PLSDA() · lines 570–937

class PLSDA:
    """
    A numerically stable implementation of PLS-DA based on Pomerantsev and Rodionova's approach.
    Supports both hard and soft classification methods.
    """

    def __init__(self, n_components=5, style='hard', scale_x=True, alpha=0.05, gamma=0.05, random_state=None):
        """
        Initialize PLS-DA model.

        Parameters:
        -----------
        n_components : int, default=5
            Number of PLS components to use
        style : str, default='hard'
            Classification style: 'hard' or 'soft'
        scale_x : bool, default=True
            Whether to scale X features
        alpha : float, default=0.05
            Type I error rate for soft classification
        gamma : float, default=0.05
            Outlier detection rate for soft classification
        random_state : int, default=None
            Random state for reproducibility
        """
        self.n_components = n_components
        self.style = style
        self.scale_x = scale_x
        self.alpha = alpha
        self.gamma = gamma
        self.random_state = random_state

        # Initialize components
        self.pls = None
        self.pca = None
        self.x_scaler = None
        self.y_scaler = None
        self.encoder = None
        self.class_centers = None
        self.cov_matrix = None
        self.within_class_cov = None
        self.class_masks = None
        self.is_fitted = False

        # Store training data for proper transformation
        self._X_train_scaled = None
        self._y_train_encoded = None

    def _stable_scale(self, X, fit=False):
        """Numerically stable scaling with condition number check."""
        if fit:
            # Use the same scaler as pychemauth - center always, scale optionally
            self.x_scaler = StandardScaler(with_mean=True, with_std=self.scale_x)
            X_scaled = self.x_scaler.fit_transform(X)
        else:
            X_scaled = self.x_scaler.transform(X)

        # Check for potential numerical issues
        cond_num = np.linalg.cond(X_scaled)
        if cond_num > 1e10:
            warnings.warn(f"High condition number detected: {cond_num:.2e}. Consider regularization.")

        return X_scaled

    def _encode_y(self, y, fit=False):
        """One-hot encode y following pychemauth approach."""
        if fit:
            self.encoder = OneHotEncoder(sparse_output=False, handle_unknown="error")
            y_encoded = self.encoder.fit_transform(y.reshape(-1, 1))

            # Center Y but never scale (following pychemauth)
            self.y_scaler = StandardScaler(with_mean=True, with_std=False)
            y_encoded = self.y_scaler.fit_transform(y_encoded)
        else:
            y_encoded = self.encoder.transform(y.reshape(-1, 1))
            y_encoded = self.y_scaler.transform(y_encoded)

        return y_encoded

    def _safe_inverse(self, matrix, method='pinv', rcond=1e-12):
        """
        Safely compute matrix inverse with regularization if needed.
        """
        if method == 'pinv':
            return np.linalg.pinv(matrix, rcond=rcond, hermitian=True)
        else:
            cond_num = np.linalg.cond(matrix)
            if cond_num > 1e10:
                warnings.warn(f"High condition number {cond_num:.2e} in matrix inversion. Using pseudo-inverse.")
                return np.linalg.pinv(matrix, rcond=rcond, hermitian=True)
            try:
                return np.linalg.inv(matrix)
            except np.linalg.LinAlgError:
                warnings.warn("Matrix inversion failed. Using pseudo-inverse.")
                return np.linalg.pinv(matrix, rcond=rcond, hermitian=True)

    def fit(self, X, y):
        """
        Fit PLS-DA model to training data following pychemauth algorithm.
        """
        X = np.array(X, dtype=np.float64)
        y = np.array(y)

        # Step 1: Preprocess data (following pychemauth exactly)
        self._X_train_scaled = self._stable_scale(X, fit=True)
        self._y_train_encoded = self._encode_y(y, fit=True)

        # Step 2: Perform PLS2 regression
        k = len(self.encoder.categories_[0])  # number of classes

        # Bound n_components properly
        upper_bound = min(self._X_train_scaled.shape[0] - 1, self._X_train_scaled.shape[1])
        self.n_components = min(self.n_components, upper_bound)

        self.pls = PLSRegression(
            n_components=self.n_components,
            max_iter=10000,
            tol=1e-9,
            scale=False  # Already scaled
        )
        self.pls.fit(self._X_train_scaled, self._y_train_encoded)

        # Step 3: Get predictions and perform PCA (CRITICAL: follow pychemauth exactly)
        # Get PLS predictions and inverse transform Y scaling to get back to one-hot scale
        y_hat_train = self.y_scaler.inverse_transform(
            self.pls.predict(self._X_train_scaled)
        )

        # Perform PCA on y_hat_train with k-1 components
        n_pca_components = k - 1
        self.pca = PCA(n_components=n_pca_components, random_state=self.random_state)
        T_train = self.pca.fit_transform(y_hat_train)

        # Step 4: Compute class centers in PCA space (CRITICAL: pychemauth method)
        # Class centers are projections of identity matrix (one-hot vectors)
        identity_matrix = np.eye(k)
        self.class_centers = self.pca.transform(identity_matrix)

        # Step 5: Store class masks for training data
        self.class_masks = {}
        for i in range(k):
            class_label = self.encoder.categories_[0][i]
            self.class_masks[i] = (y == class_label)

        # Step 6: Compute covariance matrices based on classification style
        if self.style == 'hard':
            # For hard classification: L = diag(explained_variance)
            # This matches pychemauth: L = np.cov(T_train.T) which equals diag(explained_variance)
            self.cov_matrix = np.eye(len(self.pca.explained_variance_)) * self.pca.explained_variance_

            # Handle 1D case (binary classification)
            if self.cov_matrix.ndim == 0:
                self.cov_matrix = np.array([[self.cov_matrix]])

        else:  # soft classification
            # Compute within-class scatter matrices (following pychemauth exactly)
            self.within_class_cov = {}

            for i in range(k):
                mask = self.class_masks[i]
                if np.sum(mask) > 0:
                    # Get class samples in T space and center around class center
                    t_class = T_train[mask] - self.class_centers[i]

                    # Compute scatter matrix as in pychemauth
                    S_i = np.zeros((T_train.shape[1], T_train.shape[1]), dtype=np.float64)
                    for j in range(t_class.shape[0]):
                        t_vec = t_class[j, :].reshape(-1, 1)
                        S_i += np.dot(t_vec, t_vec.T)
                    S_i /= t_class.shape[0]  # Divide by number of samples

                    # Check if positive definite and regularize if needed
                    try:
                        np.linalg.cholesky(S_i)
                        self.within_class_cov[i] = S_i
                    except np.linalg.LinAlgError:
                        # Add small regularization to diagonal
                        eigenvals = np.linalg.eigvals(S_i)
                        min_eigenval = np.min(eigenvals)
                        if min_eigenval <= 0:
                            reg_param = abs(min_eigenval) + 1e-6
                            S_i += np.eye(S_i.shape[0]) * reg_param
                        self.within_class_cov[i] = S_i
                else:
                    # Fallback for empty classes
                    self.within_class_cov[i] = np.eye(T_train.shape[1]) * np.mean(self.pca.explained_variance_)

        self.is_fitted = True
        self.n_classes_ = k
        self.classes_ = self.encoder.categories_[0]

        return self

    def transform(self, X):
        """
        Transform data to PLS-PCA space following pychemauth algorithm.
        """
        if not self.is_fitted:
            raise ValueError("Model must be fitted before transformation")

        # Follow pychemauth transformation exactly
        X_scaled = self._stable_scale(X, fit=False)

        # Get PLS predictions and inverse transform Y scaling
        y_hat = self.y_scaler.inverse_transform(
            self.pls.predict(X_scaled)
        )

        # Apply PCA transformation
        T = self.pca.transform(y_hat)

        return T

    def mahalanobis(self, X, class_idx=None):
        """
        Compute Mahalanobis distances following pychemauth algorithm.
        """
        T = self.transform(X)

        if self.style == 'hard':
            # Use pooled covariance matrix L
            cov_inv = self._safe_inverse(self.cov_matrix)

            if class_idx is not None:
                distances = np.array([
                    np.dot(np.dot((t - self.class_centers[class_idx]), cov_inv),
                           (t - self.class_centers[class_idx]).T)
                    for t in T
                ])
            else:
                distances = np.array([
                    [np.dot(np.dot((t - center), cov_inv), (t - center).T)
                     for center in self.class_centers]
                    for t in T
                ])
        else:
            # Use within-class covariance matrices
            if class_idx is not None:
                cov_inv = self._safe_inverse(self.within_class_cov[class_idx])
                distances = np.array([
                    np.dot(np.dot((t - self.class_centers[class_idx]), cov_inv),
                           (t - self.class_centers[class_idx]).T)
                    for t in T
                ])
            else:
                distances = np.array([
                    [np.dot(np.dot((t - self.class_centers[i]),
                                   self._safe_inverse(self.within_class_cov[i])),
                            (t - self.class_centers[i]).T)
                     for i in range(self.n_classes_)]
                    for t in T
                ])

        # Ensure non-negative distances
        distances = np.maximum(distances, 0)
        return distances

    def predict(self, X):
        """
        Predict class labels for samples in X.
        """
        if not self.is_fitted:
            raise ValueError("Model must be fitted before prediction")

        distances = self.mahalanobis(X)

        if self.style == 'hard':
            # Assign to closest class
            class_indices = np.argmin(distances, axis=1)
            predictions = [self.classes_[idx] for idx in class_indices]
        else:
            # Soft classification with critical distance
            d_crit = stats.chi2.ppf(1 - self.alpha, self.n_classes_ - 1)

            predictions = []
            for row in distances:
                # Sort classes by distance and find those within critical distance
                class_distance_pairs = list(zip(self.classes_, row))
                class_distance_pairs.sort(key=lambda x: x[1])  # Sort by distance

                candidate_classes = [cls for cls, d in class_distance_pairs if d < d_crit]

                if candidate_classes:
                    predictions.append(candidate_classes)
                else:
                    predictions.append(["NOT_ASSIGNED"])

        return predictions

    def predict_proba(self, X):
        """
        Predict class probabilities based on Mahalanobis distances.
        """
        distances = self.mahalanobis(X)

        # Convert distances to probabilities using exponential
        p = np.exp(-np.clip(distances / 2.0, a_max=None, a_min=-500))

        if self.style == 'hard':
            # For hard classification, use softmax
            proba = (p.T / np.sum(p.T, axis=0)).T
        else:
            # For soft classification, normalize by determinant of covariance
            norm = np.zeros(self.n_classes_, dtype=np.float64)
            for i in range(self.n_classes_):
                norm[i] = np.sqrt(np.linalg.det(2.0 * np.pi * self.within_class_cov[i]))

            proba = np.array([
                [min(1.0, p_val) for p_val in row]
                for row in (p / norm)
            ], dtype=np.float64)

        return proba

    def get_loadings(self):
        """
        Get PLS loadings (weights) for feature importance visualization.
        """
        if not self.is_fitted:
            raise ValueError("Model must be fitted before getting loadings")

        return self.pls.x_weights_

    def get_correlation_loadings(self, X_original=None):
        """
        Get correlation loadings between original features and PLS components.
        """
        if not self.is_fitted:
            raise ValueError("Model must be fitted before getting correlation loadings")

        if X_original is None:
            return self.pls.x_weights_ / np.sqrt(np.sum(self.pls.x_weights_ ** 2, axis=0))

        X_scaled = self.x_scaler.transform(X_original)
        X_transformed = self.transform(X_original)

        n_features = X_scaled.shape[1]
        n_components = X_transformed.shape[1]
        correlation_loadings = np.zeros((n_features, n_components))

        for i in range(n_features):
            for j in range(n_components):
                correlation_loadings[i, j] = np.corrcoef(X_scaled[:, i], X_transformed[:, j])[0, 1]

        return correlation_loadings

    def get_feature_importance(self):
        """
        Get feature importance from PLS weights.
        """
        if not self.is_fitted:
            raise ValueError("Model must be fitted before getting feature importance")

        return np.sum(np.abs(self.pls.x_weights_), axis=1)

    def score(self, X, y):
        """
        Score model accuracy.
        """
        predictions = self.predict(X)
        y = np.array(y)

        if self.style == 'hard':
            correct = sum(1 for pred, true in zip(predictions, y) if pred == true)
        else:
            correct = sum(1 for pred, true in zip(predictions, y) if true in pred)

        return correct / len(y)

How run_plsda_analysis searches for the number of components and fits the final model:

analysis.py · run_plsda_analysis() · lines 3419–3518

# 4. Search for optimal n_components?
search_params = parameters.get('search_optimal_n_components', False)
    # ...
    for n_comp in range(min_n, max_n + 1):
        cv_res = _run_cv_evaluation(
            X_processed, y, n_comp, parameters, cv_method, cv_folds,
            groups_cv, overall_classes, parameters.get('random_state', 42)
        )
        score = cv_res[metric_map[metric]]
        if score is None:
            # If metric unavailable, skip this n (e.g., sensitivity for multiclass)
            continue
        search_scores.append((n_comp, score))
        if score > best_score:
            best_score = score
            best_n = n_comp
            best_cv_results = cv_res
# ...
plsda = PLSDA(
    n_components=n_components,
    style=style,
    alpha=alpha,
    gamma=gamma,
    scale_x=scale_x,
    random_state=random_state
)
plsda.fit(X_processed, y)

Full source of run_plsda_analysis()

analysis.py · run_plsda_analysis() · lines 3071–3721

def run_plsda_analysis(X, y, parameters, cv_method, cv_folds):
    """Run PLS-DA analysis with StablePLSDA class and optional feature selection"""
    from sklearn.model_selection import StratifiedKFold, GroupKFold, KFold
    from sklearn.metrics import classification_report, confusion_matrix, accuracy_score, roc_auc_score
    import numpy as np

    # ---------- Helper: run cross-validation for a given n_components ----------
    def _run_cv_evaluation(X_data, y_data, n_comp, params, cv_meth, cv_folds_num,
                           groups_data, overall_classes, random_state_val):
        """Perform cross-validation and return a dict with all CV metrics, including AUC and Q²."""
        # Extract relevant parameters
        style = params.get('style', 'hard')
        alpha = params.get('alpha', 0.05)
        gamma = params.get('gamma', 0.01)
        scale_x = params.get('scale_x', True)

        cv_scores = []
        cv_predictions = []
        cv_true_labels = []
        cv_class_reports = []
        cv_confusion_matrices = []
        cv_sensitivity_scores = []
        cv_specificity_scores = []

        # AUC initialisation
        cv_auc_scores = []
        cv_scores_all = []  # will store padded probabilities for all overall classes

        is_binary = len(overall_classes) == 2
        overall_classes_list = overall_classes.tolist()
        n_overall = len(overall_classes_list)

        # Set up CV splitter
        if cv_meth == 'group':
            if groups_data is None:
                raise ValueError("Groups parameter required for group cross-validation")
            cv = GroupKFold(n_splits=cv_folds_num)
            splitter_groups = groups_data
        elif cv_meth == 'stratified':
            cv = StratifiedKFold(n_splits=cv_folds_num, shuffle=True, random_state=random_state_val)
            splitter_groups = None
        else:
            cv = KFold(n_splits=cv_folds_num, shuffle=True, random_state=random_state_val)
            splitter_groups = None

        for train_idx, test_idx in cv.split(X_data, y_data, groups=splitter_groups):
            X_train, X_test = X_data[train_idx], X_data[test_idx]
            y_train, y_test = y_data[train_idx], y_data[test_idx]

            # Fit PLS-DA
            plsda_cv = PLSDA(
                n_components=n_comp,
                style=style,
                alpha=alpha,
                gamma=gamma,
                scale_x=scale_x,
                random_state=random_state_val
            )
            plsda_cv.fit(X_train, y_train)

            # Predict (classification)
            if style == 'hard':
                y_pred = plsda_cv.predict(X_test)
                y_pred_array = np.array(y_pred)
                accuracy = accuracy_score(y_test, y_pred_array)
                cv_predictions.extend(y_pred)
            else:
                y_pred_list = plsda_cv.predict(X_test)
                y_pred = []
                for pred in y_pred_list:
                    if pred and pred[0] != "NOT_ASSIGNED":
                        y_pred.append(pred[0])
                    else:
                        y_pred.append("UNKNOWN")
                y_pred_array = np.array(y_pred)
                accuracy = accuracy_score(y_test, y_pred_array)
                cv_predictions.extend(y_pred)

            cv_scores.append(accuracy)
            cv_true_labels.extend(y_test)

            # ---------- AUC computation (REVISED: always reorder to overall_classes) ----------
            y_score = None
            model_classes = plsda_cv.classes_

            # Try to get probabilities first (preferred)
            if hasattr(plsda_cv, "predict_proba"):
                y_prob = plsda_cv.predict_proba(X_test)  # shape (n_samples, len(model_classes))
                n_samples = X_test.shape[0]
                # Create padded matrix in overall_classes order
                y_score_padded = np.zeros((n_samples, n_overall))
                class_to_idx = {cls: i for i, cls in enumerate(overall_classes_list)}
                for i, cls in enumerate(model_classes):
                    if cls in class_to_idx:
                        col = class_to_idx[cls]
                        y_score_padded[:, col] = y_prob[:, i]
                # For binary, take the positive class (second column)
                if is_binary:
                    y_score = y_score_padded[:, 1] if y_score_padded.shape[1] == 2 else None
                else:
                    # Multi-class: keep only columns for classes present in y_test (for roc_auc_score)
                    classes_in_test = np.unique(y_test)
                    test_indices = [class_to_idx[cls] for cls in classes_in_test]
                    y_score = y_score_padded[:, test_indices] if test_indices else None
                # Store padded probabilities for aggregated AUC
                cv_scores_all.extend(y_score_padded.tolist())

            elif hasattr(plsda_cv, "decision_function"):
                # Fallback: use decision function (negative distances as scores)
                y_score_raw = -plsda_cv.decision_function(X_test)
                n_samples = X_test.shape[0]
                y_score_padded = np.full((n_samples, n_overall), -np.inf)
                class_to_idx = {cls: i for i, cls in enumerate(overall_classes_list)}
                if is_binary:
                    # For binary, decision_function may return 1D (positive class) or 2D
                    if y_score_raw.ndim == 1:
                        # Single column: treat as positive class score (if model_classes[0] is negative)
                        # We assume binary classification: overall_classes[0] = negative, overall_classes[1] = positive
                        pos_col = class_to_idx.get(overall_classes_list[1], 1)
                        neg_col = class_to_idx.get(overall_classes_list[0], 0)
                        y_score_padded[:, pos_col] = y_score_raw
                        y_score_padded[:, neg_col] = -y_score_raw  # approximate negative score
                    else:
                        # 2D: map columns
                        for i, cls in enumerate(model_classes):
                            if cls in class_to_idx:
                                col = class_to_idx[cls]
                                y_score_padded[:, col] = y_score_raw[:, i]
                    # For binary, positive score is second column
                    if y_score_padded.shape[1] == 2:
                        y_score = y_score_padded[:, 1]
                    else:
                        y_score = None
                else:
                    # Multi-class: map columns
                    for i, cls in enumerate(model_classes):
                        if cls in class_to_idx:
                            col = class_to_idx[cls]
                            y_score_padded[:, col] = y_score_raw[:, i] if y_score_raw.ndim == 2 else y_score_raw
                    # Keep columns for classes present in y_test
                    classes_in_test = np.unique(y_test)
                    test_indices = [class_to_idx[cls] for cls in classes_in_test]
                    y_score = y_score_padded[:, test_indices] if test_indices else None
                cv_scores_all.extend(y_score_padded.tolist())
            else:
                y_score = None

            # Compute AUC for this fold
            auc_fold = None
            if y_score is not None and len(np.unique(y_test)) >= 2:
                try:
                    if is_binary:
                        auc_fold = roc_auc_score(y_test, y_score)
                    else:
                        auc_fold = roc_auc_score(y_test, y_score, multi_class='ovr', average='macro')
                except ValueError:
                    auc_fold = None
            cv_auc_scores.append(auc_fold)

            # ---- Classification report and confusion matrix ----
            fold_report = classification_report(y_test, y_pred_array, output_dict=True, zero_division=0)
            cv_class_reports.append(fold_report)

            cm = confusion_matrix(y_test, y_pred_array, labels=overall_classes_list)
            cv_confusion_matrices.append(cm)

            if is_binary:
                tn, fp, fn, tp = cm.ravel()
                sens = tp / (tp + fn) if (tp + fn) > 0 else 0
                spec = tn / (tn + fp) if (tn + fp) > 0 else 0
                cv_sensitivity_scores.append(sens)
                cv_specificity_scores.append(spec)
            else:
                cv_sensitivity_scores.append(None)
                cv_specificity_scores.append(None)

        # --- Aggregate AUC ---
        auc_aggregated = None
        if cv_scores_all:
            scores_matrix = np.array(cv_scores_all)
            # For binary, scores_matrix has 2 columns (negative, positive) in overall_classes order
            if scores_matrix.shape[1] == n_overall or (is_binary and scores_matrix.shape[1] == 2):
                try:
                    if is_binary:
                        pos_scores = scores_matrix[:, 1]
                        auc_aggregated = roc_auc_score(cv_true_labels, pos_scores)
                    else:
                        auc_aggregated = roc_auc_score(cv_true_labels, scores_matrix,
                                                       multi_class='ovr', average='macro')
                except ValueError:
                    auc_aggregated = None

        # Filter valid AUC scores for mean/std
        valid_auc = [a for a in cv_auc_scores if a is not None]
        cv_auc_mean = np.mean(valid_auc) if valid_auc else None
        cv_auc_std = np.std(valid_auc) if valid_auc else None

        # --- Compute Q² (cross-validated R²) for binary ---
        cv_r2 = None
        if is_binary:
            class_to_num = {overall_classes_list[0]: 0, overall_classes_list[1]: 1}
            y_true_num = np.array([class_to_num.get(label, np.nan) for label in cv_true_labels])
            y_pred_num = np.array([class_to_num.get(pred, 0.5) for pred in cv_predictions])
            mask = ~np.isnan(y_true_num) & ~np.isnan(y_pred_num)
            y_true_num = y_true_num[mask]
            y_pred_num = y_pred_num[mask]
            if len(y_true_num) > 0:
                mean_y = np.mean(y_true_num)
                ss_tot = np.sum((y_true_num - mean_y) ** 2)
                if ss_tot > 0:
                    ss_res = np.sum((y_true_num - y_pred_num) ** 2)
                    cv_r2 = 1 - ss_res / ss_tot

        # --- Other aggregated metrics (precision, recall, F1, etc.) ---
        cv_precision_scores = []
        cv_recall_scores = []
        cv_f1_scores = []
        for report in cv_class_reports:
            if 'macro avg' in report:
                cv_precision_scores.append(report['macro avg']['precision'])
                cv_recall_scores.append(report['macro avg']['recall'])
                cv_f1_scores.append(report['macro avg']['f1-score'])

        # Aggregated confusion matrix
        cv_confusion_matrix_aggregated = confusion_matrix(
            cv_true_labels, cv_predictions, labels=overall_classes_list
        ).tolist()

        # Aggregated sensitivity/specificity (binary)
        if is_binary:
            cm_agg = np.array(cv_confusion_matrix_aggregated)
            tn, fp, fn, tp = cm_agg.ravel()
            cv_sensitivity_aggregated = tp / (tp + fn) if (tp + fn) > 0 else 0
            cv_specificity_aggregated = tn / (tn + fp) if (tn + fp) > 0 else 0
        else:
            cv_sensitivity_aggregated = None
            cv_specificity_aggregated = None

        # Average confusion matrix (via external helper)
        if cv_confusion_matrices:
            cv_confusion_matrix_avg, avg_matrix_classes = create_average_confusion_matrix(
                cv_confusion_matrices, [overall_classes_list] * len(cv_confusion_matrices)
            )
            if is_binary:
                cm_avg = np.array(cv_confusion_matrix_avg)
                tn, fp, fn, tp = cm_avg.ravel()
                cv_sensitivity_avg = tp / (tp + fn) if (tp + fn) > 0 else 0
                cv_specificity_avg = tn / (tn + fp) if (tn + fp) > 0 else 0
            else:
                cv_sensitivity_avg = None
                cv_specificity_avg = None
        else:
            cv_confusion_matrix_avg = []
            avg_matrix_classes = []
            cv_sensitivity_avg = None
            cv_specificity_avg = None

        # Mean and std of per‑fold metrics
        cv_scores_filtered = [s for s in cv_scores if s is not None]
        cv_mean = np.mean(cv_scores_filtered) if cv_scores_filtered else 0.0
        cv_std = np.std(cv_scores_filtered) if cv_scores_filtered else 0.0

        cv_precision_mean = np.mean(cv_precision_scores) if cv_precision_scores else 0.0
        cv_precision_std = np.std(cv_precision_scores) if cv_precision_scores else 0.0
        cv_recall_mean = np.mean(cv_recall_scores) if cv_recall_scores else 0.0
        cv_recall_std = np.std(cv_recall_scores) if cv_recall_scores else 0.0
        cv_f1_mean = np.mean(cv_f1_scores) if cv_f1_scores else 0.0
        cv_f1_std = np.std(cv_f1_scores) if cv_f1_scores else 0.0

        # Sensitivity / specificity per‑fold mean and std (filter None)
        sens_vals = [s for s in cv_sensitivity_scores if s is not None]
        spec_vals = [s for s in cv_specificity_scores if s is not None]
        cv_sensitivity_mean = np.mean(sens_vals) if sens_vals else None
        cv_sensitivity_std = np.std(sens_vals) if sens_vals else None
        cv_specificity_mean = np.mean(spec_vals) if spec_vals else None
        cv_specificity_std = np.std(spec_vals) if spec_vals else None

        # Balanced accuracy (for binary only)
        if is_binary and cv_sensitivity_mean is not None and cv_specificity_mean is not None:
            cv_balanced_accuracy = (cv_sensitivity_mean + cv_specificity_mean) / 2
        else:
            cv_balanced_accuracy = None

        # Return all metrics including AUC and Q²
        return {
            'cv_scores': cv_scores,
            'cv_predictions': cv_predictions,
            'cv_true_labels': cv_true_labels,
            'cv_class_reports': cv_class_reports,
            'cv_confusion_matrices': cv_confusion_matrices,
            'cv_sensitivity_scores': cv_sensitivity_scores,
            'cv_specificity_scores': cv_specificity_scores,
            'cv_confusion_matrix_aggregated': cv_confusion_matrix_aggregated,
            'cv_confusion_matrix_avg': cv_confusion_matrix_avg,
            'avg_matrix_classes': avg_matrix_classes,
            'cv_mean': cv_mean,
            'cv_std': cv_std,
            'cv_precision': cv_precision_mean,
            'cv_precision_std': cv_precision_std,
            'cv_recall': cv_recall_mean,
            'cv_recall_std': cv_recall_std,
            'cv_f1': cv_f1_mean,
            'cv_f1_std': cv_f1_std,
            'cv_sensitivity': cv_sensitivity_mean,
            'cv_sensitivity_std': cv_sensitivity_std,
            'cv_sensitivity_aggregated': cv_sensitivity_aggregated,
            'cv_sensitivity_avg': cv_sensitivity_avg,
            'cv_specificity': cv_specificity_mean,
            'cv_specificity_std': cv_specificity_std,
            'cv_specificity_aggregated': cv_specificity_aggregated,
            'cv_specificity_avg': cv_specificity_avg,
            'cv_balanced_accuracy': cv_balanced_accuracy,
            'is_binary': is_binary,
            # AUC fields
            'cv_auc': sanitize_for_json(cv_auc_mean),
            'cv_auc_std': sanitize_for_json(cv_auc_std),
            'cv_auc_scores': sanitize_for_json(cv_auc_scores),
            'cv_auc_aggregated': sanitize_for_json(auc_aggregated),
            # Q² (cross-validated R²)
            'cv_r2': cv_r2,
        }

    # ---------- Main function body ----------
    # 1. Handle average_replicates
    if parameters.get('average_replicates', False):
        groups = parameters.get('groups')
        if groups is None:
            raise ValueError("'groups' must be provided when averaging replicates is enabled.")
        X, y, new_groups = average_replicates(X, y, groups, method='mean')
        parameters['groups'] = new_groups

    # 2. Get overall classes and other base parameters
    overall_classes = np.unique(y)
    overall_classes_list = overall_classes.tolist()
    n_overall_classes = len(overall_classes_list)
    is_binary = n_overall_classes == 2

    # 3. Feature selection (if any)
    feature_selection_params = get_feature_selection_params(parameters)
    if feature_selection_params:
        X_processed, selected_features, feature_scores = select_features(
            X, y, feature_selection_params, 'classification'
        )
    else:
        X_processed = X
        selected_features = list(range(X.shape[1]))
        feature_scores = None

    # 4. Search for optimal n_components?
    search_params = parameters.get('search_optimal_n_components', False)
    if search_params is not False:
        # Validate search_params
        if not isinstance(search_params, dict):
            raise ValueError("search_optimal_n_components must be a dict or False")
        required_keys = {'metric', 'min_n', 'max_n'}
        if not required_keys.issubset(search_params.keys()):
            raise ValueError(f"search_optimal_n_components must contain keys: {required_keys}")
        metric = search_params['metric']
        min_n = int(search_params['min_n'])
        max_n = int(search_params['max_n'])
        if min_n > max_n:
            raise ValueError("min_n must be <= max_n")

        # Map metric name to a key in cv_results
        metric_map = {
            'accuracy': 'cv_mean',
            'f1_macro': 'cv_f1',
            'f1_weighted': 'cv_f1_weighted',  # not computed yet; we can add if needed
            'precision_macro': 'cv_precision',
            'recall_macro': 'cv_recall',
            'sensitivity': 'cv_sensitivity',
            'specificity': 'cv_specificity',
            'balanced_accuracy': 'cv_balanced_accuracy',
            'r2': 'cv_r2'
        }
        if metric not in metric_map:
            raise ValueError(f"Unsupported metric: {metric}. Supported: {list(metric_map.keys())}")

        # Ensure metric is appropriate for binary/multiclass
        if metric in ('sensitivity', 'specificity', 'balanced_accuracy', 'r2') and not is_binary:
            raise ValueError(f"Metric '{metric}' is only defined for binary classification.")

        # Prepare groups for CV (if needed)
        groups_cv = parameters.get('groups') if cv_method == 'group' else None

        # Loop over n_components
        best_score = -np.inf
        best_n = None
        best_cv_results = None
        search_scores = []

        # Cap max_n to number of features (if X_processed has fewer features)
        max_n = min(max_n, X_processed.shape[1])

        for n_comp in range(min_n, max_n + 1):
            cv_res = _run_cv_evaluation(
                X_processed, y, n_comp, parameters, cv_method, cv_folds,
                groups_cv, overall_classes, parameters.get('random_state', 42)
            )
            score = cv_res[metric_map[metric]]
            if score is None:
                # If metric unavailable, skip this n (e.g., sensitivity for multiclass)
                continue
            search_scores.append((n_comp, score))
            if score > best_score:
                best_score = score
                best_n = n_comp
                best_cv_results = cv_res

        if best_n is None:
            raise ValueError("No valid n_components found in the search range.")

        # Use best n and its CV results
        n_components = best_n
        cv_results = best_cv_results
        # Store search info for output
        search_info = {
            'best_n_components': best_n,
            'metric': metric,
            'scores': search_scores
        }
    else:
        # No search: use given n_components (default 7)
        n_components = parameters.get('n_components', 7)
        groups_cv = parameters.get('groups') if cv_method == 'group' else None
        cv_results = _run_cv_evaluation(
            X_processed, y, n_components, parameters, cv_method, cv_folds,
            groups_cv, overall_classes, parameters.get('random_state', 42)
        )
        search_info = None

    # 5. Fit final model on all data with the selected n_components
    style = parameters.get('style', 'hard')
    alpha = parameters.get('alpha', 0.05)
    gamma = parameters.get('gamma', 0.01)
    scale_x = parameters.get('scale_x', True)
    random_state = parameters.get('random_state', 42)
    loading_type = parameters.get('loading_type', 'weights')

    plsda = PLSDA(
        n_components=n_components,
        style=style,
        alpha=alpha,
        gamma=gamma,
        scale_x=scale_x,
        random_state=random_state
    )
    plsda.fit(X_processed, y)

    # 6. Compute loadings, predictions, probabilities, etc.
    if loading_type == 'correlation':
        loadings = plsda.get_correlation_loadings(X_processed)
    else:
        loadings = plsda.get_loadings()

    if style == 'hard':
        y_pred_full = plsda.predict(X_processed)
        y_pred_array_full = np.array(y_pred_full)
        y_pred_proba = plsda.predict_proba(X_processed)
    else:
        y_pred_list = plsda.predict(X_processed)
        y_pred_full = []
        for pred in y_pred_list:
            if pred and pred[0] != "NOT_ASSIGNED":
                y_pred_full.append(pred[0])
            else:
                y_pred_full.append("UNKNOWN")
        y_pred_array_full = np.array(y_pred_full)
        y_pred_proba = plsda.predict_proba(X_processed)

    final_classes = plsda.classes_.tolist()

    # Full model confusion matrix and sensitivity/specificity (binary)
    full_confusion_matrix = confusion_matrix(y, y_pred_array_full, labels=final_classes)
    if is_binary and full_confusion_matrix.shape == (2, 2):
        tn, fp, fn, tp = full_confusion_matrix.ravel()
        full_sensitivity = tp / (tp + fn) if (tp + fn) > 0 else 0
        full_specificity = tn / (tn + fp) if (tn + fp) > 0 else 0
    else:
        full_sensitivity = None
        full_specificity = None

    # 7. Compute ordinary R² on the full training set (binary only)
    r2 = None
    if is_binary:
        class_to_num = {overall_classes_list[0]: 0, overall_classes_list[1]: 1}
        y_true_num = np.array([class_to_num.get(label, np.nan) for label in y])
        y_pred_num = np.array([class_to_num.get(pred, 0.5) for pred in y_pred_full])
        mask = ~np.isnan(y_true_num) & ~np.isnan(y_pred_num)
        y_true_num = y_true_num[mask]
        y_pred_num = y_pred_num[mask]
        if len(y_true_num) > 0:
            mean_y = np.mean(y_true_num)
            ss_tot = np.sum((y_true_num - mean_y) ** 2)
            if ss_tot > 0:
                ss_res = np.sum((y_true_num - y_pred_num) ** 2)
                r2 = 1 - ss_res / ss_tot

    # 8. Build results dict using cv_results from the (possibly searched) best model
    # Reorder average confusion matrix to match final_classes
    cv_confusion_matrix_avg = cv_results['cv_confusion_matrix_avg']
    avg_matrix_classes = cv_results['avg_matrix_classes']
    if cv_confusion_matrix_avg and avg_matrix_classes != final_classes:
        avg_to_final = {cls: idx for idx, cls in enumerate(avg_matrix_classes)}
        n_classes = len(final_classes)
        reordered_matrix = np.zeros((n_classes, n_classes))
        for i, true_cls in enumerate(final_classes):
            for j, pred_cls in enumerate(final_classes):
                if true_cls in avg_to_final and pred_cls in avg_to_final:
                    orig_i = avg_to_final[true_cls]
                    orig_j = avg_to_final[pred_cls]
                    if (orig_i < len(cv_confusion_matrix_avg) and
                            orig_j < len(cv_confusion_matrix_avg[0])):
                        reordered_matrix[i, j] = cv_confusion_matrix_avg[orig_i][orig_j]
        cv_confusion_matrix_avg = reordered_matrix.tolist()
        avg_matrix_classes = final_classes

    # Validation message for binary
    validation_message = None
    if is_binary and full_sensitivity is not None and full_specificity is not None:
        avg_sens_spec = (full_sensitivity + full_specificity) / 2
        accuracy = accuracy_score(y, y_pred_array_full)
        if abs(accuracy - avg_sens_spec) > 0.01:
            validation_message = f"Note: Accuracy ({accuracy:.3f}) differs from (sensitivity + specificity)/2 ({avg_sens_spec:.3f})"

    accuracy = accuracy_score(y, y_pred_array_full)

    # After y_pred_proba is defined (around line after fit)
    auc = None
    if is_binary:
        try:
            # For binary, positive class probability is second column (index 1) if classes are sorted
            # We assume overall_classes[1] is positive
            pos_proba = y_pred_proba[:, 1] if y_pred_proba.shape[1] == 2 else None
            if pos_proba is not None:
                auc = roc_auc_score(y, pos_proba)
        except Exception:
            auc = None

    # Explained variance
    try:
        X_transformed = plsda.transform(X_processed)
        total_variance = np.var(X_transformed, axis=0, ddof=1)
        explained_variance_ratio = (total_variance / np.sum(total_variance)).tolist()
        cumulative_variance = np.cumsum(explained_variance_ratio).tolist()
    except Exception as e:
        print(f"Warning: Could not calculate explained variance ratio: {e}")
        n_comp = min(n_components, X_processed.shape[1])
        explained_variance_ratio = [1.0 / n_comp] * n_comp
        cumulative_variance = np.cumsum(explained_variance_ratio).tolist()

    # Classification report
    try:
        class_report = classification_report(y, y_pred_array_full, output_dict=True, zero_division=0)
    except ValueError:
        unique_train = np.unique(y)
        unique_pred = np.unique(y_pred_array_full)
        all_classes = np.union1d(unique_train, unique_pred)
        class_report = classification_report(y, y_pred_array_full, labels=all_classes, output_dict=True,
                                             zero_division=0)

    # Prepare final results
    results = {
        'method': f'PLS-DA ({style})',
        'cv_scores': cv_results['cv_scores'],
        'cv_mean': cv_results['cv_mean'],
        'cv_std': cv_results['cv_std'],
        'cv_precision': cv_results['cv_precision'],
        'cv_precision_std': cv_results['cv_precision_std'],
        'cv_recall': cv_results['cv_recall'],
        'cv_recall_std': cv_results['cv_recall_std'],
        'cv_f1': cv_results['cv_f1'],
        'cv_f1_std': cv_results['cv_f1_std'],
        'cv_sensitivity': cv_results['cv_sensitivity'],
        'cv_sensitivity_std': cv_results['cv_sensitivity_std'],
        'cv_sensitivity_aggregated': cv_results['cv_sensitivity_aggregated'],
        'cv_sensitivity_avg': cv_results['cv_sensitivity_avg'],
        'cv_sensitivity_scores': cv_results['cv_sensitivity_scores'],
        'cv_specificity': cv_results['cv_specificity'],
        'cv_specificity_std': cv_results['cv_specificity_std'],
        'cv_specificity_aggregated': cv_results['cv_specificity_aggregated'],
        'cv_specificity_avg': cv_results['cv_specificity_avg'],
        'cv_specificity_scores': cv_results['cv_specificity_scores'],
        # AUC fields
        'cv_auc': cv_results['cv_auc'],
        'cv_auc_std': cv_results['cv_auc_std'],
        'cv_auc_scores': cv_results['cv_auc_scores'],
        'cv_auc_aggregated': cv_results['cv_auc_aggregated'],
        # R² and Q²
        'r2': r2,
        'cv_r2': cv_results['cv_r2'],
        'is_binary': is_binary,
        'validation_message': validation_message,
        'sensitivity': full_sensitivity,
        'specificity': full_specificity,
        'auc': auc,
        'accuracy': accuracy,
        'cv_accuracy': cv_results['cv_mean'],
        'cv_confusion_matrix': cv_confusion_matrix_avg,
        'cv_confusion_matrices': [cm.tolist() for cm in cv_results['cv_confusion_matrices']],
        'cv_confusion_matrix_aggregated': cv_results['cv_confusion_matrix_aggregated'],
        'cv_class_reports': cv_results['cv_class_reports'],
        'classification_report': class_report,
        'confusion_matrix': full_confusion_matrix.tolist(),
        'transformed_data': sanitize_for_json(plsda.transform(X_processed).tolist()),
        'labels': y.tolist(),
        'loadings': sanitize_for_json(loadings.tolist()),
        'loading_type': loading_type,
        'predictions': y_pred_full,
        'prediction_probabilities': sanitize_for_json(y_pred_proba.tolist()),
        'classes': final_classes,
        'model_info': {
            'n_components': n_components,
            'style': style,
            'alpha': alpha,
            'gamma': gamma,
            'scale_x': scale_x,
            'random_state': random_state
        },
        'feature_importance': sanitize_for_json(plsda.get_feature_importance().tolist()),
        'explained_variance_ratio': explained_variance_ratio,
        'cumulative_variance': cumulative_variance
    }

    # Add binary consistency check
    if is_binary and full_sensitivity is not None and full_specificity is not None:
        results['metric_consistency_check'] = {
            'accuracy': accuracy_score(y, y_pred_array_full),
            'sensitivity': full_sensitivity,
            'specificity': full_specificity,
            'avg_sens_spec': (full_sensitivity + full_specificity) / 2,
            'difference': abs(accuracy_score(y, y_pred_array_full) - (full_sensitivity + full_specificity) / 2)
        }

    # Add feature selection info
    if feature_selection_params:
        results['feature_selection'] = {
            'selected_features': selected_features,
            'feature_scores': feature_scores,
            'method': feature_selection_params.get('method', 'anova'),
            'k': feature_selection_params.get('k', 'all')
        }

    # Add search info if used
    if search_info is not None:
        results['optimal_n_components_search'] = search_info

    results = sanitize_for_json(results)

    # Final fitted model for inference on unknown samples (popped in run_analysis)
    return attach_model_bundle(results, plsda, selected_features, feature_selection_params)

Limitations and common pitfalls

• Too many components fit noise: training accuracy grows while CV accuracy stays flat or drops. Keep the search range reasonable.
• Hard rule and new classes. The hard rule always picks one of the known classes, even for a sample that looks like none of them. Use the soft rule if foreign samples are possible.
• Soft rule and small classes. The per-class covariance needs enough samples in every class; with very few, the acceptance areas are unstable.
• gamma has no effect in the current version; outliers are not detected separately.

References

• Barker M., Rayens W. Partial least squares for discrimination. Journal of Chemometrics, 17, 166–173 (2003). doi:10.1002/cem.785
• Pomerantsev A. L., Rodionova O. Ye. Multiclass partial least squares discriminant analysis: taking the right way — a critical tutorial. Journal of Chemometrics, 32, e3030 (2018). doi:10.1002/cem.3030
• Wold S., Sjöström M., Eriksson L. PLS-regression: a basic tool of chemometrics. Chemometrics and Intelligent Laboratory Systems, 58, 109–130 (2001). doi:10.1016/S0169-7439(01)00155-1
• PyChemAuth: github.com/mahynski/pychemauth.