PLS Regression: Zero to Hero Tutorial

This comprehensive tutorial takes you from the foundational concepts of Partial Least Squares (PLS) Regression all the way through advanced component interpretation, model validation, and practical usage within the DataStatPro application. Whether you are a complete beginner or an experienced analyst, this guide is structured to build your understanding step by step.

Table of Contents

1. Prerequisites and Background Concepts
2. What is PLS Regression?
3. The Mathematics Behind PLS Regression
4. Types of PLS Methods
5. Assumptions of PLS Regression
6. Data Preprocessing
7. PLS Components and Latent Variables
8. Choosing the Number of Components
9. Model Fit and Evaluation
10. Interpretation of PLS Results
11. Validation Methods
12. Comparison with Related Methods
13. Using the PLS Regression Component
14. Computational and Formula Details
15. Worked Examples
16. Common Mistakes and How to Avoid Them
17. Troubleshooting
18. Quick Reference Cheat Sheet

1. Prerequisites and Background Concepts

Before diving into PLS regression, it is helpful to be familiar with the following foundational concepts. Do not worry if you are not — each concept is briefly explained here.

1.1 Vectors and Matrices

A vector is an ordered list of numbers (a one-dimensional array). A matrix is a two-dimensional array of numbers with $n$ rows and $p$ columns, denoted $\mathbf{X}_{n \times p}$.

Key matrix operations used throughout this tutorial:

• Transpose: $\mathbf{X}^T$ swaps rows and columns of $\mathbf{X}$.
• Matrix multiplication: $\mathbf{A}_{m \times k} \mathbf{B}_{k \times n} = \mathbf{C}_{m \times n}$.
• Matrix inverse: $\mathbf{A}^{-1}$ exists only if $\mathbf{A}$ is square and non-singular.
• Dot product: $\mathbf{a}^T \mathbf{b} = \sum_{i} a_i b_i$ (a scalar).
• Norm: $\|\mathbf{a}\| = \sqrt{\mathbf{a}^T \mathbf{a}} = \sqrt{\sum_i a_i^2}$ (the length of a vector).

1.2 Projection and Orthogonality

The projection of vector $\mathbf{y}$ onto vector $\mathbf{t}$ is:

$\text{proj}_{\mathbf{t}} \mathbf{y} = \frac{\mathbf{t}^T \mathbf{y}}{\mathbf{t}^T \mathbf{t}} \mathbf{t}$

Two vectors are orthogonal if their dot product equals zero: $\mathbf{a}^T \mathbf{b} = 0$. Orthogonality is a central property of PLS components.

1.3 Variance and Covariance

The variance of a variable $X$ measures its spread:

$\text{Var}(X) = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar{x})^2$

The covariance between variables $X$ and $Y$ measures how they vary together:

$\text{Cov}(X, Y) = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})$

The covariance matrix $\mathbf{S}_{p \times p}$ of a matrix $\mathbf{X}$ contains pairwise covariances of all $p$ variables. After mean-centring $\mathbf{X}$:

$\mathbf{S} = \frac{1}{n-1}\mathbf{X}^T \mathbf{X}$

1.4 Correlation and Multicollinearity

Correlation is the standardised covariance:

$r_{XY} = \frac{\text{Cov}(X, Y)}{\sqrt{\text{Var}(X) \cdot \text{Var}(Y)}}$

Multicollinearity occurs when predictor variables are highly correlated with each other. It causes serious problems in ordinary least squares (OLS) regression:

• Standard errors of coefficients become very large.
• Coefficient estimates become unstable and unreliable.
• The matrix $\mathbf{X}^T \mathbf{X}$ becomes nearly singular (non-invertible).

PLS regression was specifically designed to handle multicollinearity.

1.5 Eigenvalues and Eigenvectors

For a square matrix $\mathbf{A}$, an eigenvector $\mathbf{v}$ and its associated eigenvalue $\lambda$ satisfy:

$\mathbf{A}\mathbf{v} = \lambda \mathbf{v}$

Eigenvalues indicate the amount of variance explained by each eigenvector direction. They are the foundation of Principal Component Analysis (PCA) and are closely related to PLS decompositions.

1.6 Ordinary Least Squares (OLS) Regression

OLS regression models the response $\mathbf{y}$ as a linear function of predictors $\mathbf{X}$:

$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}$

The OLS estimator minimises the sum of squared residuals:

$\hat{\boldsymbol{\beta}}_{OLS} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}$

This requires $\mathbf{X}^T \mathbf{X}$ to be invertible — which fails when:

• $p > n$ (more predictors than observations).
• Predictors are highly collinear.

PLS provides a solution to both of these failure modes.

2. What is PLS Regression?

Partial Least Squares (PLS) Regression is a multivariate statistical method that models the relationship between a matrix of predictor variables $\mathbf{X}$ and one or more response variables $\mathbf{Y}$ by extracting a small number of latent components (factors) that simultaneously:

1. Explain as much variance in $\mathbf{X}$ as possible.
2. Have maximum covariance (predictive power) with $\mathbf{Y}$.

Unlike OLS regression, PLS does not require $\mathbf{X}^T \mathbf{X}$ to be invertible and performs well even when predictors far outnumber observations or are highly collinear.

2.1 The Core Idea

The name "Partial Least Squares" reflects the method's origins:

• Partial: Each step involves projecting (regressing) the data onto a component that captures part of the variation in $\mathbf{X}$ and $\mathbf{Y}$.
• Least Squares: At each step, the component is found to minimise a least-squares criterion.

In essence, PLS finds a compressed, low-dimensional representation of $\mathbf{X}$ (the latent components or scores) that is maximally informative about $\mathbf{Y}$. This is what distinguishes PLS from Principal Component Regression (PCR), which only maximises variance in $\mathbf{X}$ without considering $\mathbf{Y}$.

2.2 Real-World Applications

PLS regression is the workhorse of many applied sciences. Common applications include:

• Chemometrics: Predicting chemical concentrations (e.g., protein content, moisture, pH) from near-infrared (NIR) or Raman spectroscopy data — the quintessential PLS application with hundreds or thousands of correlated spectral wavelengths as predictors.
• Genomics & Bioinformatics: Relating gene expression profiles (tens of thousands of genes) to clinical outcomes, phenotypes, or disease states.
• Food Science: Predicting sensory quality attributes (taste, texture) from instrumental measurements.
• Pharmaceutical Sciences: Relating molecular descriptors (QSAR/QSPR) to drug activity, toxicity, or physicochemical properties.
• Process Monitoring & Control: Modelling industrial process variables to predict product quality and detect deviations.
• Neuroscience: Relating brain imaging data (fMRI voxels) to behavioural or cognitive measures.
• Social & Behavioural Sciences: Predicting complex outcomes (job performance, wellbeing) from many correlated survey items.
• Environmental Science: Relating environmental variables (soil, water, climate) to ecological outcomes.

2.3 When to Use PLS Regression

PLS regression is particularly appropriate when:

2.4 PLS vs. Related Methods: An Overview

3. The Mathematics Behind PLS Regression

3.1 The PLS Model Structure

PLS simultaneously decomposes both $\mathbf{X}$ ($n \times p$) and $\mathbf{Y}$ ($n \times q$) into score matrices and loading matrices:

$\mathbf{X} = \mathbf{T}\mathbf{P}^T + \mathbf{E}$

$\mathbf{Y} = \mathbf{U}\mathbf{Q}^T + \mathbf{F}$

Where:

• $\mathbf{T}$ ($n \times A$) = X-scores matrix (latent components from $\mathbf{X}$).
• $\mathbf{P}$ ($p \times A$) = X-loadings matrix (how $\mathbf{X}$ variables relate to $\mathbf{X}$-components).
• $\mathbf{U}$ ($n \times A$) = Y-scores matrix (latent components from $\mathbf{Y}$).
• $\mathbf{Q}$ ($q \times A$) = Y-loadings matrix (how $\mathbf{Y}$ variables relate to $\mathbf{Y}$-components).
• $\mathbf{E}$ ($n \times p$) = X-residual matrix.
• $\mathbf{F}$ ($n \times q$) = Y-residual matrix.
• $A$ = number of retained latent components.

The relationship between the $\mathbf{X}$-scores and $\mathbf{Y}$-scores is modelled as:

$\mathbf{U} = \mathbf{T}\mathbf{B}_{diag} + \mathbf{H}$

Where $\mathbf{B}_{diag}$ is a diagonal matrix of inner relation coefficients and $\mathbf{H}$ is residual.

3.2 The Objective: Maximum Covariance

The core objective of PLS is to find weight vectors $\mathbf{w}$ (for $\mathbf{X}$) and $\mathbf{c}$ (for $\mathbf{Y}$) such that the covariance between the resulting scores is maximised:

$\max_{\mathbf{w}, \mathbf{c}} \text{Cov}(\mathbf{X}\mathbf{w},\ \mathbf{Y}\mathbf{c})^2 = \max_{\mathbf{w}, \mathbf{c}} (\mathbf{w}^T \mathbf{X}^T \mathbf{Y} \mathbf{c})^2$

Subject to the normalisation constraints: $\|\mathbf{w}\| = 1$ and $\|\mathbf{c}\| = 1$.

This is equivalent to finding the first singular value decomposition (SVD) of the cross-product matrix $\mathbf{X}^T \mathbf{Y}$:

$\mathbf{X}^T \mathbf{Y} = \mathbf{W} \mathbf{S} \mathbf{C}^T$

Where $\mathbf{S}$ is diagonal with singular values (covariances), $\mathbf{W}$ contains X-weight vectors, and $\mathbf{C}$ contains Y-weight vectors.

3.3 Score and Loading Computation

For each component $a = 1, 2, \dots, A$:

X-scores (latent variable scores):

$\mathbf{t}_a = \mathbf{X}_a \mathbf{w}_a$

Where $\mathbf{X}_a$ is the deflated (residualised) $\mathbf{X}$ matrix at step $a$, and $\mathbf{w}_a$ is the X-weight vector for component $a$.

Y-scores:

$\mathbf{u}_a = \mathbf{Y}_a \mathbf{c}_a$

X-loadings (regression of $\mathbf{X}_a$ on $\mathbf{t}_a$):

$\mathbf{p}_a = \frac{\mathbf{X}_a^T \mathbf{t}_a}{\mathbf{t}_a^T \mathbf{t}_a}$

Y-loadings (regression of $\mathbf{Y}_a$ on $\mathbf{u}_a$):

$\mathbf{q}_a = \frac{\mathbf{Y}_a^T \mathbf{u}_a}{\mathbf{u}_a^T \mathbf{u}_a}$

Inner relation coefficient:

$b_a = \frac{\mathbf{u}_a^T \mathbf{t}_a}{\mathbf{t}_a^T \mathbf{t}_a}$

3.4 Deflation

After extracting each component, the data matrices are deflated (the component's contribution is removed) to ensure that subsequent components capture new, orthogonal information:

$\mathbf{X}_{a+1} = \mathbf{X}_a - \mathbf{t}_a \mathbf{p}_a^T$

For PLS1 (single response):

$\mathbf{y}_{a+1} = \mathbf{y}_a - b_a \mathbf{t}_a$

For PLS2 (multiple responses):

$\mathbf{Y}_{a+1} = \mathbf{Y}_a - b_a \mathbf{t}_a \mathbf{q}_a^T$

This sequential deflation ensures the X-scores $\mathbf{t}_1, \mathbf{t}_2, \dots, \mathbf{t}_A$ are mutually orthogonal: $\mathbf{t}_a^T \mathbf{t}_b = 0$ for $a \neq b$.

3.5 The PLS Regression Coefficients

After extracting $A$ components, the final PLS regression coefficients relating $\mathbf{X}$ directly to $\mathbf{Y}$ are:

$\hat{\mathbf{B}}_{PLS} = \mathbf{W}^* \mathbf{Q}^T$

Where $\mathbf{W}^*$ ($p \times A$) is the matrix of modified weight vectors (also called the W-star or $\mathbf{R}$ matrix), which accounts for the sequential deflation:

$\mathbf{W}^* = \mathbf{W}(\mathbf{P}^T \mathbf{W})^{-1}$

The predicted values are then:

$\hat{\mathbf{Y}} = \mathbf{X} \hat{\mathbf{B}}_{PLS} = \mathbf{X} \mathbf{W}^* \mathbf{Q}^T = \mathbf{T} \mathbf{Q}^T$

Where $\mathbf{T} = \mathbf{X} \mathbf{W}^*$ are the X-scores computed directly from $\mathbf{X}$ without deflation.

For a single response variable $\mathbf{y}$:

$\hat{\mathbf{y}} = \mathbf{X} \hat{\boldsymbol{\beta}}_{PLS}, \quad \hat{\boldsymbol{\beta}}_{PLS} = \mathbf{W}^* \mathbf{q}$

The predicted value for a new observation $\mathbf{x}_{new}$:

$\hat{y}_{new} = \mathbf{x}_{new}^T \hat{\boldsymbol{\beta}}_{PLS} = \mathbf{x}_{new}^T \mathbf{W}^* \mathbf{q}$

3.6 Relationship Between PLS and SVD

The core step of PLS — finding weight vectors that maximise the covariance between $\mathbf{t} = \mathbf{X}\mathbf{w}$ and $\mathbf{u} = \mathbf{Y}\mathbf{c}$ — reduces to finding the leading left and right singular vectors of $\mathbf{X}^T \mathbf{Y}$:

$\mathbf{X}^T \mathbf{Y} \mathbf{Y}^T \mathbf{X} \mathbf{w} = \lambda^2 \mathbf{w}$

$\mathbf{Y}^T \mathbf{X} \mathbf{X}^T \mathbf{Y} \mathbf{c} = \lambda^2 \mathbf{c}$

The maximum covariance $\lambda$ is the first singular value of $\mathbf{X}^T \mathbf{Y}$, and $\mathbf{w}$, $\mathbf{c}$ are the corresponding left and right singular vectors.

4. Types of PLS Methods

Several algorithmic variants of PLS exist. The most important ones are described below.

4.1 PLS1

PLS1 handles a single continuous response variable $\mathbf{y}$ ($q = 1$). It is the most common form of PLS regression in practice.

• Input: $\mathbf{X}$ ($n \times p$) and $\mathbf{y}$ ($n \times 1$).
• Output: $A$ latent components, regression coefficients $\hat{\boldsymbol{\beta}}_{PLS}$ ($p \times 1$).
• The NIPALS algorithm (see Section 14) is used for computation.

4.2 PLS2

PLS2 handles multiple response variables simultaneously ($q > 1$). It extracts components that explain $\mathbf{X}$ variance while maximising covariance with the entire $\mathbf{Y}$ matrix.

• Input: $\mathbf{X}$ ($n \times p$) and $\mathbf{Y}$ ($n \times q$).
• Output: $A$ latent components, regression coefficient matrix $\hat{\mathbf{B}}_{PLS}$ ($p \times q$).
• A single set of X-components is found for all responses simultaneously.

💡 PLS2 is more parsimonious than running separate PLS1 models for each response, but if the responses are very different in nature, separate PLS1 models may give better individual predictions.

4.3 PLS-DA (PLS Discriminant Analysis)

PLS-DA applies PLS regression to a categorical response variable by encoding the class membership as a binary (or dummy-coded) $\mathbf{Y}$ matrix. It is widely used for classification in chemometrics, metabolomics, and genomics.

• For a two-class problem: Encode $\mathbf{y} \in \{0, 1\}$ and apply PLS1. Classify using a threshold on $\hat{y}$.
• For a multi-class problem ($K$ classes): Encode as a $K$-column binary matrix (one-hot) and apply PLS2. Assign to the class with the highest predicted value.

⚠️ PLS-DA can overfit when many components are used, especially with small or imbalanced datasets. Cross-validation is essential for assessing classification performance.

4.4 OPLS (Orthogonal PLS)

OPLS (Trygg & Wold) separates the $\mathbf{X}$ variation into:

• Predictive variation: Correlated with $\mathbf{Y}$ (the PLS component).
• Orthogonal variation: Uncorrelated with $\mathbf{Y}$ (structured noise in $\mathbf{X}$).

This produces a model with a single predictive component (for PLS1) plus orthogonal components that account for systematic $\mathbf{X}$ variation not related to $\mathbf{Y}$. OPLS produces simpler, more interpretable loading plots.

4.5 Kernel PLS

Kernel PLS extends PLS to handle non-linear relationships between $\mathbf{X}$ and $\mathbf{Y}$ by implicitly mapping the data into a high-dimensional feature space using a kernel function $k(\mathbf{x}_i, \mathbf{x}_j)$.

Common kernels:

• Linear: $k(\mathbf{x}_i, \mathbf{x}_j) = \mathbf{x}_i^T \mathbf{x}_j$
• Polynomial: $k(\mathbf{x}_i, \mathbf{x}_j) = (\mathbf{x}_i^T \mathbf{x}_j + c)^d$
• Radial Basis Function (RBF): $k(\mathbf{x}_i, \mathbf{x}_j) = \exp(-\gamma \|\mathbf{x}_i - \mathbf{x}_j\|^2)$

4.6 Sparse PLS

Sparse PLS introduces variable selection by imposing $L_1$ (Lasso-type) penalties on the weight vectors, driving some weights to exactly zero. This simultaneously performs dimensionality reduction and variable selection, producing more interpretable models in high-dimensional settings ($p \gg n$).

4.7 Summary of PLS Variants

The DataStatPro application implements PLS1 and PLS2 (standard PLS regression) and PLS-DA, which are the focus of this tutorial.

5. Assumptions of PLS Regression

PLS regression is a relatively assumption-light method compared to OLS. However, certain conditions should be met for valid results.

5.1 Linearity

PLS assumes a linear relationship between the latent components (scores $\mathbf{T}$) and the response $\mathbf{Y}$. Non-linear relationships between the original $\mathbf{X}$ variables and $\mathbf{Y}$ may be partially captured if the non-linearity is reflected in the latent structure, but Kernel PLS is preferred for strongly non-linear data.

5.2 Continuous (or Appropriately Encoded) Variables

• Predictor variables ($\mathbf{X}$): Should be continuous or ordinal. Binary variables can be included but must be mean-centred and scaled.
• Response variable ($\mathbf{y}$): Should be continuous for PLS1/PLS2. For categorical responses, use PLS-DA with appropriate dummy coding.

5.3 No Requirement for Multivariate Normality

Unlike some classical multivariate methods, PLS does not require the predictors or response to follow a multivariate normal distribution. This makes PLS robust to skewed or non-normal variables (though severe non-normality may still affect inference).

5.4 No Perfect Redundancy (Degenerate Cases)

If a predictor variable $X_j$ is a perfect linear combination of other predictor variables (i.e., the column is a deterministic function of others), it carries no additional information and should be removed before fitting PLS. Similarly, a response variable that is a perfect linear combination of $\mathbf{X}$ columns will cause a degenerate model.

5.5 Sufficient Sample Size

While PLS handles $p > n$ settings, model quality improves with more observations. General guidelines:

• Minimum: $n \geq 5A$ where $A$ is the number of components (to avoid severe overfitting).
• Preferred: $n \geq 10A$ for reliable cross-validated results.
• For PLS-DA: At least 5–10 observations per class per component.

5.6 Representativeness of Calibration Set

The calibration (training) samples should span the full range of variation expected in future prediction samples. PLS is an interpolation method — predictions for samples outside the calibration space (extrapolation) are unreliable.

5.7 No Gross Outliers

Extreme outliers in $\mathbf{X}$ or $\mathbf{Y}$ can disproportionately influence the extracted components and distort the model. Outliers should be detected (using score plots and leverage/residual diagnostics) and investigated before finalising the model.

6. Data Preprocessing

Data preprocessing is arguably the most critical step in PLS analysis. The choice of preprocessing can have a profound impact on the extracted components and the resulting model.

6.1 Mean Centring

Mean centring subtracts the column mean from each variable:

$\tilde{x}_{ij} = x_{ij} - \bar{x}_j, \quad \bar{x}_j = \frac{1}{n}\sum_{i=1}^n x_{ij}$

Mean centring is almost always required for PLS. Without it, the first component tends to describe the mean of the data rather than the variance structure, and subsequent components are distorted.

After mean centring, the PLS model is fitted to $\tilde{\mathbf{X}}$ and $\tilde{\mathbf{Y}}$, and predictions are adjusted using the response mean:

$\hat{y}_{new} = \bar{y} + \tilde{\mathbf{x}}_{new}^T \hat{\boldsymbol{\beta}}_{PLS}$

6.2 Autoscaling (Mean Centring + Unit Variance Scaling)

Autoscaling (also called standardisation or z-score scaling) both mean-centres and scales each variable to unit variance:

$\tilde{x}_{ij} = \frac{x_{ij} - \bar{x}_j}{s_j}, \quad s_j = \sqrt{\frac{1}{n-1}\sum_{i=1}^n (x_{ij} - \bar{x}_j)^2}$

When to use autoscaling:

• When predictor variables are measured in different units (e.g., age in years, income in dollars, height in centimetres). Without scaling, variables with larger numerical ranges will dominate the components.
• As the default choice when there is no strong reason to prefer otherwise.

When not to autoscale:

• When variables are measured in the same units and differences in variance are meaningful (e.g., spectroscopic data where high-variance spectral regions are genuinely more informative).
• When some variables are near-zero variance (noise-dominated) — autoscaling would amplify noise.

6.3 Other Scaling Methods

💡 The response variable $\mathbf{y}$ should also be mean-centred (and scaled if using PLS2 with multiple responses on different scales). For PLS1, scaling $\mathbf{y}$ to unit variance is optional but often beneficial.

6.4 Handling Missing Data

PLS with missing data requires special treatment. Options include:

• Complete case analysis: Remove observations with any missing values (only if missingness is minimal and random).
• Mean imputation: Replace missing values with the column mean (simple but ignores covariance structure).
• NIPALS-based imputation: The NIPALS algorithm can handle missing data iteratively — missing values are replaced with their model estimates at each iteration.
• Multiple imputation: Generate multiple complete datasets and pool the PLS results.

⚠️ Missing data in $\mathbf{X}$ is more tractable than missing data in $\mathbf{Y}$. If $\mathbf{y}$ values are missing, those observations typically cannot be used in model calibration.

6.5 Outlier Detection Before Modelling

Before fitting PLS, screen for outliers using:

• Univariate checks: Box plots, histograms, z-scores ($|z| > 3$ as a flag).
• Multivariate checks: Mahalanobis distance, PCA score plots.
• Domain knowledge: Values that are physically impossible (negative concentrations, ages > 150) should be corrected or removed.

7. PLS Components and Latent Variables

Understanding what PLS components represent is essential for interpreting PLS models.

7.1 X-Scores ($\mathbf{T}$)

The X-scores matrix $\mathbf{T}$ ($n \times A$) contains the coordinates of each observation in the low-dimensional latent space:

$\mathbf{t}_a = \mathbf{X} \mathbf{w}_a^*$

Where $\mathbf{w}_a^*$ is the $a$-th column of $\mathbf{W}^*$.

• Each row $\mathbf{t}_i$ is the score vector for observation $i$ — its position in the latent space.
• Scores are mutually orthogonal: $\mathbf{t}_a^T \mathbf{t}_b = 0$ for $a \neq b$.
• Score plots ($t_1$ vs. $t_2$, etc.) reveal the structure of the observations: clusters, trends, outliers.

7.2 X-Loadings ($\mathbf{P}$)

The X-loadings matrix $\mathbf{P}$ ($p \times A$) describes how the original $\mathbf{X}$ variables contribute to each latent component during deflation:

$\mathbf{p}_a = \frac{\mathbf{X}_a^T \mathbf{t}_a}{\mathbf{t}_a^T \mathbf{t}_a}$

• Loading $p_{ja}$ indicates how strongly variable $X_j$ is associated with component $a$.
• Loading plots reveal which $\mathbf{X}$ variables are most important for each component and which variables are correlated (variables with similar loading vectors are collinear).

7.3 X-Weights ($\mathbf{W}$ and $\mathbf{W}^*$)

The X-weights $\mathbf{W}$ describe how $\mathbf{X}$ variables are weighted to form the X-scores:

$\mathbf{t}_a = \mathbf{X}_a \mathbf{w}_a \quad \text{(during deflation)}$

The modified X-weights $\mathbf{W}^* = \mathbf{W}(\mathbf{P}^T \mathbf{W})^{-1}$ relate the original (undeflated) $\mathbf{X}$ directly to the scores:

$\mathbf{T} = \mathbf{X} \mathbf{W}^*$

⚠️ There is a subtle but important distinction between $\mathbf{W}$ and $\mathbf{W}^*$. For interpreting the relationship between $\mathbf{X}$ variables and the PLS components, use $\mathbf{W}^*$ (not $\mathbf{W}$), as $\mathbf{W}^*$ accounts for the deflation steps and relates to the original $\mathbf{X}$ directly.

7.4 Y-Loadings ($\mathbf{Q}$) and Y-Weights ($\mathbf{C}$)

The Y-loadings $\mathbf{Q}$ ($q \times A$) describe how the $\mathbf{Y}$ variables are reconstructed from the X-scores:

$\hat{\mathbf{Y}} = \mathbf{T} \mathbf{Q}^T$

Loading $q_{ja}$ indicates how strongly response variable $Y_j$ is associated with component $a$ in the final prediction equation. For PLS1 ($q = 1$), $\mathbf{Q}$ reduces to a scalar $q_a$ for each component.

7.5 Y-Scores ($\mathbf{U}$)

The Y-scores $\mathbf{U}$ ($n \times A$) are the latent components extracted from $\mathbf{Y}$:

$\mathbf{u}_a = \mathbf{Y}_a \mathbf{c}_a$

In the inner relation $\mathbf{u}_a \approx b_a \mathbf{t}_a$, the Y-scores should be close to the X-scores (scaled by $b_a$). The tightness of this relationship is diagnostic of model quality:

• If $\mathbf{t}_a$ and $\mathbf{u}_a$ are strongly correlated, the component is predictively powerful.
• If they are weakly correlated, the component explains $\mathbf{X}$ variance but not $\mathbf{Y}$ — suggesting the component is not useful for prediction.

7.6 The Biplot

The PLS biplot overlays the score plot (observations) and loading plot (variables) in the same space:

• Observations (scores $t_1, t_2$) are plotted as points.
• Variables (loadings $p_{a1}, p_{a2}$ or weights $w_{a1}^*, w_{a2}^*$) are plotted as vectors.
• A variable vector pointing toward a cluster of observations indicates those observations have high values for that variable.
• Variables with loading vectors in the same direction are positively correlated; opposite directions indicate negative correlation.

8. Choosing the Number of Components

The number of latent components $A$ is the primary hyperparameter in PLS regression. Too few components underfits (high bias, misses important structure); too many overfits (low bias but high variance, memorises noise).

8.1 Cross-Validation (Primary Method)

Cross-validation (CV) is the gold-standard method for selecting $A$ in PLS. The most common approach is $k$-fold cross-validation:

1. Divide the $n$ observations into $k$ roughly equal folds.
2. For each fold $v = 1, \dots, k$: a. Fit PLS with $a = 1, 2, \dots, A_{\max}$ components on the data excluding fold $v$. b. Predict fold $v$ observations using each fitted model. c. Compute prediction errors for fold $v$.
3. Average the prediction errors across all folds.
4. Select $A$ minimising the cross-validated prediction error.

Leave-One-Out Cross-Validation (LOOCV): Special case of $k$-fold CV with $k = n$. Computationally intensive but uses maximum data for fitting at each step.

Recommended $k$: $k = 5$ or $k = 10$ is standard. For small datasets ($n < 30$), LOOCV is preferred.

8.2 PRESS Statistic (Predicted Residual Error Sum of Squares)

The PRESS statistic summarises the cross-validated prediction error for $A$ components:

$PRESS(A) = \sum_{v=1}^k \sum_{i \in \text{fold } v} \left(y_i - \hat{y}_{i,-v}(A)\right)^2$

Where $\hat{y}_{i,-v}(A)$ is the prediction for observation $i$ when it was in the held-out fold $v$, using a model with $A$ components.

Select $A^* = \arg\min_A PRESS(A)$.

💡 To guard against overfitting while maintaining parsimony, some guidelines recommend choosing the smallest $A$ for which $PRESS(A)$ is within 1 standard error of the minimum PRESS (the "one-standard-error rule").

8.3 $Q^2$ (Cross-Validated $R^2$)

$Q^2$ is the cross-validated analogue of $R^2$:

$Q^2(A) = 1 - \frac{PRESS(A)}{SS_{Y,total}}$

Where:

$SS_{Y,total} = \sum_{i=1}^n (y_i - \bar{y})^2$

$Q^2$ ranges from $-\infty$ (model is worse than predicting the mean) to 1 (perfect cross-validated prediction).

Guidelines for $Q^2$:

The optimal number of components is where $Q^2$ reaches a maximum (or where adding more components does not meaningfully increase $Q^2$).

8.4 The $R^2_Y$ vs. $Q^2$ Plot

A standard diagnostic in PLS is to plot both $R^2_Y$ (variance explained in $\mathbf{Y}$ — a training set metric) and $Q^2$ (cross-validated metric) against the number of components $A$:

• $R^2_Y$ always increases with more components (never decreases).
• $Q^2$ first increases then decreases (or plateaus) as the model starts to overfit.
• The optimal $A$ is where $Q^2$ peaks or where the gap between $R^2_Y$ and $Q^2$ begins to widen rapidly.

⚠️ A model with $R^2_Y$ much higher than $Q^2$ is overfitting. Aim for $Q^2$ close to $R^2_Y$.

8.5 Scree Plot of Eigenvalues / Variance Explained

A scree plot of the variance explained in $\mathbf{X}$ ($R^2_X$) by each successive component can also guide component selection: look for an "elbow" where additional components explain diminishing amounts of variance. However, $R^2_X$ alone ignores predictive relevance for $\mathbf{Y}$ — always prioritise $Q^2$.

8.6 Permutation Testing

Permutation tests provide a rigorous null-hypothesis test for whether a PLS model with $A$ components explains more variance in $\mathbf{Y}$ than expected by chance:

1. Randomly permute (shuffle) the $\mathbf{y}$ vector $B$ times (e.g., $B = 999$).
2. Fit PLS with $A$ components to each permuted dataset and record $R^2_Y$ and $Q^2$.
3. The p-value is the proportion of permuted $R^2_Y$ (or $Q^2$) values that exceed the observed value.

A significant result ($p < 0.05$) confirms the model captures a real relationship, not a spurious one.

9. Model Fit and Evaluation

9.1 $R^2_X$ (Variance Explained in $\mathbf{X}$)

The proportion of total $\mathbf{X}$ variance explained by component $a$ is:

$R^2_{X,a} = \frac{\|\mathbf{t}_a \mathbf{p}_a^T\|^2}{\|\mathbf{X}\|^2} = \frac{\mathbf{t}_a^T \mathbf{t}_a \cdot \mathbf{p}_a^T \mathbf{p}_a}{\|\mathbf{X}\|^2}$

Cumulative $R^2_X$ after $A$ components:

$R^2_X(A) = \sum_{a=1}^A R^2_{X,a}$

A high $R^2_X$ indicates the components capture most of the systematic variation in $\mathbf{X}$.

9.2 $R^2_Y$ (Variance Explained in $\mathbf{Y}$)

The proportion of total $\mathbf{Y}$ variance explained by $A$ components — the training set $R^2$:

$R^2_Y(A) = 1 - \frac{\|\mathbf{Y} - \hat{\mathbf{Y}}\|^2}{\|\mathbf{Y} - \bar{\mathbf{Y}}\|^2} = 1 - \frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{\sum_{i=1}^n (y_i - \bar{y})^2}$

For PLS1, this reduces to:

$R^2_Y(A) = 1 - \frac{RSS(A)}{SS_{total}}$

Where $RSS(A) = \sum_{i=1}^n (y_i - \hat{y}_i)^2$ is the residual sum of squares.

⚠️ $R^2_Y$ alone is an optimistic (biased) estimate of model quality because it is computed on the training data. Always report $Q^2$ alongside $R^2_Y$.

9.3 Root Mean Squared Error of Calibration (RMSEC)

$RMSEC = \sqrt{\frac{1}{n}\sum_{i=1}^n (y_i - \hat{y}_i)^2}$

This is the training set prediction error in the same units as $\mathbf{y}$.

9.4 Root Mean Squared Error of Cross-Validation (RMSECV)

$RMSECV = \sqrt{\frac{PRESS}{n}} = \sqrt{\frac{1}{n}\sum_{i=1}^n (y_i - \hat{y}_{i,-v})^2}$

RMSECV is the cross-validated prediction error and is the primary model selection criterion alongside $Q^2$.

$Q^2 = 1 - \frac{PRESS}{SS_{total}} = 1 - \frac{n \cdot RMSECV^2}{SS_{total}}$

9.5 Root Mean Squared Error of Prediction (RMSEP)

When a separate independent test set of $n_{test}$ observations is available (not used in model fitting or cross-validation):

$RMSEP = \sqrt{\frac{1}{n_{test}}\sum_{i=1}^{n_{test}} (y_i - \hat{y}_i)^2}$

RMSEP is the most honest estimate of true prediction error and should always be reported when a genuine external test set exists.

Hierarchy of prediction error estimates:

$RMSEC \leq RMSECV \leq RMSEP$

The gap between RMSEC and RMSECV/RMSEP indicates the degree of overfitting.

9.6 Bias and Slope of Predicted vs. Observed

A predicted vs. observed plot should ideally show points scattered symmetrically around the 1:1 line (slope = 1, intercept = 0). Formally test for systematic bias using a regression of observed on predicted:

$y_i = b_0 + b_1 \hat{y}_i + \epsilon_i$

• $b_0 \neq 0$: Systematic bias (model consistently over- or under-predicts).
• $b_1 \neq 1$: Scale bias (model predictions are proportionally stretched or compressed).

9.7 Summary of Model Fit Statistics

10. Interpretation of PLS Results

10.1 Regression Coefficients ($\hat{\boldsymbol{\beta}}_{PLS}$)

The PLS regression coefficients $\hat{\boldsymbol{\beta}}_{PLS} = \mathbf{W}^* \mathbf{q}$ have the same interpretation as OLS coefficients: $\hat{\beta}_j$ is the expected change in $\hat{y}$ for a one-unit increase in $X_j$, with all other variables held constant.

However, because PLS implicitly regularises through dimensionality reduction:

• Coefficients tend to be shrunk toward zero compared to OLS (especially for components not retained).
• They are more stable under multicollinearity than OLS coefficients.
• The number of components $A$ acts as the regularisation parameter: fewer components = more shrinkage.

⚠️ When $\mathbf{X}$ variables are autoscaled, the coefficients are on the standardised scale. Multiply by $s_j$ (the SD of $X_j$) and divide by $s_y$ to interpret on the original scales, or report the raw (unstandardised) coefficients after back-scaling.

10.2 Variable Importance in Projection (VIP)

The VIP score (Wold, 1994) quantifies the contribution of each predictor variable to the PLS model, accounting for all $A$ components:

$VIP_j = \sqrt{\frac{p \sum_{a=1}^A \left(R^2_Y(a) - R^2_Y(a-1)\right) w_{ja}^{*2}}{\sum_{a=1}^A \left(R^2_Y(a) - R^2_Y(a-1)\right) \sum_{j=1}^p w_{ja}^{*2}}}$

Where:

• $p$ = total number of predictor variables.
• $R^2_Y(a) - R^2_Y(a-1)$ = marginal variance explained by component $a$.
• $w_{ja}^{*2}$ = squared modified weight of variable $j$ in component $a$.

Properties of VIP:

• $VIP_j \geq 0$ always.
• The average squared VIP across all variables equals 1: $\frac{1}{p}\sum_{j=1}^p VIP_j^2 = 1$.
• Therefore, the average VIP = 1 (approximately, depending on scaling).

Interpretation:

💡 VIP is the most widely used variable selection criterion in PLS. Variables with VIP < 0.8 (or a domain-specific threshold) are candidates for removal, which can improve model parsimony and sometimes prediction performance.

10.3 Loadings Plot

The loadings plot displays the X-loadings $\mathbf{p}_a$ (or weights $\mathbf{w}_a^*$) for two components simultaneously, revealing:

• Which variables are most strongly associated with each component.
• Which variables are correlated with each other (similar loading vectors).
• The direction and magnitude of each variable's contribution.

For spectroscopic data, a plot of the loadings as a function of wavelength (a "loading spectrum") is particularly informative, as it reveals which spectral regions are important.

10.4 Score Plot

The score plot displays the X-scores ($t_1$ vs. $t_2$, etc.) for all observations, revealing:

• Clustering of observations (groups with similar $\mathbf{X}$ profiles).
• Trends (e.g., samples ordered along $t_1$ by the response value).
• Outliers (isolated points far from the main cluster).

When colour-coded by $\mathbf{y}$ values, the score plot shows whether the latent structure in $\mathbf{X}$ aligns with $\mathbf{Y}$ — the hallmark of a good PLS model.

10.5 Weight-Loading Biplot ($w^*$-$q$ Plot)

The correlation loadings plot or $w^*$-$q$ biplot plots the modified X-weights $\mathbf{w}^*$ and Y-loadings $\mathbf{q}$ together in the same space. Variables (X and Y) that appear close together in this plot are positively correlated; variables on opposite sides are negatively correlated. This provides a comprehensive view of the $\mathbf{X}$–$\mathbf{Y}$ relationship structure.

10.6 Leverage and Residuals (Influence Analysis)

Leverage measures how influential observation $i$ is in determining the model:

$h_i = t_{i\cdot}(\mathbf{T}^T\mathbf{T})^{-1}t_{i\cdot}^T = \frac{1}{n} + \sum_{a=1}^A \frac{t_{ia}^2}{\sum_{j=1}^n t_{ja}^2}$

Where $t_{i\cdot}$ is the $i$-th row of $\mathbf{T}$.

Standardised Y-residuals:

$e_i^* = \frac{y_i - \hat{y}_i}{\hat{\sigma}\sqrt{1 - h_i}}$

A leverage-residual plot (also called a Williams plot) displays $h_i$ vs. $e_i^*$:

• High leverage + small residual: Influential, well-fitting observation (good leverage).
• High leverage + large residual: Outlier with high influence — investigate carefully.
• Low leverage + large residual: Poorly predicted outlier but low influence.
• Low leverage + small residual: Normal observation.

10.7 Hotelling's $T^2$ and SPE (DModX)

Two complementary multivariate control statistics are used to identify unusual samples:

Hotelling's $T^2$: Measures the distance of a sample from the centre of the model in the latent space:

$T^2_i = \mathbf{t}_i^T (\mathbf{T}^T \mathbf{T} / n)^{-1} \mathbf{t}_i = n \sum_{a=1}^A \frac{t_{ia}^2}{\mathbf{t}_a^T \mathbf{t}_a}$

An approximate $F$-distribution critical value:

$T^2_{crit} = \frac{A(n^2-1)}{n(n-A)} F_{\alpha,A,n-A}$

Observations with $T^2_i > T^2_{crit}$ are outliers within the model space (unusual combination of latent components).

SPE (Squared Prediction Error) / DModX: Measures the distance of a sample from the PLS model plane (residual variability not captured by the model):

$SPE_i = \|\mathbf{x}_i - \hat{\mathbf{x}}_i\|^2 = \sum_{j=1}^p (x_{ij} - \hat{x}_{ij})^2$

A high SPE indicates the observation does not conform to the $\mathbf{X}$-covariance structure of the calibration set.

💡 Use both $T^2$ and SPE together. A sample can be unusual inside the model ($T^2$ high) or outside the model (SPE high), or both.

11. Validation Methods

Rigorous validation is essential to ensure a PLS model genuinely predicts new observations, rather than memorising noise in the training data.

11.1 Internal Validation: Cross-Validation

$k$-fold cross-validation and LOOCV (described in Section 8.1) provide internal validation estimates ($Q^2$, RMSECV). Internal validation is mandatory but can still be optimistic if the same data were used to select preprocessing, model type, and components.

Monte Carlo Cross-Validation (MCCV): Randomly splits the data into training and validation sets $B$ times (e.g., $B = 100$), each time using a different random split (e.g., 80% train / 20% validate). The distribution of RMSECV values across splits provides uncertainty estimates.

11.2 External Validation: Independent Test Set

The most rigorous validation uses a truly independent test set — samples not used in any stage of model building (not in training, not in cross-validation, not in component selection):

$RMSEP = \sqrt{\frac{1}{n_{test}}\sum_{i=1}^{n_{test}} (y_i - \hat{y}_i)^2}$

How to split data for external validation:

• Random split: Randomly allocate (e.g., 70–80% training, 20–30% test). Only valid if samples are representative.
• Kennard-Stone algorithm: Selects a maximally representative subset of the data for the training set based on Euclidean distances in $\mathbf{X}$-space, ensuring the calibration set spans the full range of the data.
• DUPLEX algorithm: Simultaneously selects representative training and test sets.
• Systematic split: For time-ordered data, use earlier samples for training and later ones for testing.

11.3 Permutation Test for Model Validity

(Described in Section 8.6.) A permutation test confirms that the model's $Q^2$ is significantly above what would be obtained by chance.

A permutation plot shows the distribution of $Q^2$ values from permuted models overlaid with the observed $Q^2$. If the observed $Q^2$ falls well above the permutation distribution, the model is valid.

11.4 Y-Randomisation Test

Related to the permutation test, Y-randomisation (also called response permutation) repeatedly randomises $\mathbf{y}$, fits PLS with the same $A$ components, and records $R^2_Y$ and $Q^2$:

• The permuted $Q^2$ values should be near zero or negative (since a randomly shuffled $\mathbf{y}$ has no true relationship with $\mathbf{X}$).
• The observed $Q^2$ should be far above the permuted distribution.
• If permuted models give high $Q^2$ values, the original result is likely spurious (due to chance correlations).

11.5 The Ratio $R^2_Y / Q^2$ as a Validity Check

A simple heuristic from chemometrics practice:

• If $Q^2 / R^2_Y > 0.5$: Model is not overfitting (good).
• If $Q^2 / R^2_Y < 0.5$ (or $Q^2 \ll R^2_Y$): Model may be overfitting — reduce number of components or collect more data.

12. Comparison with Related Methods

12.1 PLS vs. Principal Component Regression (PCR)

Both PLS and PCR reduce dimensionality before regression. The key difference:

💡 PLS is generally preferred over PCR for prediction tasks because it ensures the extracted components are relevant for predicting $\mathbf{Y}$. PCR may extract components that explain much $\mathbf{X}$ variance but are useless for prediction.

12.2 PLS vs. Ridge Regression

12.3 PLS vs. Lasso

12.4 PLS vs. OLS Multiple Regression

12.5 When to Choose Which Method

13. Using the PLS Regression Component

The PLS Regression component in the DataStatPro application provides a full end-to-end workflow for fitting, validating, and interpreting PLS models.

Step-by-Step Guide

Step 1 — Select Dataset Choose the dataset from the "Dataset" dropdown. The dataset should have at least one response variable and two or more predictor variables.

Step 2 — Select Analysis Type Choose the PLS analysis type:

• PLS1 (single continuous response)
• PLS2 (multiple continuous responses)
• PLS-DA (discriminant analysis; categorical response)

Step 3 — Select Predictor Variables (X) Select one or more predictor variables from the "Predictor Variables (X)" dropdown. These should be continuous or ordinal numeric variables.

💡 You can select all available numeric predictors and rely on VIP scores for post-hoc variable selection.

Step 4 — Select Response Variable(s) (Y)

• For PLS1: Select a single continuous response from the "Response Variable (Y)" dropdown.
• For PLS2: Select two or more continuous response variables.
• For PLS-DA: Select a categorical variable. You will be prompted to confirm the class coding.

Step 5 — Configure Preprocessing Select the preprocessing method for $\mathbf{X}$ (and optionally $\mathbf{Y}$):

• Mean Centring Only (recommended when variables are in the same units)
• Autoscaling (UV) (recommended default)
• Pareto Scaling
• No Preprocessing (not recommended unless data are already preprocessed)

Step 6 — Configure Number of Components Choose how to determine the number of components $A$:

• Automatic (Cross-Validation): The application selects $A$ by minimising RMSECV.
• Manual: Specify $A$ directly (e.g., based on domain knowledge or prior analysis).
• Set the maximum number of components to evaluate (default: $\min(n-1, p, 10)$).

Step 7 — Configure Cross-Validation Select the cross-validation scheme:

• $k$-Fold CV: Specify $k$ (default: 10).
• Leave-One-Out CV (LOOCV): Use for small datasets.
• Monte Carlo CV: Specify the number of iterations and train/test split ratio.

Step 8 — Select Display Options Choose which outputs to display:

• ✅ Score Plot ($t_1$ vs. $t_2$)
• ✅ Loading Plot ($p_1$ vs. $p_2$)
• ✅ Biplot ($w^*$-$q$ plot)
• ✅ $R^2_Y$ and $Q^2$ vs. Number of Components Plot
• ✅ RMSEC and RMSECV vs. Number of Components Plot
• ✅ Predicted vs. Observed Plot
• ✅ VIP Scores Plot
• ✅ Regression Coefficients Plot
• ✅ Residuals Plot
• ✅ Leverage-Residual (Williams) Plot
• ✅ $T^2$ and SPE (DModX) Control Charts
• ✅ Component Summary Table ($R^2_X$, $R^2_Y$, $Q^2$ per component)
• ✅ Permutation Plot

Step 9 — Run the Analysis Click "Run PLS Regression". The application will:

1. Apply the selected preprocessing to $\mathbf{X}$ and $\mathbf{Y}$.
2. Run cross-validation to determine the optimal number of components (if automatic).
3. Fit the final PLS model with the selected $A$ components using NIPALS or SIMPLS.
4. Compute scores ($\mathbf{T}$, $\mathbf{U}$), loadings ($\mathbf{P}$, $\mathbf{Q}$), weights ($\mathbf{W}$, $\mathbf{W}^*$).
5. Compute VIP scores and regression coefficients.
6. Compute model fit statistics ($R^2_X$, $R^2_Y$, $Q^2$, RMSEC, RMSECV).
7. Compute leverage, residuals, $T^2$, and SPE for all observations.
8. Generate all selected visualisations and tables.
9. Run permutation tests (if selected).

14. Computational and Formula Details

14.1 The NIPALS Algorithm for PLS1

The Non-linear Iterative Partial Least Squares (NIPALS) algorithm is the classical iterative procedure for PLS decomposition. For PLS1 (single response):

Inputs: Mean-centred (and scaled) $\mathbf{X}_1 = \tilde{\mathbf{X}}$ ($n \times p$) and $\mathbf{y}_1 = \tilde{\mathbf{y}}$ ($n \times 1$).

For each component $a = 1, 2, \dots, A$:

1. Initialise: $\mathbf{u}_a = \mathbf{y}_a$ (or any non-zero starting vector).

2. Compute X-weight vector: $\mathbf{w}_a = \frac{\mathbf{X}_a^T \mathbf{u}_a}{\|\mathbf{X}_a^T \mathbf{u}_a\|}$ (Normalise to unit length.)

3. Compute X-score vector: $\mathbf{t}_a = \mathbf{X}_a \mathbf{w}_a$

4. Compute Y-weight (scalar for PLS1): $c_a = \frac{\mathbf{y}_a^T \mathbf{t}_a}{\mathbf{t}_a^T \mathbf{t}_a}$

5. Compute Y-score vector: $\mathbf{u}_a = \mathbf{y}_a \cdot c_a / \|c_a\|^2 = \frac{\mathbf{y}_a^T \mathbf{t}_a}{\mathbf{t}_a^T \mathbf{t}_a} \mathbf{y}_a \cdot \frac{\mathbf{t}_a^T \mathbf{t}_a}{\mathbf{y}_a^T \mathbf{t}_a}$

For PLS1, since $q = 1$: $\mathbf{u}_a = \mathbf{y}_a$ (no iteration needed — convergence is immediate).

6. Check convergence (PLS2 only — for PLS1, skip to step 7).

7. Compute X-loadings: $\mathbf{p}_a = \frac{\mathbf{X}_a^T \mathbf{t}_a}{\mathbf{t}_a^T \mathbf{t}_a}$

8. Compute inner relation coefficient: $b_a = \frac{\mathbf{u}_a^T \mathbf{t}_a}{\mathbf{t}_a^T \mathbf{t}_a}$

9. Deflate $\mathbf{X}$: $\mathbf{X}_{a+1} = \mathbf{X}_a - \mathbf{t}_a \mathbf{p}_a^T$

10. Deflate $\mathbf{y}$: $\mathbf{y}_{a+1} = \mathbf{y}_a - b_a \mathbf{t}_a$

After all $A$ components, compute $\mathbf{W}^*$:

$\mathbf{W}^* = \mathbf{W}(\mathbf{P}^T \mathbf{W})^{-1}$

Where $\mathbf{W} = [\mathbf{w}_1, \mathbf{w}_2, \dots, \mathbf{w}_A]$ and $\mathbf{P} = [\mathbf{p}_1, \mathbf{p}_2, \dots, \mathbf{p}_A]$.

Compute PLS regression coefficients:

$\hat{\boldsymbol{\beta}}_{PLS} = \mathbf{W}^* \mathbf{b}$

Where $\mathbf{b} = [b_1, b_2, \dots, b_A]^T$ (inner relation coefficients, equal to Y-loadings $\mathbf{q}$ in PLS1).

14.2 The NIPALS Algorithm for PLS2

For PLS2 ($q > 1$), an inner iteration is needed at each component to converge to the dominant covariance direction:

1. Initialise: $\mathbf{u}_a =$ first column of $\mathbf{Y}_a$ (or any non-zero column).

2. Outer iteration (repeat until convergence of $\mathbf{t}_a$):

   a. Compute X-weight: $\mathbf{w}_a = \mathbf{X}_a^T \mathbf{u}_a / \|\mathbf{X}_a^T \mathbf{u}_a\|$

   b. Compute X-score: $\mathbf{t}_a = \mathbf{X}_a \mathbf{w}_a$

   c. Compute Y-weight: $\mathbf{c}_a = \mathbf{Y}_a^T \mathbf{t}_a / \|\mathbf{Y}_a^T \mathbf{t}_a\|$

   d. Compute Y-score: $\mathbf{u}_a = \mathbf{Y}_a \mathbf{c}_a$

   e. Check: if $\|\mathbf{t}_a^{(new)} - \mathbf{t}_a^{(old)}\| / \|\mathbf{t}_a^{(old)}\| < \epsilon$ (e.g., $\epsilon = 10^{-10}$), converged.

3. Compute loadings: $\mathbf{p}_a = \mathbf{X}_a^T \mathbf{t}_a / (\mathbf{t}_a^T \mathbf{t}_a)$

4. Compute inner coefficient: $b_a = \mathbf{u}_a^T \mathbf{t}_a / (\mathbf{t}_a^T \mathbf{t}_a)$

5. Deflate both matrices: $\mathbf{X}_{a+1} = \mathbf{X}_a - \mathbf{t}_a \mathbf{p}_a^T$ $\mathbf{Y}_{a+1} = \mathbf{Y}_a - b_a \mathbf{t}_a \mathbf{c}_a^T$