Partial Least Squares Discriminant Analysis

Partial Least Squares Discriminant Analysis is a supervised classification approach based on Partial Least Squares. Westerhuis et al. describe PLS-DA as a classification tool based on the PLS model in which the dependent variable is chosen to represent class membership. Szymańska et al. write the same idea as PLS regression with a special dummy response.

The response is therefore not a continuous chemical quantity such as concentration. Class labels are encoded as numbers, PLS relates $X$ to that coded response, and the continuous PLS prediction is then converted into a class decision by an explicit rule. The coding step and the decision step are part of the method. They are not optional commentary.

Standard PLS-DA is not PCA followed by LDA. Barker and Rayens discuss a statistical relationship between PLS, canonical correlation, and linear discriminant analysis. That relationship is not the computational procedure used on this page. The procedure here is class encoding, PLS regression, and a stated assignment rule.

PCA does not use class labels. PLS-DA does. Class membership enters the fit through the encoded response, so the latent directions are supervised. Barker and Rayens argue that PLS is to be preferred over PCA when discrimination is the goal and dimension reduction is needed. That statement is not a claim that PLS-DA is superior to LDA, PCR, or other classifiers.

The same supervision that can emphasize group structure also creates the central risk. Westerhuis et al. show that PLS-DA can separate randomly assigned labels in high-dimensional small- sample data. Apparent class structure in a fitted score plot is therefore not independent evidence of discrimination. Predictive validation is required.

The PLS latent-variable extraction is the same family of two-block methods documented on the PLSR page. What changes is the response and the layer that turns a continuous prediction into a class.

1. Encode the class labels as a numeric response. This page uses $0/1$ coding for two classes.
2. Learn preprocessing from training observations only, then fit PLS of $X$ on the encoded response with $A$ components.
3. Predict a continuous dummy response for new observations after applying the training preprocessing.
4. Convert that continuous value into a class with the stated decision rule.

The continuous prediction is not a calibrated class probability unless a separate probability model is constructed. Szymańska et al. emphasize that the PLS prediction can take any real value, not only the coded class numbers, and that translation into membership is a distinct classification step.

Let $X\in {\mathbb{R}}^{n\times p}$ have observations as rows. For two classes this page encodes membership as $y\in \{0,1{\}}^n$, following the $0/1$ convention used by Westerhuis et al. in their experimental section. Westerhuis et al. also allow the other class to be coded as $-1$ when the positive class is $1$. Szymańska et al. use $-1/+1$ coding. Those conventions are not mixed in one derivation.

For $C>2$ classes, Westerhuis et al. state that dummy variables are defined and a PLS2 algorithm is used. The dummy matrix used here as the direct extension of Equation 1 is

Algorithm 1 is the classification workflow. The PLS solver may be NIPALS, SIMPLS, or another documented PLS implementation. PLS-DA is not itself a new latent decomposition.

Binary PLS-DA on this page is PLS1 on a single coded vector. Multiclass PLS-DA, as named by Westerhuis et al., uses dummy variables and PLS2. That is a joint model of the dummy block. It is not automatically a set of independent one-versus-rest PLS models.

If an explicit dummy matrix is built in software, construct it as a response matrix. scikit-learn documents LabelBinarizer for target indicator coding. OneHotEncoder is documented as a feature encoder. The educational NumPy dummy matrix matches Equation 5 without treating labels as ordinary feature columns.

The number of PLS components is a hyperparameter. Do not choose it from training accuracy, from a fitted score plot, or from a full-data $Q^2$. Westerhuis et al. require that the complete modelling procedure be cross-validated. The predicted observation must not participate in model development, including pretreatment, scaling, and component selection.

Their FIT, 1CV, and 2CV comparison is the same nested structure documented on the Cross-Validation page. FIT resubstitutes after selecting $A$. Single CV still uses the predicted samples to choose $A$. Cross-model validation (2CV) holds the test observation out of both fitting and component selection. Only 2CV recovered the expected permutation null for the number of misclassifications in their study.

Score plots of fitted PLS-DA scores are supervised representations. Westerhuis et al. obtained clear score separation after randomly assigning class labels, including a case with $Q^2=-0.18$. They advocate against using PLS-DA score plots to infer between-class differences. Within- class structure may still be visible because the model is not forced to show it. Classification evidence must come from predictions.

Permutation testing, as used by Westerhuis et al. and Szymańska et al., builds a reference distribution of a diagnostic under random class labeling. It can ask whether an observed number of misclassifications, AUROC, or $Q^2$ exceeds what occurs by chance on that dataset. It does not prove biological validity, and it does not replace held-out prediction. This page does not prescribe a universal permutation count.

Westerhuis et al. define

$Q^2=1-{\sum }_i(y_i-{\hat{y}}_i{)}^2/{\sum }_i(y_i-\bar{y}{)}^2$. They warn that the numerical meaning of a given $Q^2$ for class discrimination is not obvious without a suitable null distribution. Szymańska et al. likewise advise against reading a single $Q^2$ value as proof of discrimination, and they recommend number of misclassifications and AUROC as more direct two-class diagnostics. No threshold such as $Q^2>0.5$ is used here.

UCI Wine, cultivars 0 and 1 only (59 and 71 samples, 13 chemical variables). Cultivar 2 is excluded. The figure is a PLS-DA score plot of $\mathrm{LV}1$ versus $\mathrm{LV}2$ from a two-component fit after column standardization with StandardScaler (population standard deviation). That is not the sample-standard- deviation autoscaling defined on the Normalization and Scaling page. The plotted coordinates are latent-variable scores $T$, not the dummy-response predictions $\hat{y}$. These are supervised scores. Apparent class structure in this plot is not independent validation.

Classification numbers in the caption are nested out-of-fold assignments under the midpoint rule, not a reading of the score plot. Inner 4-fold CV selects $A$ by lowest number of misclassifications. Outer 5-fold CV produces the held-out dummy-response predictions. This two-cultivar subset is not a difficult discrimination.

scikit-learn does not provide a PLSDA estimator. The Python listing fits PLSRegression on a $0/1$ response inside a Pipeline with StandardScaler and scale=False. MATLAB plsregress is SIMPLS regression, not a dedicated classifier. The same coding, column scaling of $X$, and midpoint rule are applied after $\hat{y}=[1X_{\mathrm{scaled}}]\mathrm{BETA}$. Compare predicted dummy responses and class assignments, not latent weights. sklearn NIPALS and MATLAB SIMPLS need not share intermediate matrices.

Sensitivity and specificity follow Altman and Bland as used by Szymańska et al.: $\mathrm{Se}=\mathrm{TP}/(\mathrm{TP}+\mathrm{FN})$ and $\mathrm{Sp}=\mathrm{TN}/(\mathrm{TN}+\mathrm{FP})$. Raw accuracy can be misleading under class imbalance. Broader classification metrics are on the Classification Metrics page. Undefined sensitivity or specificity is reported as NaN there. The educational listing here raises if a class is absent.

import numpy as npfrom sklearn.cross_decomposition import PLSRegressionfrom sklearn.model_selection import StratifiedKFoldfrom sklearn.pipeline import Pipelinefrom sklearn.preprocessing import StandardScaler  def encode_binary_01(labels, class0, class1):    labels = np.asarray(labels)    mask = (labels == class0) | (labels == class1)    if not np.any(mask):        raise ValueError("No observations belong to the requested classes.")    y = np.where(labels[mask] == class1, 1.0, 0.0)    return y, mask  def assign_binary_midpoint(y_hat, threshold=0.5):    y_hat = np.asarray(y_hat, dtype=float).reshape(-1)    if not np.all(np.isfinite(y_hat)):        raise ValueError("y_hat must contain only finite values.")    return (y_hat >= threshold).astype(int)  def dummy_class_matrix(labels):    labels = np.asarray(labels)    classes, inverse = np.unique(labels, return_inverse=True)    Y = np.zeros((labels.shape[0], classes.shape[0]), dtype=float)    Y[np.arange(labels.shape[0]), inverse] = 1.0    return Y, classes  def make_plsda_pipeline(n_components):    n_components = int(n_components)    if n_components < 1:        raise ValueError("n_components must be a positive integer.")    return Pipeline(        [            ("scale", StandardScaler()),            ("pls", PLSRegression(n_components=n_components, scale=False)),        ]    )  def _validate_xy(X, y):    X = np.asarray(X, dtype=float)    y = np.asarray(y, dtype=float).reshape(-1)    if X.ndim != 2:        raise ValueError("X must be a 2D array of samples by variables.")    if X.shape[0] != y.shape[0]:        raise ValueError("X and y must have the same number of observations.")    if min(X.shape) == 0:        raise ValueError("X must have at least one row and one column.")    if not np.all(np.isfinite(X)) or not np.all(np.isfinite(y)):        raise ValueError("X and y must contain only finite values.")    return X, y  def fit_plsda_binary(X, y, n_components):    X, y = _validate_xy(X, y)    if n_components > min(X.shape):        raise ValueError("n_components cannot exceed min(n_samples, n_features).")    model = make_plsda_pipeline(n_components)    model.fit(X, y)    return model  def predict_plsda_response(model, X_new):    X_new = np.asarray(X_new, dtype=float)    if X_new.ndim != 2:        raise ValueError("X_new must be a 2D array of samples by variables.")    if not np.all(np.isfinite(X_new)):        raise ValueError("X_new must contain only finite values.")    return np.asarray(model.predict(X_new), dtype=float).reshape(-1)  def predict_plsda_class(model, X_new, threshold=0.5):    y_hat = predict_plsda_response(model, X_new)    return assign_binary_midpoint(y_hat, threshold=threshold), y_hat  def confusion_binary(y_true, y_class):    y_true = np.asarray(y_true, dtype=int).reshape(-1)    y_class = np.asarray(y_class, dtype=int).reshape(-1)    if y_true.shape[0] != y_class.shape[0]:        raise ValueError("y_true and y_class must have the same length.")    tp = int(np.sum((y_true == 1) & (y_class == 1)))    tn = int(np.sum((y_true == 0) & (y_class == 0)))    fp = int(np.sum((y_true == 0) & (y_class == 1)))    fn = int(np.sum((y_true == 1) & (y_class == 0)))    return {"tp": tp, "tn": tn, "fp": fp, "fn": fn}  def sensitivity(cm):    den = cm["tp"] + cm["fn"]    if den == 0:        raise ValueError("Sensitivity is undefined when there are no positive samples.")    return cm["tp"] / den  def specificity(cm):    den = cm["tn"] + cm["fp"]    if den == 0:        raise ValueError("Specificity is undefined when there are no negative samples.")    return cm["tn"] / den  def nested_oof_plsda(    X,    y,    a_max=8,    outer_splits=5,    inner_splits=4,    random_state=0,    threshold=0.5,):    X, y = _validate_xy(X, y)    y_int = y.astype(int)    outer = StratifiedKFold(        n_splits=outer_splits, shuffle=True, random_state=random_state    )    inner = StratifiedKFold(        n_splits=inner_splits, shuffle=True, random_state=random_state    )    y_hat = np.full(y.shape[0], np.nan)    a_star = np.zeros(y.shape[0], dtype=int)    folds = []    for train, test in outer.split(X, y_int):        a_limit = min(int(a_max), train.shape[0] - 1, X.shape[1])        best_a = 1        best_nmc = np.inf        for a in range(1, a_limit + 1):            oof = np.full(train.shape[0], np.nan)            X_tr, y_tr = X[train], y[train]            y_tr_int = y_int[train]            for inner_tr, inner_va in inner.split(X_tr, y_tr_int):                model = make_plsda_pipeline(a)                model.fit(X_tr[inner_tr], y_tr[inner_tr])                oof[inner_va] = predict_plsda_response(model, X_tr[inner_va])            nmc = int(np.sum(y_tr_int != assign_binary_midpoint(oof, threshold)))            if nmc < best_nmc:                best_nmc = nmc                best_a = a        model = make_plsda_pipeline(best_a)        model.fit(X[train], y[train])        y_hat[test] = predict_plsda_response(model, X[test])        a_star[test] = best_a        folds.append((train.copy(), test.copy()))    y_class = assign_binary_midpoint(y_hat, threshold)    scaler_means = []    for train, _test in folds:        scaler_means.append(X[train].mean(axis=0))    return {        "y_hat": y_hat,        "y_class": y_class,        "a_star": a_star,        "folds": folds,        "scaler_means": np.vstack(scaler_means),    } 

1. 1.

   Barker, M., & Rayens, W. (2003). Partial least squares for discrimination. Journal of Chemometrics, 17(3), 166-173.

   doi:10.1002/cem.785
2. 2.

   Westerhuis, J. A., Hoefsloot, H. C. J., Smit, S., Vis, D. J., Smilde, A. K., van Velzen, E. J. J., van Duijnhoven, J. P. M., & van Dorsten, F. A. (2008). Assessment of PLSDA cross validation. Metabolomics, 4, 81-89.

   doi:10.1007/s11306-007-0099-6
3. 3.

   Szymańska, E., Saccenti, E., Smilde, A. K., & Westerhuis, J. A. (2012). Double-check: validation of diagnostic statistics for PLS-DA models in metabolomics studies. Metabolomics, 8(Suppl 1), 3-16.

   doi:10.1007/s11306-011-0330-3
4. 4.

   Altman, D. G., & Bland, J. M. (1994). Diagnostic tests. 1: Sensitivity and specificity. BMJ, 308(6943), 1552.

   doi:10.1136/bmj.308.6943.1552
5. 5.

   scikit-learn Developers (n.d.). PLSRegression, Pipeline, StandardScaler, and LabelBinarizer. scikit-learn 1.6 documentation.

6. 6.

   The MathWorks, Inc. (n.d.). plsregress. MATLAB documentation.