Centering and Scaling

The idea of standardizing measurements to make them comparable goes back to the beginnings of modern statistics. In 1893 and 1894, Karl Pearson introduced the term standard deviation, building on Francis Galton's work on regression and on how far individuals deviate from the average. Expressing a value as "how many standard deviations away from the mean" turns any variable, whatever its units, into a dimensionless number. Today we call that a z-score, and it lets us compare heights in centimetres with weights in kilograms, or temperatures with pressures.

When chemometrics emerged in the 1970s and 1980s, centering and scaling took on a new urgency. The workhorse methods of the field, Principal Component Analysis (PCA) and Partial Least Squares (PLS), are built on variances and covariances [1],[2]. Without centering, the first principal component of a set of spectra mostly describes the mean spectrum instead of the differences between samples, which are usually the whole point of the analysis. And as laboratories started to combine measurements of very different kinds (spectra, temperatures, pH values, concentrations) in one data matrix, scaling became unavoidable.

Over the decades, several scaling approaches were developed for different kinds of data. Autoscaling (unit variance scaling) became the standard for mixed-unit data. Pareto scaling, a softer alternative associated with the Umetrics school of multivariate analysis [3], became especially popular in metabolomics [4]. Range scaling found its niche where variables have natural operating ranges. The choice is not a minor technical detail: it decides which patterns a multivariate model is able to see [5].

Before applying PCA, PLS or almost any multivariate method, ask a simple question: are my variables on comparable scales?

In analytical chemistry, the answer is often no. A single data matrix might contain:

• NIR absorbances between 0.01 and 2.5,
• temperatures between 20 and 180 °C,
• moisture contents between 0.1 and 15 %,
• Raman intensities between 100 and 50 000 counts.

PCA looks for the directions in which the data vary the most. If Raman intensity varies over tens of thousands of counts while moisture varies over a few percent, the first principal component will describe Raman intensity almost exclusively. Not because it is more important chemically, but simply because its numbers are bigger.

Centering and scaling put the variables on a common footing, so that the model reflects the chemistry rather than the units someone happened to choose.

Mean centering is the most basic and most universally applied preprocessing step. For each variable (column), subtract its mean across all samples:

where

is the mean of variable $j$. In matrix notation:

where $1$ is a column of ones and $\bar{\mathbf{x}}$ is the vector of column means. The product $1{\bar{\mathbf{x}}}^T$ is simply a matrix in which every row is the mean.

After centering, every variable has mean zero. Geometrically, the whole cloud of samples is moved so that its centre sits at the origin. It sounds trivial, but it changes what PCA finds.

PCA, in its basic form, looks for lines through the origin that pass as close as possible to the data. If the data cloud is far away from the origin, the best such line is the one that points from the origin towards the cloud. That line describes where the data are (the mean), not how they vary.

Without centering, PC1 still "explains" almost all of the sum of squares, but all it tells you is the mean absorbance level: it captures only a few percent of the variation between samples. After centering, PC1 lines up with that variation and captures almost all of it. The demo shows both numbers, so you can see that a high "percentage explained" on uncentred data says nothing about how well the model describes the differences between samples.

In spectroscopy, this is exactly what happens with raw spectra. Without centering, the first component is essentially the mean spectrum and the differences between samples are pushed into later components. After centering, the first component describes the most important source of variation between samples, the second the next one, and so on [1],[2].

Autoscaling combines mean centering with division by the standard deviation. It is the z-score applied to every column:

with

the standard deviation of variable $j$. After autoscaling, every variable has mean zero and variance one, so each one starts with the same weight in the analysis, whatever its units.

Autoscaling is usually the right choice when:

• variables have different units, for example temperature, pressure and pH (when a whole spectrum is combined with a few such variables, add block scaling, see the decision guide below),
• variables have very different magnitudes, one between 0 and 1 and another between 0 and 10 000,
• all variables are potentially important and none should dominate just because of its scale.

Autoscaling has a well-known pitfall: it inflates variables that hardly vary. If a variable is almost constant, its standard deviation is tiny, and dividing by a tiny number magnifies whatever noise it contains. A variable that was practically irrelevant becomes as important as any other.

In spectroscopy this is a real problem. Flat spectral regions with no chemical information still contain noise. After autoscaling, those noisy regions get the same weight as the information-rich bands.

Pareto scaling is a compromise between no scaling and autoscaling. Instead of dividing by the standard deviation, divide by its square root:

Large variables are shrunk, but not all the way. The intuition is easy to see with numbers. If variable A has a standard deviation 100 times larger than variable B:

• no scaling: A has 100 times the spread of B, and 10 000 times its variance, which is what PCA sees,
• Pareto scaling: 10 times the spread ($\sqrt{100}=10$), 100 times the variance,
• autoscaling: equal spread and equal variance.

Pareto scaling became particularly popular in metabolomics, helped by the systematic comparison of scaling methods published by van den Berg and co-workers in 2006 [4]. In metabolomic data, large signals often belong to abundant, well-measured metabolites, while very small signals sit close to the noise. Autoscaling would give a peak at the noise floor the same weight as a major metabolite. Pareto scaling keeps large signals more influential, but stops them from completely drowning out the rest.

Range scaling divides each centred variable by its observed range:

It is related to min-max normalization (which maps each variable to the interval from 0 to 1), but keeps the variable centred at zero [4].

Range scaling is most natural in process monitoring, where variables have well-defined operating windows. A reactor that runs between 150 and 200 °C has a meaningful span of 50 °C; a pressure that runs between 1 and 5 bar has a span of 4 bar. Dividing by the span expresses each variable relative to its practical range.

Its weakness is that the range depends on only two samples: the largest and the smallest. A single outlier can stretch the range of a variable and shrink all its other values. Autoscaling is less sensitive, because the standard deviation uses every sample, but it is not robust either: a gross outlier inflates it too. Robust alternatives divide by the median absolute deviation (MAD) or the interquartile range (see outlier detection).

A few more approaches appear in the literature, mainly in metabolomics [4]:

• VAST scaling (variable stability scaling) autoscales each variable and then multiplies it by the inverse of its coefficient of variation, ${\bar{x}}_j/s_j$. Variables that are stable relative to their mean get more weight [6].
• Level scaling divides each centred variable by its mean, so the values become relative deviations from the average. Useful when relative changes matter more than absolute ones.
• Log transformation is not a scaling method but is often combined with one. It compresses data that span several orders of magnitude and turns multiplicative effects into additive ones. It requires positive values.
• Power transformations such as Box-Cox [8] form a family, indexed by a parameter $\lambda$, that includes the logarithm ($\lambda =0$) and, up to a linear rescaling, the square root ($\lambda =1/2$). They are used to stabilise variance or make distributions more symmetric.

The demo below brings everything together. Twenty-four samples from two production batches are described by six variables in different units. The difference between the batches is real, but it lives in variables with small numbers (pH, moisture, an absorbance). Conductivity, with numbers in the hundreds, varies a lot for reasons unrelated to the batch. One last variable is an absorbance in a flat spectral region: almost constant, pure noise.

Switch between the preprocessing options and watch three things: the common-axis strip plot at the top, the PCA scores, and the PC1 loadings (the weight of each variable in PC1). The "batch separation" number is the distance between the two batch centroids in the PC1-PC2 plane divided by the typical within-batch spread (the root-mean-square distance of each sample to its own batch centroid): below 1 the batches overlap, above 2 they are clearly apart.

What you should see:

• None. The strip plot shows conductivity far away from everything else. PCA without centering points PC1 at the mean, which is dominated by conductivity. The 100 % shown for PC1 is a share of the uncentred sum of squares, not of the variance. Batch separation is about 0.4.
• Mean centering. Better: the model now describes variation. But conductivity still owns PC1 (loading 1.00), and the batches overlap (separation still about 0.4).
• Autoscaling. Every variable has the same spread. The batch difference in pH, moisture and absorbance becomes the main direction of the data, and the batches separate (about 2.2). Notice also that the noisy flat-region absorbance now has as much variance as anything else: that is the price of autoscaling. Here it does not line up with the batch difference (its PC1 loading is about zero), so it adds noise to PC2 and later components rather than to PC1. In other datasets it can blur the main components too.
• Pareto. Conductivity is tamed (its values now span tens instead of hundreds), but it still dominates PC1, so in this dataset the batches stay mixed (about 0.4). Pareto is a compromise, and here the compromise is not enough.
• Range scaling. Close to autoscaling in this example (separation about 2.4).

The lesson is not "always autoscale". It is that preprocessing decides what the model can see, and that the right choice depends on what the variables mean.

There is no universal answer, but these guidelines cover most practical situations [1],[5].

The example builds three variables in different units: an absorbance, a temperature in °C and a moisture content in %. Absorbance and moisture are correlated; temperature varies a lot, but independently. With centering only, temperature owns PC1 (over 99 % of the variance). Autoscaling lets the moisture-absorbance correlation appear as PC1 (about 64 %). Pareto scaling sits in between, and here temperature still takes PC1 (about 94 %).

import numpy as np
import matplotlib.pyplot as plt

def fit_scaler(X, method="auto"):
    """Learn centering and scaling parameters from the training data."""
    mean = X.mean(axis=0)
    if method == "center":
        scale = np.ones(X.shape[1])
    elif method == "auto":
        scale = X.std(axis=0, ddof=1)
    elif method == "pareto":
        scale = np.sqrt(X.std(axis=0, ddof=1))
    elif method == "range":
        scale = X.max(axis=0) - X.min(axis=0)
    else:
        raise ValueError(method)
    scale[scale == 0] = 1.0          # constant variables: avoid division by zero
    return mean, scale

def apply_scaler(X, mean, scale):
    return (X - mean) / scale

def pca_scores(X, n_components=2):
    U, S, Vt = np.linalg.svd(X, full_matrices=False)
    explained = 100 * S**2 / np.sum(S**2)
    return U[:, :n_components] * S[:n_components], explained

# Mixed-unit example: absorbance, temperature (°C) and moisture (%)
rng = np.random.default_rng(42)
n = 50
moisture = 3 + 1.2 * rng.standard_normal(n)
absorbance = 0.3 + 0.02 * (moisture - 3) + 0.01 * rng.standard_normal(n)
temperature = 175 + 15 * rng.standard_normal(n)
X = np.column_stack([absorbance, temperature, moisture])

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
for ax, method in zip(axes, ["center", "auto", "pareto"]):
    mean, scale = fit_scaler(X, method)
    scores, expl = pca_scores(apply_scaler(X, mean, scale))
    ax.scatter(scores[:, 0], scores[:, 1], alpha=0.7)
    ax.set_title(method)
    ax.set_xlabel(f"PC1 ({expl[0]:.1f} %)")
    ax.set_ylabel(f"PC2 ({expl[1]:.1f} %)")
plt.tight_layout()
plt.show()

% Mixed-unit example: absorbance, temperature (°C) and moisture (%)
rng(42);
n = 50;
moisture = 3 + 1.2 * randn(n, 1);
absorbance = 0.3 + 0.02 * (moisture - 3) + 0.01 * randn(n, 1);
temperature = 175 + 15 * randn(n, 1);
X = [absorbance, temperature, moisture];

methods = {'center', 'auto', 'pareto'};
for k = 1:3
    [mu, sc] = fit_scaler(X, methods{k});
    Xs = (X - mu) ./ sc;
    [U, S, ~] = svd(Xs, 'econ');
    scores = U(:, 1:2) * S(1:2, 1:2);
    expl = 100 * diag(S).^2 / sum(diag(S).^2);
    subplot(1, 3, k);
    scatter(scores(:, 1), scores(:, 2), 'filled');
    title(methods{k});
    xlabel(sprintf('PC1 (%.1f %%)', expl(1)));
    ylabel(sprintf('PC2 (%.1f %%)', expl(2)));
end

function [mu, sc] = fit_scaler(X, method)
    % Learn centering and scaling parameters from the training data
    mu = mean(X, 1);
    switch method
        case 'center', sc = ones(1, size(X, 2));
        case 'auto',   sc = std(X, 0, 1);
        case 'pareto', sc = sqrt(std(X, 0, 1));
        case 'range',  sc = max(X, [], 1) - min(X, [], 1);
    end
    sc(sc == 0) = 1;
end

fit_scaler <- function(X, method = "auto") {
  # Learn centering and scaling parameters from the training data
  mu <- colMeans(X)
  sc <- switch(method,
    center = rep(1, ncol(X)),
    auto   = apply(X, 2, sd),
    pareto = sqrt(apply(X, 2, sd)),
    range  = apply(X, 2, function(v) diff(range(v))))
  sc[sc == 0] <- 1
  list(mean = mu, scale = sc)
}

apply_scaler <- function(X, s) sweep(sweep(X, 2, s$mean), 2, s$scale, "/")

# Mixed-unit example: absorbance, temperature (°C) and moisture (%)
set.seed(42)
n <- 50
moisture <- 3 + 1.2 * rnorm(n)
absorbance <- 0.3 + 0.02 * (moisture - 3) + 0.01 * rnorm(n)
temperature <- 175 + 15 * rnorm(n)
X <- cbind(absorbance, temperature, moisture)

par(mfrow = c(1, 3))
for (m in c("center", "auto", "pareto")) {
  Xs <- apply_scaler(X, fit_scaler(X, m))
  p <- prcomp(Xs, center = FALSE, scale. = FALSE)
  expl <- 100 * p$sdev^2 / sum(p$sdev^2)
  plot(p$x[, 1], p$x[, 2], main = m, pch = 19,
       xlab = sprintf("PC1 (%.1f %%)", expl[1]), ylab = sprintf("PC2 (%.1f %%)", expl[2]))
}

When you build a model, the means and standard deviations (or ranges) are computed from the training data and become part of the model. New samples must be preprocessed with those same numbers, never with numbers recomputed from the new data.

# schematic: fit_pls stands for any PLS implementation
# Calibration
mean, scale = fit_scaler(X_train, "auto")
model = fit_pls(apply_scaler(X_train, mean, scale), y_train)

# Prediction: reuse the calibration mean and scale
y_pred = model.predict(apply_scaler(X_new, mean, scale))   # NOT fit_scaler(X_new)

If you recomputed the means and standard deviations from each new batch, every batch would be shifted and stretched differently, and the model coefficients would no longer mean what they meant during calibration. The same rule applies inside cross-validation: compute the preprocessing parameters on the training part of each split only.

In PLS the response $\mathbf{y}$ is usually centred too (and scaled when there are several responses in different units). Its calibration mean must then be added back to the predictions.

Always do

• Mean-centre before PCA, PLS or any covariance-based method.
• Store the preprocessing parameters (means, standard deviations, ranges) from the calibration set and reuse them for every new sample.
• Document your choices: which preprocessing, in which order, and why.

Watch out for

• Autoscaling noisy variables. Nearly constant variables get inflated. Consider removing them, or use Pareto scaling.
• Inconsistent preprocessing between calibration, validation and prediction.
• Autoscaling spectra by reflex. For data in one unit, mean centering alone usually keeps more of the chemical structure.

Centering moves the data cloud to the origin, so PCA and PLS describe how samples differ instead of where they are. It is almost always needed.

Centering and scaling come last in the preprocessing chain, after row-wise corrections such as SNV, MSC or derivatives.

Scaling decides how much each variable is allowed to contribute. Autoscaling gives every variable the same weight, which is essential for mixed units (with block scaling when one block has far more variables than the others) but dangerous for nearly constant, noisy variables. Pareto and range scaling are compromises with their own strengths.

Preprocessing parameters belong to the model: compute them on the training data, and reuse them, unchanged, for everything that comes after.