Principal Component Analysis (PCA)

If you hand a chemometrician a new data set, the first thing they will almost certainly do is run a principal component analysis. PCA takes a table with many correlated variables, such as a few hundred spectra measured at a thousand wavelengths, and redraws it on a handful of new axes chosen so that the most important differences between the samples become visible. It does not need labels, concentrations or any prior knowledge. It simply asks: in which directions do these samples differ the most?

The idea is more than a century old. In 1901, the English mathematician and biometrician Karl Pearson published "On Lines and Planes of Closest Fit to Systems of Points in Space" [1]. His question was geometric: given a cloud of points, which straight line (or plane) passes as close as possible to all of them at once? He measured closeness perpendicular to the line, not vertically as in ordinary regression, and found that the best line runs through the centre of the cloud along the direction in which the points are most spread out.

Three decades later, the American statistician and economist Harold Hotelling reached the same mathematics from a different starting point. Working with batteries of psychological test scores, he looked for a small number of linear combinations of the tests that captured as much of their variance as possible. In 1933 he published "Analysis of a Complex of Statistical Variables into Principal Components", which gave the method its name and the eigenvalue framework we still use [2]. Pearson thought in terms of closest-fitting subspaces; Hotelling thought in terms of maximum variance. As you will see in a moment, these are two descriptions of the same thing.

For chemistry, the bridge came through computation. In 1966, Herman Wold described an iterative algorithm, later known as NIPALS, that extracts principal components one at a time using simple regressions [3]. It needs very little memory and it tolerates missing values, which made it well suited to the computers and data of the time. (Wold first called it NILES, for nonlinear iterative least squares; the name NIPALS came later.) His son Svante Wold and colleagues such as Kim Esbensen and Paul Geladi then made PCA a cornerstone of chemical data analysis during the 1970s and 1980s. Their 1987 tutorial [4] remains one of the most cited papers in chemometrics. Today PCA sits at the start of almost every chemometric workflow.

A modern near-infrared (NIR) spectrometer can record absorbance at more than a thousand wavelengths. Suppose you measure 200 olive oil samples, each at 1050 wavelengths between 900 and 2500 nm. Your data form a $200\times 1050$ matrix: 200 rows (samples) and 1050 columns (variables). Mathematically, every sample is a point in a space with 1050 dimensions.

But do you really need 1050 numbers to describe how these oils differ? Almost certainly not. The absorbance at 1720 nm is nearly the same information as the absorbance at 1721 nm, because both come from the same molecular vibration. Whole regions move together: the C-H overtone bands near 1210 nm and 1730 nm rise and fall with the fatty acid chains, and the O-H band near 1940 nm follows the water content. The 1050 wavelengths are not 1050 independent pieces of information. They are redundant measurements of a much smaller number of underlying sources of variation: a few chemical constituents, plus physical effects such as light scattering.

PCA answers two questions about such data:

• How many independent directions of variation does the data set really contain?
• What do they look like, both in terms of the samples (which are similar, which are different) and in terms of the variables (which wavelengths are involved)?

If the answer to the first question is "about five", PCA can compress 1050 dimensions into five with very little loss. That makes the data visible (you can plot two or three dimensions, not 1050), separates systematic variation from noise, and provides a stable starting point for regression and classification.

Imagine the samples as a cloud of points. Through that cloud, some directions show a lot of spread and others very little. The directions with the most spread carry the most information about how the samples differ.

The first principal component (PC1) is the direction along which the samples have the largest variance. Mathematically it is a weighted sum of the original variables,

where $x_1,\dots ,x_J$ are the $J$ original variables (for example, wavelengths) and the weights $p_{11},\dots ,p_{J1}$ define the direction (the first index is the variable, the second the component). PCA chooses the weights so that the values $t_1$ have the largest possible variance across the samples, with the constraint that the weight vector has unit length, ${\sum }_j{p}_{j1}^2=1$. Without that constraint we could make the variance as large as we liked just by multiplying all the weights by ten.

The second principal component (PC2) is the direction of largest variance that is perpendicular to PC1. It collects the largest part of what PC1 missed. The third is perpendicular to both, and so on. The result is a new set of axes, tailored to the data, in which:

• The first few axes capture most of the systematic variation.
• Later axes capture less and less, and eventually mostly noise.
• All axes are perpendicular, and the sample coordinates on different axes are uncorrelated.

The panel below lets you search for PC1 by hand. The points are a mean-centred cloud measured on two variables. Rotate the blue axis and watch how much of the variance it captures.

This little experiment shows why Pearson and Hotelling agree. For every sample, the squared distance from the centre splits, by Pythagoras, into a part along the axis and a part perpendicular to it. The total is fixed by the data. So making the captured part as large as possible is exactly the same as making the perpendicular residuals as small as possible.

Before computing PCA you must decide how to preprocess the data, and this choice determines what kind of variation PCA will find. Here is the short version; centring and scaling are covered in detail in their own article.

Mean centring subtracts the average of each variable (each column) from every sample:

Without centring, PC1 points roughly towards the average spectrum. It then "explains" a huge share of the total sum of squares, but what it describes is the fact that all samples absorb light, not how they differ. Centring moves the origin to the middle of the cloud so that the components describe differences between samples. For exploratory PCA, it is almost always done.

After centring there is a second choice:

• Centring only (PCA on the covariance matrix). Variables with larger variance have more influence. For spectra, where all variables share the same unit, this is usually what you want: strong, informative bands naturally weigh more than flat regions.
• Autoscaling (PCA on the correlation matrix). Each variable is also divided by its standard deviation, so every variable has unit variance and the same say. This is essential when the variables have different units (temperature, pH and absorbance in the same table), but for spectra it can backfire, because noisy regions with no chemical information get promoted to the same importance as the real bands.
• Pareto scaling divides each centred variable by the square root of its standard deviation. It is a compromise between the two, popular in metabolomics: large variables still weigh more, but not as overwhelmingly as without scaling.

The next panel runs a real PCA on 36 simulated NIR spectra of three products. Switch between the preprocessing options and see how the scores, loadings and scree plot change.

Now the maths behind the pictures [5],[6],[7]. Start from the centred (and possibly scaled) data matrix $\mathbf{X}$ with $I$ rows (samples) and $J$ columns (variables), as in most chemometrics texts. (Most other entries in this course, such as data reduction and the regression models, use $n$ and $p$ for the same two numbers.)

The covariance matrix summarises how every pair of variables varies together:

It is a symmetric $J\times J$ matrix. Its diagonal contains the variance of each variable, and the off-diagonal elements contain the covariances between pairs of variables. Two wavelengths that rise and fall together have a large positive covariance.

An eigenvector of a matrix is a direction that the matrix only stretches (or shrinks), without changing it; the stretching factor is the eigenvalue. The covariance matrix can be written entirely in terms of its eigenvectors and eigenvalues:

where:

• $\mathbf{P}$ is a $J\times J$ matrix whose columns are the eigenvectors. These are the loadings: the directions of the principal components. They are orthonormal, meaning perpendicular to each other and of unit length.
• $\mathbf{\Lambda }=\text{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_J)$ holds the eigenvalues, sorted so that ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge 0$.

The link with the rotation demo is direct: the eigenvalue ${\lambda }_a$ is exactly the variance of the samples projected onto the eigenvector ${\mathbf{p}}_a$. The eigenvector with the largest eigenvalue is the direction you found by turning the axis.

In practice, PCA is usually computed with the singular value decomposition (SVD) of the data matrix itself:

where $\mathbf{U}$ and $\mathbf{V}$ have orthonormal columns and $\mathbf{\Sigma }$ is diagonal with the singular values ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge 0$. The columns of $\mathbf{V}$ are the loadings, the scores are $\mathbf{U}\mathbf{\Sigma }$, and the eigenvalues are ${\lambda }_a={\sigma }_a^2/(I-1)$. Like the normal equations, forming ${\mathbf{X}}^T\mathbf{X}$ explicitly can lose numerical precision, so most software works with the SVD directly.

The SVD also tells us why PCA is such a good compressor: keeping the first $A$ components gives the best possible approximation of $\mathbf{X}$ by a matrix of rank $A$, in the least squares sense [8]. No other choice of $A$ directions leaves a smaller sum of squared residuals.

Each eigenvalue says how much variance its component captures. The proportion of variance explained by component $a$ is

and the cumulative proportion for the first $A$ components is

In spectroscopic data, the first few components often explain well over 90% of the variance, and the remaining hundreds of components share the rest, mostly noise. Be careful with these percentages, though: they depend heavily on preprocessing, as the spectra demo shows.

Keeping $A$ components, PCA writes the data as

where:

• $\mathbf{T}$ ($I\times A$) is the scores matrix: the new coordinates of each sample.
• $\mathbf{P}$ ($J\times A$) is the loadings matrix: the directions, one column per component.
• $\mathbf{E}$ ($I\times J$) is the residual matrix: what is left after $A$ components.

This compact equation is the heart of PCA interpretation. Each sample is rebuilt as a weighted sum of a few loading vectors, and the weights are its scores.

The scores are the projections of the samples onto the loadings:

Plotting the scores of one component against another gives a scores plot, the most important picture in exploratory chemometrics. In a scores plot:

• Samples close together have similar profiles.
• Clusters suggest groups with different compositions or origins.
• Gradients, smooth trends across the plot, suggest a quantity that varies continuously, such as a concentration.
• Isolated points may be outliers, worth a closer look.
• Ordered paths of samples measured over time reveal drift or ageing.

Colouring the points by known information (class, batch, date, concentration) is the quickest way to see which of these patterns the components are picking up.

The loadings tell you how each original variable contributes to each component. For spectra, a loading vector has one value per wavelength and can be plotted like a spectrum. Large positive or negative values mark the wavelengths that drive that component. If the PC1 loading has a large positive peak at 1730 nm, then samples with high PC1 scores absorb more at 1730 nm than the average sample, which you might relate to the C-H overtone of fatty acid chains.

Remember that a loading describes differences from the average (when the data are centred), and that its sign is arbitrary: flipping the signs of both the loading and the scores gives exactly the same model.

A biplot shows scores (as points) and loadings (as arrows) in the same figure, so you can see at once which samples are similar and which variables are responsible [9]. Biplots work well with a moderate number of variables, such as a dozen chemical measurements. With hundreds of wavelengths, the arrows turn into a hedgehog, and separate scores and loading plots are clearer.

Choosing the number of components is one of the main decisions in PCA. Too few, and real information is lost; too many, and you start modelling noise.

The scree plot shows the eigenvalues (or the percentage of explained variance) against the component number. Typically it drops steeply and then levels off. The "elbow" where it flattens suggests how many components to keep [10]. The name, coined by Raymond Cattell in 1966, comes from geology: scree is the rubble that piles up at the foot of a cliff, and the flat tail of the plot looks like it.

Another simple rule keeps as many components as needed to reach a threshold, such as 95% or 99% of the variance. It is easy to apply but arbitrary: the right threshold depends on how noisy the data are. In the compression panel of the data reduction article, a 99% threshold keeps one component too many and a 95% threshold one too few.

A more rigorous approach leaves out part of the data, builds the model on the rest, and checks how well the left-out part is reconstructed [11]. While a new component improves the reconstruction of unseen data, it captures structure; when it stops helping, it is fitting noise. Doing this correctly for PCA needs some care: the left-out values must not be used to compute their own scores. The simple "row-wise" scheme of many generic tools, which leaves out whole samples and then projects them onto the model, breaks that rule, and its error keeps falling as components are added. Element-wise schemes, which leave out individual values and predict them from the rest, avoid the problem [12]. Dedicated chemometrics software implements them; general-purpose libraries such as scikit-learn do not.

The broken stick rule compares each eigenvalue with what it would be if the total variance were shared at random. If a stick of unit length is broken at random into $J$ pieces, the expected length of the $a$-th longest piece is

and a component is kept if its share of the total variance exceeds $b_a$. The broken-stick model comes from ecology, where it was used for the relative abundance of species [13], and was later applied to eigenvalues [14]. In comparisons with other stopping rules it is one of the more reliable [15].

Two statistics tell you how each sample relates to the model. Hotelling's $T^2$ measures how far a sample lies from the centre within the model:

and the Q residual (also called the squared prediction error) measures how much of the sample the model does not describe [16]:

A high $T^2$ flags an extreme but familiar sample; a high $Q$ flags a sample with variation the model has never seen. For $T^2$ a common 95% limit, assuming normally distributed scores, is

where $F_{0.95}$ is the 95th percentile of the F distribution; the limits for $Q$ were derived by Jackson and Mudholkar [16]. Both statistics are explained, with limits and an interactive example, in the article on outlier detection.

The SVD computes all components at once. NIPALS (nonlinear iterative partial least squares) [3] extracts them one at a time, alternating between two small regressions. It remains useful because you can stop after the components you need, and because it can skip missing values.

Each round is a step of the classical power method, so NIPALS converges to the eigenvector with the largest eigenvalue of what remains. Convergence is fast when consecutive eigenvalues are well separated and slow when they are nearly equal.

Missing values are handled by leaving them out of the sums in steps 2 and 3, both in the numerator and in the normalising denominator. In step 3 that denominator is then no longer 1, because the observed part of $\mathbf{p}$ is shorter than the whole vector:

where the sums run only over the values that were actually measured. The article on missing data discusses this and other strategies.

The most common use of PCA in chemistry is the first look at a set of spectra or chromatographic profiles. The workflow is almost always the same: preprocess (scatter correction, smoothing, perhaps derivatives), centre, compute PCA, and study the scores and loadings. Typical questions are food authenticity (origin of olive oils or wines, adulteration of honey), pharmaceutical quality (identity of raw materials, consistency between batches), environmental monitoring (sources of pollution) and forensics (comparison of inks, fibres or paints). If the classes separate in the scores plot, PCA has found discriminating information. If they do not, the differences may be small compared to other sources of variation, and a supervised method may be needed.

PCA is also the basis of multivariate statistical process control (MSPC). A PCA model is built from data collected while a process runs normally, and every new measurement is projected onto it. Its $T^2$ and $Q$ values are followed in control charts. When one of them crosses its limit, contribution plots show which original variables are responsible, helping operators find the cause. It builds on control limits for the residuals derived by J. Edward Jackson and Govind Mudholkar in 1979 [16]. The approach was developed in the 1990s by John MacGregor, Theodora Kourti and colleagues at McMaster University [17], and in the process chemometrics community by Barry Wise and Neal Gallagher [18]. It is now a standard tool in the process analytical technology (PAT) framework used in pharmaceutical manufacturing.

Because the first components carry the structure and the last ones mostly noise, rebuilding the data from $A$ components gives a cleaner version of it:

and storing ${\mathbf{T}}_A$ and ${\mathbf{P}}_A$ takes far less space than $\mathbf{X}$ when $A$ is small. The article on data reduction shows this trade-off interactively.

The examples simulate NIR-like spectra of three olive oil cultivars, mean-centre them, and produce a scores plot, a scree plot and a loadings plot. Within each cultivar the composition varies a little from sample to sample, so PCA has real chemical variation to find, not only noise. Finally they compute $T^2$ and $Q$ for a four-component model.

import numpy as np
from scipy.stats import f
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# Simulate NIR spectra of olive oil samples from 3 cultivars
rng = np.random.default_rng(42)
n_per_class = 50
n_wavelengths = 500
wavelengths = np.linspace(900, 2500, n_wavelengths)

def peak(center, width2):
    return np.exp(-((wavelengths - center)**2) / width2)

# Each cultivar: typical amplitudes and positions of three bands
cultivars = {
    'Cultivar A': ([0.8, 0.4, 0.30], [1200, 1730, 1940]),
    'Cultivar B': ([0.6, 0.6, 0.35], [1210, 1725, 1935]),
    'Cultivar C': ([0.7, 0.5, 0.25], [1195, 1735, 1945]),
}

blocks = []
for amps, centres in cultivars.values():
    bands = np.vstack([peak(centres[0], 8000), peak(centres[1], 3000), peak(centres[2], 2000)])
    # composition varies from sample to sample (8% around the typical amplitudes)
    weights = np.array(amps) * (1 + 0.08 * rng.standard_normal((n_per_class, 3)))
    blocks.append(weights @ bands + rng.normal(0, 0.005, (n_per_class, n_wavelengths)))
X = np.vstack(blocks)
labels = np.repeat(list(cultivars.keys()), n_per_class)

# Mean centring (standard for spectroscopic data).
# sklearn's PCA would centre anyway; we do it explicitly to make the step visible.
X_centered = X - X.mean(axis=0)

pca = PCA(n_components=10)
scores = pca.fit_transform(X_centered)
loadings = pca.components_.T          # (n_wavelengths, n_components)
explained = pca.explained_variance_ratio_ * 100

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))

# Scores plot
colors = {'Cultivar A': '#2563eb', 'Cultivar B': '#dc2626', 'Cultivar C': '#059669'}
for label, color in colors.items():
    mask = labels == label
    axes[0].scatter(scores[mask, 0], scores[mask, 1], c=color, label=label,
                    alpha=0.7, edgecolors='k', linewidths=0.4)
axes[0].set_xlabel(f'PC1 ({explained[0]:.1f}%)')
axes[0].set_ylabel(f'PC2 ({explained[1]:.1f}%)')
axes[0].set_title('Scores')
axes[0].legend()

# Scree plot
comps = np.arange(1, 11)
axes[1].bar(comps, explained, color='steelblue')
axes[1].plot(comps, np.cumsum(explained), 'o-', color='#d97706', label='Cumulative')
axes[1].set_xlabel('Principal component')
axes[1].set_ylabel('Explained variance (%)')
axes[1].set_title('Scree plot')
axes[1].legend()

# Loadings
axes[2].plot(wavelengths, loadings[:, 0], label='PC1')
axes[2].plot(wavelengths, loadings[:, 1], label='PC2')
axes[2].axhline(0, color='gray', linewidth=0.5)
axes[2].set_xlabel('Wavelength (nm)')
axes[2].set_ylabel('Loading')
axes[2].set_title('Loadings')
axes[2].legend()

plt.tight_layout()
plt.show()

# Hotelling's T² and Q residuals for a 4-component model (the elbow of the scree plot)
A = 4
I = X.shape[0]
T2 = np.sum(scores[:, :A]**2 / pca.explained_variance_[:A], axis=1)
T2_lim = A * (I - 1) / (I - A) * f.ppf(0.95, A, I - A)      # 95% limit
E = X_centered - scores[:, :A] @ pca.components_[:A, :]
Q = np.sum(E**2, axis=1)
print(f"Variance explained (%): {np.round(explained[:5], 1)}")
print(f"Cumulative variance ({A} PCs): {explained[:A].sum():.1f}%")
print(f"T² limit (95%): {T2_lim:.2f}; samples above it: {np.sum(T2 > T2_lim)} of {I}")
print(f"Q range: {Q.min():.4f} to {Q.max():.4f}")

% Simulate NIR spectra of olive oil samples from 3 cultivars
rng(42);
n_per_class = 50;
n_wavelengths = 500;
wavelengths = linspace(900, 2500, n_wavelengths);
peak = @(c, w2) exp(-((wavelengths - c).^2) / w2);

% Each cultivar: typical amplitudes and positions of three bands
amps    = [0.8 0.4 0.30; 0.6 0.6 0.35; 0.7 0.5 0.25];
centres = [1200 1730 1940; 1210 1725 1935; 1195 1735 1945];

X = [];
for k = 1:3
    bands = [peak(centres(k, 1), 8000); peak(centres(k, 2), 3000); peak(centres(k, 3), 2000)];
    % composition varies from sample to sample (8% around the typical amplitudes)
    weights = amps(k, :) .* (1 + 0.08 * randn(n_per_class, 3));
    X = [X; weights * bands + 0.005 * randn(n_per_class, n_wavelengths)];
end
labels = repelem((1:3)', n_per_class);

% pca() centres the data by default and uses the SVD
[coeff, score, latent, ~, explained] = pca(X);

figure('Position', [100 100 1500 420]);

subplot(1, 3, 1);
gscatter(score(:, 1), score(:, 2), labels, [0.15 0.39 0.92; 0.86 0.15 0.15; 0.02 0.59 0.40]);
xlabel(sprintf('PC1 (%.1f%%)', explained(1)));
ylabel(sprintf('PC2 (%.1f%%)', explained(2)));
legend('Cultivar A', 'Cultivar B', 'Cultivar C'); title('Scores'); grid on;

subplot(1, 3, 2);
bar(1:10, explained(1:10)); hold on;
plot(1:10, cumsum(explained(1:10)), 'o-', 'LineWidth', 1.5);
xlabel('Principal component'); ylabel('Explained variance (%)');
title('Scree plot'); grid on;

subplot(1, 3, 3);
plot(wavelengths, coeff(:, 1:2), 'LineWidth', 1.5); yline(0, 'Color', [0.5 0.5 0.5]);
xlabel('Wavelength (nm)'); ylabel('Loading');
legend('PC1', 'PC2'); title('Loadings'); grid on;

% Hotelling's T² and Q residuals for a 4-component model (the elbow of the scree plot)
A = 4;
I = size(X, 1);
T2 = sum(score(:, 1:A).^2 ./ latent(1:A)', 2);
T2_lim = A * (I - 1) / (I - A) * finv(0.95, A, I - A);     % 95% limit
E = (X - mean(X)) - score(:, 1:A) * coeff(:, 1:A)';
Q = sum(E.^2, 2);
fprintf('Variance explained (%%): %s\n', mat2str(round(explained(1:5)', 1)));
fprintf('Cumulative variance (%d PCs): %.1f%%\n', A, sum(explained(1:A)));
fprintf('T² limit (95%%): %.2f; samples above it: %d of %d\n', T2_lim, sum(T2 > T2_lim), I);
fprintf('Q range: %.4f to %.4f\n', min(Q), max(Q));

# Simulate NIR spectra of olive oil samples from 3 cultivars
set.seed(42)
n_per_class <- 50
n_wavelengths <- 500
wavelengths <- seq(900, 2500, length.out = n_wavelengths)
peak <- function(center, width2) exp(-((wavelengths - center)^2) / width2)

# Each cultivar: typical amplitudes and positions of three bands
amps    <- rbind(c(0.8, 0.4, 0.30), c(0.6, 0.6, 0.35), c(0.7, 0.5, 0.25))
centres <- rbind(c(1200, 1730, 1940), c(1210, 1725, 1935), c(1195, 1735, 1945))

X <- NULL
for (k in 1:3) {
  bands <- rbind(peak(centres[k, 1], 8000), peak(centres[k, 2], 3000), peak(centres[k, 3], 2000))
  # composition varies from sample to sample (8% around the typical amplitudes)
  weights <- matrix(amps[k, ], n_per_class, 3, byrow = TRUE) *
             (1 + 0.08 * matrix(rnorm(n_per_class * 3), n_per_class, 3))
  X <- rbind(X, weights %*% bands +
                matrix(rnorm(n_per_class * n_wavelengths, sd = 0.005), n_per_class))
}
labels <- rep(c("Cultivar A", "Cultivar B", "Cultivar C"), each = n_per_class)

# prcomp centres by default; no scaling for spectra
pca_result <- prcomp(X, center = TRUE, scale. = FALSE)
scores <- pca_result$x
loadings <- pca_result$rotation
explained <- 100 * pca_result$sdev^2 / sum(pca_result$sdev^2)

par(mfrow = c(1, 3), mar = c(4.5, 4.5, 3, 1))
cols <- c("Cultivar A" = "#2563eb", "Cultivar B" = "#dc2626", "Cultivar C" = "#059669")

plot(scores[, 1], scores[, 2], col = cols[labels], pch = 16,
     xlab = sprintf("PC1 (%.1f%%)", explained[1]),
     ylab = sprintf("PC2 (%.1f%%)", explained[2]), main = "Scores")
legend("topright", legend = names(cols), col = cols, pch = 16)

bar_x <- barplot(explained[1:10], names.arg = 1:10, col = "steelblue",
                 xlab = "Principal component", ylab = "Explained variance (%)",
                 main = "Scree plot", ylim = c(0, 100))
lines(bar_x, cumsum(explained[1:10]), type = "b", pch = 16, col = "#d97706")

matplot(wavelengths, loadings[, 1:2], type = "l", lty = 1, lwd = 1.5,
        col = c("#2563eb", "#dc2626"), xlab = "Wavelength (nm)",
        ylab = "Loading", main = "Loadings")
abline(h = 0, col = "gray")
legend("topright", legend = c("PC1", "PC2"), col = c("#2563eb", "#dc2626"), lwd = 1.5)

# Hotelling's T² and Q residuals for a 4-component model (the elbow of the scree plot)
A <- 4
I <- nrow(X)
T2 <- rowSums(sweep(scores[, 1:A]^2, 2, pca_result$sdev[1:A]^2, "/"))
T2_lim <- A * (I - 1) / (I - A) * qf(0.95, A, I - A)     # 95% limit
E <- scale(X, center = TRUE, scale = FALSE) - scores[, 1:A] %*% t(loadings[, 1:A])
Q <- rowSums(E^2)
cat(sprintf("Variance explained (%%): %s\n", paste(round(explained[1:5], 1), collapse = " ")))
cat(sprintf("Cumulative variance (%d PCs): %.1f%%\n", A, sum(explained[1:A])))
cat(sprintf("T² limit (95%%): %.2f; samples above it: %d of %d\n", T2_lim, sum(T2 > T2_lim), I))
cat(sprintf("Q range: %.4f to %.4f\n", min(Q), max(Q)))

In Python the first components explain about 76.5%, 13.2%, 7.4% and 1.1% of the variance, and then the scree plot falls flat: four components carry the chemistry, and the rest is noise (the other languages draw different random numbers, so their percentages differ slightly). With four components, $Q$ is about the same for every sample, close to $500\times 0.005^2=0.0125$, the sum of squares of the simulated noise over 500 wavelengths: what is left really is noise. Only a few samples exceed the 95% $T^2$ limit, fewer than the 5% (7 or 8 samples) expected by chance, so none stands out as an outlier.

PCA is the right first tool when you want to:

• Explore a new data set: see groups, trends and outliers before any modelling.
• Compress many correlated variables into a few uncorrelated scores.
• Reduce noise by rebuilding the data from the first components.
• Monitor a process against a model of normal behaviour with $T^2$ and $Q$.
• Prepare for regression or classification, for example in principal component regression.

It is not the right tool when:

• You need to predict a property. PCA ignores the response; the direction of largest variance is not necessarily the most predictive. Use PCR or, usually better, PLS.
• You need to separate known classes. PCA maximises variance, not class separation. If classes overlap in the scores plot, consider supervised methods such as PLS-DA (partial least squares discriminant analysis, outlined in the PLS article) or LDA (linear discriminant analysis), which look for the directions that best separate known classes.
• The structure is strongly nonlinear. PCA finds straight directions; data that lie on a curved surface may need nonlinear methods.
• You have very few samples. With a handful of samples the components can describe noise. Cross-validation becomes essential.

PCA rotates the coordinate system so that the first axis follows the largest variation in the samples, the second the largest remaining variation, and so on. The scores tell you where each sample lies on those axes; the loadings tell you which variables build each axis. Looking for the largest variance and looking for the closest-fitting line are the same thing, which is why PCA is both a visualisation tool and the best low-rank approximation of the data.

From here, several methods build on PCA:

• Factor analysis adds an explicit model of latent factors and noise, and uses rotation to make the factors easier to interpret.
• Clustering assigns samples to groups algorithmically, often working on PCA scores.
• Principal component regression regresses a response on the PCA scores.
• Partial least squares builds components that are guided by the response, which makes them better suited to prediction.