Principal Component Regression (PCR)

A spectrum may have a thousand wavelengths, but it rarely has a thousand independent things to say. If a set of samples varies in three constituents, a baseline and a little scatter, the spectra are, up to noise, combinations of perhaps five underlying patterns. The thousand numbers are highly redundant.

Principal component regression (PCR) takes that observation literally. First it compresses the spectra into a handful of new variables, the principal component scores, which summarize the main patterns of variation and are uncorrelated with each other. Then it runs an ordinary regression on those few scores. Compress first, regress second: the collinearity that ruined multiple linear regression simply disappears.

PCR combines two older tools. Principal component analysis goes back to Karl Pearson in 1901 [2] and Harold Hotelling in 1933 [3]; least squares regression is older still. The idea of regressing on principal components had been suggested in 1957 by Maurice Kendall [12] and by Hotelling himself [13]. In 1965, William F. Massy, then at Stanford, published the first systematic treatment, Principal Components Regression in Exploratory Statistical Research [1], proposing to regress the response on the leading principal components when there are too many correlated predictors to handle directly.

In chemistry, PCR became important with the rise of near-infrared spectroscopy in the 1980s. NIR spectra have hundreds or thousands of strongly correlated wavelengths and are often collected on relatively few reference samples. Principal components were already being used to analyse NIR spectra of cereals in the mid-1980s [14], and papers by Næs and Martens [4] and by Haaland and Thomas [5] in 1988 established PCR, together with the newer partial least squares (PLS) method, as one of the two pillars of multivariate calibration.

Statisticians had also pointed out its weak spot. Ian Jolliffe noted in 1982 that the components that explain the most variance in the predictors are not necessarily the ones that predict the response, and backed it with real datasets in which low-variance components turned out to be predictive [6]. That remark is the key to understanding both PCR and its rival, PLS.

With a spectral matrix $\mathbf{X}$ of 150 samples and 1000 wavelengths, least squares needs to invert the $1000\times 1000$ matrix ${\mathbf{X}}^T\mathbf{X}$. Its rank is at most 149 once the data are centred (at most 150 without centring), so it is singular: infinitely many coefficient vectors fit the calibration perfectly, and almost all of them are useless for prediction. Even with more samples than wavelengths, neighbouring wavelengths are so correlated that ${\mathbf{X}}^T\mathbf{X}$ is nearly singular and the coefficients are wildly unstable, as we saw in the collinearity demo.

There are three broad ways out:

1. Regularization: keep all variables but restrain the coefficients (ridge, lasso).
2. Variable selection: keep only a few variables (lasso, stepwise methods).
3. Latent variables: replace the many correlated variables by a few new ones that capture the essential structure (PCR, PLS).

PCR is the simplest member of the third family.

Principal component analysis (see the PCA entry for a full treatment) looks for the direction in variable space along which the samples spread the most. That direction is the first loading vector ${\mathbf{p}}_1$, and the position of each sample along it is its first score $t_{i1}$. The second component is the direction of largest remaining spread that is perpendicular to the first, and so on. Each component is a "pattern" (the loading, which for spectra looks like a spectrum) and a set of "amounts" (the scores, one number per sample). The scores of different components are uncorrelated, and the first few components usually describe almost all of the systematic variation in the data, leaving mostly noise for the later ones.

All of the ill-conditioning is absorbed in step 2, where PCA decides which few directions to keep. Step 3 never sees the troublesome directions at all.

Start from the singular value decomposition of the centred spectra, $\mathbf{X}=\mathbf{UD}{\mathbf{V}}^T$. The columns of $\mathbf{V}$ are the loadings, the singular values $d_a$ (on the diagonal of $\mathbf{D}$) measure how much the data spread along each loading, and the scores are $\mathbf{T}=\mathbf{XV}=\mathbf{UD}$. Equivalently, the loadings are the eigenvectors of ${\mathbf{X}}^T\mathbf{X}$ and the eigenvalues are $d_a^2$.

Keeping the first $A$ components, ${\mathbf{T}}_A=\mathbf{X}{\mathbf{P}}_A$, and since the scores are orthogonal, ${\mathbf{T}}_A^T{\mathbf{T}}_A=\text{diag}(d_1^2,\dots ,d_A^2)$. The regression of $\mathbf{y}$ on the scores and its translation back to wavelengths become:

Compare the last expression with least squares and ridge regression, written in the same SVD language in the ridge entry:

PCA builds its components from $\mathbf{X}$ alone. It has no idea what you want to predict. The first component is the direction of largest spectral variation, whatever its cause. If the largest variation in your spectra comes from water content, particle size or temperature, the first components will describe exactly that, and the analyte may only appear in a later component.

The demo below makes this concrete with two wavelengths. The cloud of samples is stretched along one direction by a large interferent. Use the slider to decide where the analyte information lives, and compare the first principal component with the first PLS direction, which (as we will see in the next entry) is built using the response.

When the analyte information lies along the big variation, PC1 is a great predictor. When it lies along the small variation, PC1 points the wrong way and its scores have almost no relation to the analyte: a one-component PCR model is useless, and you need the second component as well. With real spectra and several interferents, the analyte may only appear in the fourth or fifth component, and Hadi and Ling showed that, in principle, all the predictive information can even sit in the last component [7].

The practical consequences:

• PCR often needs more components than PLS to reach the same prediction error.
• Selecting components strictly in order of variance is a choice, not a law. Some practitioners select components by their correlation with $\mathbf{y}$ instead, at the price of a more complicated (and more easily overfitted) selection step.

The number of components $A$ is the key tuning parameter of PCR, and it controls a familiar trade-off:

• Too few components leave out spectral information that the model needs: the predictions are biased.
• Too many components bring in directions that describe noise or irrelevant variation: the model starts to fit the calibration samples rather than the chemistry, and it predicts new samples worse.

The standard tool is cross-validation:

The demo below runs this procedure on twenty simulated NIR spectra. Move the number of components and watch four things at once: the regression vector, the calibration and cross-validated errors, and how much of the variance in $\mathbf{X}$ and in $\mathbf{y}$ each component explains. The green point marks the suggested number of components: the fewest whose cross-validated error is within 5 % of the minimum (the parsimony rule described below). The grey curve is the cross-validated error of PLS, the method of the next entry; you can switch to it to compare.

With a small analyte signal, the first components of PCR explain a lot of spectral variance but almost none of the analyte (look at the green bars): they describe the interferents and the baseline. The cross-validated error stays high until the component that carries the analyte finally arrives. Meanwhile the calibration error keeps going down, which is why it can never be used to choose $A$.

Now make the analyte signal dominant. Once the analyte is the main source of spectral variation, the first components already carry it, and PCR catches up with PLS.

The examples simulate NIR-like spectra built from three underlying factors, choose the number of components by 10-fold cross-validation with a parsimony rule, and plot the result.

import numpy as np
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_predict, KFold
from sklearn.pipeline import make_pipeline
import matplotlib.pyplot as plt

# --- Simulate NIR calibration data: 150 samples, 500 wavelengths ---
rng = np.random.default_rng(42)
n_samples, n_wavelengths = 150, 500
wavelengths = np.linspace(900, 1700, n_wavelengths)

loading1 = np.exp(-((wavelengths - 1200) ** 2) / 5000)
loading2 = np.exp(-((wavelengths - 1400) ** 2) / 3000)
loading3 = np.exp(-((wavelengths - 1100) ** 2) / 8000)
f1 = rng.standard_normal(n_samples)
f2 = rng.standard_normal(n_samples)
f3 = 0.5 * rng.standard_normal(n_samples)

X = np.outer(f1, loading1) + np.outer(f2, loading2) + np.outer(f3, loading3)
X += 0.05 * rng.standard_normal((n_samples, n_wavelengths))
y = 3.0 * f1 + 1.5 * f2 + 0.3 * rng.standard_normal(n_samples)

# --- Choose the number of components by cross-validation ---
# PCA centres X internally and LinearRegression fits the intercept,
# so every fold is centred with its own training data.
max_components = 15
kf = KFold(n_splits=10, shuffle=True, random_state=42)
rmsecv = []
for a in range(1, max_components + 1):
    pcr = make_pipeline(PCA(n_components=a), LinearRegression())
    y_cv = cross_val_predict(pcr, X, y, cv=kf)
    rmsecv.append(np.sqrt(np.mean((y - y_cv) ** 2)))

# Parsimony: the fewest components whose RMSECV is within 2 % of the minimum
rmsecv = np.array(rmsecv)
best_a = int(np.argmax(rmsecv <= 1.02 * rmsecv.min())) + 1
print(f"Chosen number of components: {best_a}  (RMSECV = {rmsecv[best_a - 1]:.3f})")

plt.figure(figsize=(8, 4))
plt.plot(range(1, max_components + 1), rmsecv, "o-")
plt.axvline(best_a, color="red", linestyle="--")
plt.xlabel("Number of principal components")
plt.ylabel("RMSECV")
plt.title("PCR: choosing the number of components")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# --- Final model and its regression vector in wavelength space ---
pca = PCA(n_components=best_a).fit(X)
T = pca.transform(X)
reg = LinearRegression().fit(T, y)
beta_pcr = pca.components_.T @ reg.coef_          # one coefficient per wavelength
y_fit = (X - pca.mean_) @ beta_pcr + reg.intercept_

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].bar(range(1, 11), 100 * PCA(n_components=10).fit(X).explained_variance_ratio_)
axes[0].set_xlabel("Component"); axes[0].set_ylabel("Explained X variance (%)")
axes[1].plot(wavelengths, beta_pcr)
axes[1].set_xlabel("Wavelength (nm)"); axes[1].set_ylabel("Regression coefficient")
axes[2].scatter(y, y_fit, alpha=0.6)
axes[2].plot([y.min(), y.max()], [y.min(), y.max()], "r--")
axes[2].set_xlabel("Reference"); axes[2].set_ylabel("Fitted")
for ax in axes:
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

% Principal component regression
% Requires the Statistics and Machine Learning Toolbox
rng(42);
n_samples = 150;
n_wavelengths = 500;
wavelengths = linspace(900, 1700, n_wavelengths);

loading1 = exp(-((wavelengths - 1200).^2) / 5000);
loading2 = exp(-((wavelengths - 1400).^2) / 3000);
loading3 = exp(-((wavelengths - 1100).^2) / 8000);
f1 = randn(n_samples, 1);  f2 = randn(n_samples, 1);  f3 = 0.5 * randn(n_samples, 1);

X = f1 * loading1 + f2 * loading2 + f3 * loading3 + 0.05 * randn(n_samples, n_wavelengths);
y = 3.0 * f1 + 1.5 * f2 + 0.3 * randn(n_samples, 1);

% --- Cross-validation, centring inside every fold ---
max_comp = 15;
cvp = cvpartition(n_samples, 'KFold', 10);
press = zeros(max_comp, 1);
for f = 1:cvp.NumTestSets
    tr = training(cvp, f);  te = test(cvp, f);
    mu = mean(X(tr, :));  ybar = mean(y(tr));
    Xtr = X(tr, :) - mu;  Xte = X(te, :) - mu;
    [~, ~, V] = svd(Xtr, 'econ');
    for a = 1:max_comp
        P = V(:, 1:a);
        q = (Xtr * P) \ (y(tr) - ybar);
        y_hat = ybar + Xte * P * q;
        press(a) = press(a) + sum((y(te) - y_hat).^2);
    end
end
rmsecv = sqrt(press / n_samples);
% Parsimony: the fewest components whose RMSECV is within 2 % of the minimum
best_a = find(rmsecv <= 1.02 * min(rmsecv), 1);
fprintf('Chosen number of components: %d (RMSECV = %.3f)\n', best_a, rmsecv(best_a));

figure;
plot(1:max_comp, rmsecv, 'o-', 'LineWidth', 1.5);
xline(best_a, 'r--');
xlabel('Number of principal components'); ylabel('RMSECV');
title('PCR: choosing the number of components'); grid on;

% --- Final model ---
mu = mean(X);  ybar = mean(y);
Xc = X - mu;
[~, S, V] = svd(Xc, 'econ');
P = V(:, 1:best_a);
q = (Xc * P) \ (y - ybar);
beta_pcr = P * q;                                  % one coefficient per wavelength
y_fit = ybar + Xc * beta_pcr;

explained = 100 * diag(S).^2 / sum(diag(S).^2);
figure;
subplot(1, 3, 1); bar(explained(1:10)); xlabel('Component'); ylabel('Explained X variance (%)');
subplot(1, 3, 2); plot(wavelengths, beta_pcr); xlabel('Wavelength (nm)'); ylabel('Coefficient');
subplot(1, 3, 3); scatter(y, y_fit, 'filled'); hold on;
plot([min(y) max(y)], [min(y) max(y)], 'r--'); xlabel('Reference'); ylabel('Fitted');

library(pls)  # standard R package for PCR and PLS
set.seed(42)

n_samples <- 150
n_wavelengths <- 500
wavelengths <- seq(900, 1700, length.out = n_wavelengths)

loading1 <- exp(-((wavelengths - 1200)^2) / 5000)
loading2 <- exp(-((wavelengths - 1400)^2) / 3000)
loading3 <- exp(-((wavelengths - 1100)^2) / 8000)
f1 <- rnorm(n_samples); f2 <- rnorm(n_samples); f3 <- 0.5 * rnorm(n_samples)

X <- outer(f1, loading1) + outer(f2, loading2) + outer(f3, loading3) +
     matrix(rnorm(n_samples * n_wavelengths, 0, 0.05), n_samples, n_wavelengths)
y <- 3.0 * f1 + 1.5 * f2 + rnorm(n_samples, 0, 0.3)

dat <- data.frame(y = y)
dat$X <- I(X)

# --- PCR with 10-fold cross-validation (pcr() centres X; scale = FALSE keeps spectra unscaled) ---
pcr_model <- pcr(y ~ X, data = dat, ncomp = 15, validation = "CV", segments = 10, scale = FALSE)
summary(pcr_model)

rmsecv <- RMSEP(pcr_model, estimate = "CV")
plot(rmsecv, main = "PCR: choosing the number of components")
# Parsimony: the one-standard-error rule picks the fewest components whose
# RMSECV is within one standard error of the minimum
best_a <- selectNcomp(pcr_model, method = "onesigma", plot = TRUE)
cat("Chosen number of components:", best_a, "\n")

# --- Explained variance, regression vector and fitted values ---
par(mfrow = c(1, 3))
barplot(explvar(pcr_model)[1:10], names.arg = 1:10,
        xlab = "Component", ylab = "Explained X variance (%)")
plot(wavelengths, drop(coef(pcr_model, ncomp = best_a)), type = "l",
     xlab = "Wavelength (nm)", ylab = "Regression coefficient")
plot(y, drop(fitted(pcr_model)[, , best_a]), xlab = "Reference", ylab = "Fitted")
abline(0, 1, col = "red", lty = 2)

Both the Python and the MATLAB runs choose three components, one per underlying factor. From three components on, the RMSECV curve is flat (about 0.32 in Python and 0.29 in MATLAB), and its exact minimum falls at six components in Python and four in MATLAB, only because of random wobbles of less than 1 %. Taking the plain minimum would add components that buy nothing, which is exactly what the parsimony rule prevents. The R version uses the one-standard-error rule of the pls package, which follows the same idea.

Good applications:

• Calibration of major constituents by NIR or MIR, when the analyte is itself a dominant source of spectral variation: PCR then performs about as well as PLS.
• Exploratory calibration: the same scores and loadings you already inspect in PCA become the basis of the model, which makes the model easy to explain.
• Teaching and understanding multivariate calibration: PCR is the clearest way to see the compress-then-regress idea before meeting PLS.

Consider alternatives when:

• the analyte is a minor constituent whose signal is buried under larger sources of variation: PLS finds the relevant directions sooner;
• you want the fewest possible components, for robustness or for interpretation of each component;
• the first principal components clearly describe interferences rather than the analyte (check the scores coloured by the response).

Mean-centre X and y. Without centring, the first component mostly describes the mean spectrum instead of the variation between samples, wasting a component.

Do not autoscale spectra without a reason. Scaling each wavelength to unit variance gives noisy baseline regions the same weight as informative peaks. For spectra, mean-centring is usually the right choice; autoscaling is appropriate when variables have different units.

Look at the loadings. If the loadings of the retained components show features at chemically meaningful wavelengths, the model is using real chemistry. Loadings that look like noise suggest too many components.

Colour the score plots by y. A gradient of colour across PC1 and PC2 means the first components are predictive. If the gradient only appears in PC5 versus PC6, PLS deserves a try.

Preprocess first. Scatter corrections (such as SNV or MSC), smoothing and derivatives remove physical effects that would otherwise dominate the first components.

Validate independently. Cross-validation chooses $A$; the final performance should be checked on samples that played no part in building the model, including the choice of $A$. That error is the RMSEP.

PCR replaces the original correlated variables by a few uncorrelated principal component scores and regresses the response on them. In the SVD language it keeps the $A$ largest directions of the data completely and discards the rest: a hard version of the shrinkage that ridge regression applies smoothly. That makes it stable, transparent and easy to connect with exploratory PCA.

Its weakness is that the components are chosen without looking at the response. Large spectral variation is not necessarily relevant variation, so PCR sometimes needs many components before it captures the analyte. Partial least squares keeps the compress-then-regress idea but lets the response guide the compression, and that single change is why it became the workhorse of chemometric calibration. For more on the practice of PCR in spectroscopy, see Martens and Næs [10]; for PCA itself, Wold, Esbensen and Geladi [9] and Jolliffe [11].