Outlier Detection

Every analyst has met that one measurement. The replicate that sits far from the others, the spectrum that looks a little odd, the calibration standard that refuses to sit on the line. Is it a mistake, or is it the most interesting sample in the data set? Outlier detection is the set of tools that helps us ask that question in a fair, reproducible way.

The question is old. In 1863, the American mathematician and astronomer William Chauvenet included in his Manual of Spherical and Practical Astronomy a rule for rejecting suspicious observations [1], a simpler version of a criterion proposed by Benjamin Peirce in 1852 [21]. Astronomers repeated their measurements of star positions many times, and now and then one reading drifted far from the rest: a misread scale, a vibration of the telescope, a moment of distraction. Chauvenet's criterion said that an observation should be rejected if the probability of getting a deviation that large was smaller than $1/(2n)$ for a sample of size $n$. The rule is crude by modern standards, but the principle behind it is still ours: we need a quantitative rule to decide when a value is too far from the rest, because human judgement alone is inconsistent.

Almost a century later, Frank E. Grubbs, a statistician at the Aberdeen Proving Ground in Maryland, placed the problem on firmer ground. In 1950 and in a widely cited 1969 review, he developed formal hypothesis tests for one or more outliers in a sample drawn from a normal distribution [2],[3]. Variants of Grubbs' test are still used in quality control and in standard methods such as ASTM E178, and ISO 5725-2, the standard for interlaboratory precision studies, prescribes it together with Cochran's test for variances.

The next big step came in 1977, when R. Dennis Cook introduced what we now call Cook's distance: a measure of how much a regression changes when a single observation is deleted [4]. This shifted the question from "is this point unusual?" to "is this point pulling the model?". In chemometrics that second question is crucial. A single contaminated spectrum in a PCA or PLS model can rotate the components, distort the score plot and degrade the predictions for every other sample. The tools developed for this setting (Hotelling's $T^2$, Q residuals, leverage and their combinations) are now part of every serious chemometric workflow.

An outlier is an observation that does not follow the same pattern as the majority of the data. The definition is deliberately vague, because "not following the pattern" depends on which pattern you expect.

With one variable, an outlier is simply a value far from the centre. A pH of 14.3 in a set of readings around 7.0 needs no statistics to be noticed. With many variables, outliers can be much more subtle. A sample may have perfectly ordinary values for every single variable and still be unusual in its combination of values: for example, a high absorbance at one wavelength together with a low absorbance at another wavelength that normally rises with it. Checking each variable on its own will never find that sample. This is why multivariate data need multivariate diagnostics.

Two ideas are easy to confuse:

• An outlier is far from the rest of the data. It may or may not affect a model.
• An influential point is a point that, if removed, noticeably changes the fitted model.

A point can be influential without looking like an outlier: if it sits at the edge of the range of the predictors, it can quietly drag the model towards itself. And a point can be a clear outlier without being influential: if it sits in the middle of the predictor range, it may barely change the slope of a regression line.

Most chemometric methods (least squares regression, PCA, PLS) are built on squared deviations. Squaring gives large deviations a disproportionate weight: a residual ten times larger than usual contributes a hundred times more to the sum of squares. When all the data are well behaved, this is what makes least squares efficient. When one sample is wrong, it is exactly what lets that sample dominate the model. In PCA, a single extreme sample can:

• Rotate the principal components towards itself, distorting the model for everybody else.
• Use up a whole component just to describe its own peculiarity.
• Inflate the apparent explained variance and give a misleading picture of the data.
• Bias the regression coefficients of a PLS model, degrading predictions for normal samples.

Before moving to multivariate tools, it is worth knowing the classical univariate rules. They cannot detect multivariate outliers, but they are useful for screening individual variables, such as a reference concentration or a replicate measurement.

The simplest approach is to measure how many standard deviations each value lies from the mean:

where $\bar{x}$ is the sample mean and $s$ the sample standard deviation. A common rule flags $|z_i|>3$, because under a normal distribution only about 0.3% of values should be that far out.

With small samples the rule cannot even fire: no value can be more than $(n-1)/\sqrt{n}$ standard deviations from the mean, which is below 3 for $n\le 10$. With five replicates, for example, the largest possible $|z|$ is about 1.79, however wild the suspect value is.

The weakness is that the mean and the standard deviation are themselves pulled by the outlier. An extreme value inflates $s$, which shrinks every z-score, including its own. The outlier partly hides itself. This is called the masking effect. Its twin is swamping: a distorted mean and standard deviation can make perfectly good values look suspicious.

A more robust alternative uses the quartiles. The interquartile range is

where $Q_1$ and $Q_3$ are the first and third quartiles (the values below which 25% and 75% of the data lie). A value is flagged if it is below $Q_1-1.5\text{IQR}$ or above $Q_3+1.5\text{IQR}$. This is the rule behind John Tukey's box plot [11]. For normal data the fences sit at roughly $\pm 2.7\sigma$, and using 3 instead of 1.5 identifies "far out" values.

Quartiles are barely affected by a few extreme values, so the rule does not suffer from masking in the way the z-score does.

Grubbs' test asks whether the most extreme value in a sample is consistent with all values coming from the same normal population [2],[3]. The test statistic is the largest absolute z-score:

The hypothesis that there is no outlier is rejected at significance level $\alpha$ if

where $t_{\alpha /(2n),n-2}$ is the upper critical value of Student's t distribution with $n-2$ degrees of freedom.

Analytical chemists often meet Dixon's Q test, designed for very small sets of replicates (roughly 3 to 10) [16]. Sort the values; the statistic is the gap between the suspect value and its nearest neighbour, divided by the range of all the values:

and the suspect is rejected when $Q$ exceeds a tabulated critical value for the given $n$. Dixon's $Q$ has nothing to do with the Q residual discussed below; the shared letter is a coincidence.

Here is where chemometrics needs its own tools. The methods below take into account the correlations between variables, so they can see unusual combinations that univariate rules miss.

The Mahalanobis distance generalises "how many standard deviations away" to several correlated variables [13]:

Here ${\mathbf{x}}_i$ is the vector of measurements of sample $i$, $\bar{\mathbf{x}}$ is the mean vector, and $\mathbf{S}$ is the sample covariance matrix, which describes how every pair of variables varies together.

Geometrically, the Mahalanobis distance first stretches and rotates the space so that the elongated data cloud becomes round, and then measures ordinary distance. A point far out along the direction in which the data are naturally spread gets a modest distance; a point the same number of units away in a direction where the data are thin gets a large one.

For multivariate normal data, $D_i^2$ approximately follows a ${\chi }^2$ distribution with $p$ degrees of freedom, where $p$ is the number of variables. Samples with $D_i^2>{\chi }_{p,1-\alpha }^2$ are flagged, typically with $\alpha =0.05$ or $0.01$.

There is a catch for spectroscopy. Computing ${\mathbf{S}}^{-1}$ requires more samples than variables. A data set of 60 spectra with 1000 wavelengths has a singular covariance matrix that cannot be inverted. The solution is to measure distances in a compressed space, and this is exactly what PCA offers.

After fitting a PCA model with $A$ components, each sample is summarised by its scores ${\mathbf{t}}_i=[t_{i1},\dots ,t_{iA}]$, its coordinates along the principal components. Hotelling's $T^2$ is the Mahalanobis distance computed in this score space [9]:

where ${\lambda }_a$ is the variance of the scores on component $a$ (its eigenvalue). Dividing by ${\lambda }_a$ puts every component on the same footing: a large score on a component that naturally varies a lot is less surprising than the same score on a quiet component.

The limit most chemometrics software uses for the calibration samples, at significance level $\alpha$, comes from the F distribution:

where $n$ is the number of samples. For a brand-new sample projected onto an existing model, the factor in front becomes $A(n^2-1)/(n(n-A))$, which gives a slightly wider limit. Strictly, for samples that were part of the model $T^2$ follows a scaled beta distribution, $\frac{(n-1{)}^2}{n}\text{Beta}\left(\frac{A}{2},\frac{n-A-1}{2}\right)$ [17]. The F-based limit is the common, conservative approximation; the difference is small for large calibration sets but noticeable for small ones.

A high $T^2$ means the sample is extreme within the model: it combines the known patterns in an unusual amount. If PC1 describes fat content and PC2 moisture, a sample with a large $T^2$ has an unusually high or low fat content, an unusual fat and moisture combination, or both.

$T^2$ only looks inside the model. The Q residual, also called the squared prediction error (SPE), looks at what the model leaves behind:

where ${\mathbf{e}}_i={\mathbf{x}}_i-{\hat{\mathbf{x}}}_i$ is the difference between the (preprocessed) sample and its reconstruction from $A$ components, ${\hat{\mathbf{x}}}_i={\sum }_{a=1}^A{t}_{ia}{\mathbf{p}}_a$. The vectors ${\mathbf{p}}_a$ are the loadings, the directions of the principal components. In plain words, $Q_i$ is the squared length of the part of the sample that the model cannot describe.

A limit for $Q$ was derived by Jackson and Mudholkar [5]:

where ${\theta }_k={\sum }_{a=A+1}^r{\lambda }_a^k$ for $k=1,2,3$ sums powers of the eigenvalues of the components that were not kept ($r$ is the rank of the centred data, at most $\min (n-1,p)$), $h_0=1-2{\theta }_1{\theta }_3/(3{\theta }_2^2)$, and $c_{\alpha }$ is the upper $(1-\alpha )$ standard normal quantile. It looks intimidating, but every term comes from the discarded eigenvalues. A simpler alternative, popular in multivariate statistical process control (MSPC), approximates the distribution of $Q$ by a scaled chi-squared distribution, $g{\chi }_h^2$, with $g=v/(2m)$ and $h=2m^2/v$, where $m$ and $v$ are the mean and variance of the calibration Q values [18]. When none of this is available, a practical fallback is the 95th or 99th percentile of the Q values of the calibration samples.

A high $Q$ means the sample contains variation the model was never built to describe: a new interferent, a contaminant, an instrumental artefact, or simply a different kind of sample.

The two statistics are easiest to understand by playing with them. The panel below shows a tiny data set with two variables, for example the absorbance at two wavelengths that both respond to the same analyte. A one-component PCA model is the blue line. Drag the orange sample around and watch both statistics react.

Two things are worth noticing. First, the two statistics measure different directions: $T^2$ grows along the model, $Q$ grows away from it. Second, when the suspicious sample is allowed into the model, it pulls the model towards itself and partly hides its own residual. This is the multivariate version of masking, and it is why diagnostics are usually computed on a model built from trusted reference samples, or with robust methods (see below).

Plotting $T^2$ on one axis and $Q$ on the other, with a limit line for each, divides the samples into four regions [8],[10],[18],[19]:

When the goal is a calibration model (ordinary least squares, PCR or PLS), we also care about the response. It is not enough to know whether a sample is unusual; we want to know whether it changes the model. This is the field of influence diagnostics [15].

In ordinary least squares the fitted values are a linear function of the measured responses, $\hat{\mathbf{y}}=\mathbf{Hy}$, with the hat matrix

(the same $({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T$ you met in the normal equations). Here $\mathbf{X}$ includes a leading column of ones for the intercept, so each row is ${\mathbf{x}}_i=[1,x_{i1},\dots ,x_{ip}{]}^T$. The leverage of sample $i$ is the diagonal element

It tells how strongly the fitted value of sample $i$ depends on its own measured response. With an intercept in the model, leverage lies between $1/n$ and 1, and its average is $(p+1)/n$, where $p$ is the number of predictors. A common rule of thumb flags $h_{ii}>2(p+1)/n$ or, more conservatively, $h_{ii}>3(p+1)/n$.

The key insight is that leverage depends only on the predictors, not on the response. A high-leverage point sits in a sparsely populated region of the predictor space, like a lever arm far from the pivot. It has the potential to move the model. Whether it actually does depends on whether its response agrees with the trend of the other samples.

In PCA and PLS models, leverage is computed in the score space and is directly related to $T^2$:

so leverage and $T^2$ are essentially the same diagnostic seen from two angles.

The ordinary residual is $e_i=y_i-{\hat{y}}_i$. Because high-leverage points pull the fit towards themselves, their raw residuals tend to be small. The internally studentized residual corrects for this:

where $s$ is the residual standard deviation of the model. The internally studentized residual has a hidden ceiling: it can never exceed $\sqrt{n-k}$ in absolute value, where $k$ is the number of fitted coefficients. For small calibration sets ($n=10$, $k=2$ gives about 2.8) that makes fixed cut-offs such as $\pm 3$ unreachable. The externally studentized residual goes one step further and estimates the standard deviation without sample $i$, written $s_{(i)}$, so that a large error cannot inflate its own yardstick, and it has no such ceiling:

If the model is correct, $t_i$ follows a t distribution with $n-p-2$ degrees of freedom. In practice, values beyond about $\pm 2.5$ or $\pm 3$ deserve a look.

Cook's distance combines leverage and residual into one measure of influence [4]. It quantifies how much all fitted values move when sample $i$ is removed:

where ${\hat{\mathbf{y}}}_{(i)}$ are the fitted values of the model built without sample $i$, and $k=p+1$ is the number of fitted coefficients. Happily, there is no need to refit the model $n$ times; the same quantity can be written as

so Cook's distance is large only when a sample has both a large residual and a high leverage. Common guidelines are to investigate samples with $D_i>4/n$, and to treat $D_i>1$ as strongly influential. Cook's distance can be large for a good leverage point too, because $h_{ii}/(1-h_{ii})$ becomes large as $h_{ii}$ approaches 1. In the demo below, the Good leverage preset has $D\approx 0.92$, well above $4/n\approx 0.36$, although the slope hardly changes (1.995 against 1.984 without it). Always look at the fit with and without the sample.

The Williams plot shows leverage on the horizontal axis and the externally studentized residual on the vertical axis, with a vertical line at $h^{\ast }=3(p+1)/n$ and horizontal lines at $\pm 2.5$ or $\pm 3$. It is the regression counterpart of the T² against Q plot:

• Low leverage, small residual: an ordinary, well-fitted sample.
• High leverage, small residual: a "good" leverage point. It anchors the line, and agrees with it.
• Low leverage, large residual: a response outlier that does not move the model much.
• High leverage, large residual: a "bad" leverage point, the most dangerous case.

Try it yourself. Drag the orange standard and compare the fitted line with and without it:

Detecting outliers is only half the job. The harder question is what to do about them, and the single most important principle is: investigate before you delete.

As a rough guide, these situations call for further investigation rather than removal:

• High Q, low T². The sample contains variation the model has not seen. Find out what it is before deciding.
• A chemically plausible edge case. A sample at an extreme pH, temperature or concentration may mark the boundary of the model's useful range. Keeping it widens the range; removing it narrows it.
• A single flag in a large data set. With 500 samples, a 99% limit is expected to flag about five perfectly good samples by chance, so one or two flags are no surprise.

And these are good grounds for removal, with the cause written down:

• A documented measurement failure, such as a failing lamp, a blocked probe window or a saturated detector.
• A confirmed preparation error, such as a mislabelled vial or a wrong dilution.
• Several independent diagnostics agree and the raw spectrum shows a clear problem.

All the classical diagnostics share a weakness: they define "normal" using the mean and covariance of the data, and those statistics are themselves distorted by the outliers they are trying to find. Robust methods estimate the centre and spread of the data in ways that resist contamination.

For a single variable, the median and the median absolute deviation replace the mean and the standard deviation:

Multiplying the MAD by 1.4826 turns it into an estimate of the standard deviation for normal data, which gives a robust z-score:

The MAD has a breakdown point of 50%: up to half of the data can be replaced by arbitrary values before the estimate becomes meaningless. The standard deviation has a breakdown point of zero, because a single wild value can make it as large as you like [12]. The same median and MAD can also be used to centre and scale variables robustly before a multivariate analysis (see centering and scaling).

For several variables, the minimum covariance determinant (MCD) estimator looks for the subset of samples (typically a bit more than half) whose covariance matrix has the smallest determinant, in other words the tightest ellipse, and uses that subset to define the centre and shape of the data. The estimator was proposed by Rousseeuw [22] and is computed in practice with the FAST-MCD algorithm [14]. Outliers are then judged against this clean core.

Because the MCD also needs more samples than variables, Hubert, Rousseeuw and Vanden Branden proposed ROBPCA, a robust PCA for high-dimensional data that combines projection pursuit with the MCD [6]. It provides robust scores, a robust score distance (the analogue of $T^2$) and a robust orthogonal distance (the analogue of $Q$), and it sorts the samples into regular observations, good leverage points, orthogonal outliers and bad leverage points.

In regression, least trimmed squares (LTS), introduced by Rousseeuw in 1984, minimises the sum of the $h$ smallest squared residuals instead of all of them, with $h$ typically about $(n+p+1)/2$ [7]. The worst-fitting samples simply do not take part in the fit, so they cannot drag it. Once the robust model is fitted, the excluded samples can be examined one by one.

The examples below simulate NIR-like spectra from two absorbing components, insert three suspicious samples (an extreme concentration, a new interferent, and both together), and compute Hotelling's T² and Q residuals with their limits.

To avoid the masking described above, the PCA model is built only from the first 55 samples, which we trust. The last five (the three suspects and two ordinary samples) are then projected onto it as new samples and judged with the new-sample $T^2$ limit. You should see the three suspects far beyond the limits (the axes are logarithmic) and the two ordinary samples among the calibration samples. The extreme sample is flagged mainly by $T^2$; its $Q$ is only just above the limit, because small errors in the loadings are magnified at such large scores. If you instead build the model with all 60 samples, the suspects distort it: the Q limit becomes hundreds of times larger and every ordinary sample gets a Q far above the noise level.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def fit_pca_model(X_cal, n_components=2, alpha=0.05):
    """
    Build a PCA model on trusted calibration samples and compute
    the limits for Hotelling's T² and the Q residual.

    X_cal : array (n_samples, n_variables), mean-centred internally
    n_components : number of principal components in the model
    alpha : significance level of the limits
    """
    n, p = X_cal.shape
    A = n_components
    mean = X_cal.mean(axis=0)

    # PCA via SVD of the centred calibration data
    U, s, Vt = np.linalg.svd(X_cal - mean, full_matrices=False)
    eigenvalues = s**2 / (n - 1)
    model = {'mean': mean, 'loadings': Vt[:A, :].T, 'eigenvalues': eigenvalues[:A]}

    # T² limits (F distribution): calibration samples and new samples
    F = stats.f.ppf(1 - alpha, A, n - A)
    model['T2_lim_cal'] = A * (n - 1) / (n - A) * F
    model['T2_lim_new'] = A * (n**2 - 1) / (n * (n - A)) * F

    # Q limit (Jackson-Mudholkar), percentile of the calibration Q as a fallback
    eig_res = eigenvalues[A:]
    theta1, theta2, theta3 = (np.sum(eig_res**k) for k in (1, 2, 3))
    _, Q_cal = project(model, X_cal)
    model['Q_lim'] = np.percentile(Q_cal, 100 * (1 - alpha))
    if theta1 > 0 and theta2 > 0:
        h0 = 1 - 2 * theta1 * theta3 / (3 * theta2**2)
        if h0 > 0:
            c_alpha = stats.norm.ppf(1 - alpha)
            model['Q_lim'] = theta1 * (
                c_alpha * np.sqrt(2 * theta2 * h0**2) / theta1
                + 1
                + theta2 * h0 * (h0 - 1) / theta1**2
            ) ** (1 / h0)
    return model

def project(model, X):
    """Hotelling's T² and Q residual of samples X under an existing model."""
    Xc = X - model['mean']
    P = model['loadings']
    T = Xc @ P
    T2 = np.sum(T**2 / model['eigenvalues'], axis=1)
    Q = np.sum((Xc - T @ P.T)**2, axis=1)
    return T2, Q


# Simulated NIR-like spectra
rng = np.random.default_rng(42)
n_samples, n_wavelengths = 60, 200
wavelengths = np.linspace(1000, 2500, n_wavelengths)

comp1 = np.exp(-((wavelengths - 1400)**2) / 20000)
comp2 = np.exp(-((wavelengths - 1900)**2) / 30000)
concentrations = rng.random((n_samples, 2)) * 2

X = concentrations @ np.array([comp1, comp2])
X += rng.normal(0, 0.005, X.shape)

# Three suspicious samples (Python indices start at 0)
comp3 = np.exp(-((wavelengths - 2200)**2) / 10000)
X[55] = 5 * comp1 + 4 * comp2 + rng.normal(0, 0.005, n_wavelengths)  # extreme: high T²
X[56] = 0.8 * comp1 + 1.0 * comp2 + 1.5 * comp3                      # new interferent: high Q
X[57] = 4 * comp1 + 0.5 * comp2 + 2.0 * comp3                        # both

# Build the model on the first 55 (trusted) samples only, then project
# the last five as new samples. Fitting the model with the suspects
# inside would let them mask themselves (see the text above).
X_cal, X_new = X[:55], X[55:]
model = fit_pca_model(X_cal, n_components=2, alpha=0.05)
T2_cal, Q_cal = project(model, X_cal)
T2_new, Q_new = project(model, X_new)

flag_cal = (T2_cal > model['T2_lim_cal']) | (Q_cal > model['Q_lim'])
flag_new = (T2_new > model['T2_lim_new']) | (Q_new > model['Q_lim'])
print(f"T² limit (new samples) = {model['T2_lim_new']:.2f}, Q limit = {model['Q_lim']:.2e}")
print(f"Calibration samples over a limit: {flag_cal.sum()} of {len(X_cal)}")
for k in range(len(X_new)):
    print(f"Sample {55 + k}: T² = {T2_new[k]:.2f}, Q = {Q_new[k]:.2e}"
          + ("  <- flagged" if flag_new[k] else ""))

fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(T2_cal, Q_cal, c='#3b82f6', alpha=0.7, label='Calibration samples')
ax.scatter(T2_new[~flag_new], Q_new[~flag_new], c='#10b981', marker='s', s=60,
           edgecolors='k', label='New, within limits')
ax.scatter(T2_new[flag_new], Q_new[flag_new], c='#ef4444', marker='D', s=80,
           edgecolors='k', label='New, flagged')
for k in np.where(flag_new)[0]:
    ax.annotate(str(55 + k), (T2_new[k], Q_new[k]), textcoords='offset points', xytext=(5, 5))
ax.axvline(model['T2_lim_new'], color='gray', linestyle='--', label='T² limit')
ax.axhline(model['Q_lim'], color='gray', linestyle=':', label='Q limit')
ax.set_xscale('log')   # the suspects are orders of magnitude away
ax.set_yscale('log')
ax.set_xlabel("Hotelling's T²")
ax.set_ylabel('Q residual (SPE)')
ax.set_title('T² against Q')
ax.legend()
plt.tight_layout()
plt.show()

% Simulated NIR-like spectra
% (xline and yline require MATLAB R2018b or later)
rng(42);
n_samples = 60;
n_wavelengths = 200;
wavelengths = linspace(1000, 2500, n_wavelengths);

comp1 = exp(-((wavelengths - 1400).^2) / 20000);
comp2 = exp(-((wavelengths - 1900).^2) / 30000);
concentrations = rand(n_samples, 2) * 2;

X = concentrations * [comp1; comp2] + randn(n_samples, n_wavelengths) * 0.005;

% Three suspicious samples (MATLAB indices start at 1)
comp3 = exp(-((wavelengths - 2200).^2) / 10000);
X(56, :) = 5 * comp1 + 4 * comp2 + randn(1, n_wavelengths) * 0.005;  % high T²
X(57, :) = 0.8 * comp1 + 1.0 * comp2 + 1.5 * comp3;                   % high Q
X(58, :) = 4 * comp1 + 0.5 * comp2 + 2.0 * comp3;                     % both

% Model on the first 55 (trusted) samples; the last five are new samples
X_cal = X(1:55, :);
X_new = X(56:end, :);
model = fit_pca_model(X_cal, 2, 0.05);
[T2_cal, Q_cal] = project_pca(model, X_cal);
[T2_new, Q_new] = project_pca(model, X_new);

flag_cal = (T2_cal > model.T2_lim_cal) | (Q_cal > model.Q_lim);
flag_new = (T2_new > model.T2_lim_new) | (Q_new > model.Q_lim);
fprintf('T^2 limit (new samples) = %.2f, Q limit = %.2e\n', model.T2_lim_new, model.Q_lim);
fprintf('Calibration samples over a limit: %d of %d\n', sum(flag_cal), size(X_cal, 1));
for k = 1:size(X_new, 1)
    fprintf('Sample %d: T^2 = %.2f, Q = %.2e%s\n', 55 + k, T2_new(k), Q_new(k), ...
            repmat('  <- flagged', 1, flag_new(k)));
end

figure;
loglog(T2_cal, Q_cal, 'o', 'MarkerFaceColor', [0.23 0.51 0.96], 'MarkerEdgeColor', 'none'); hold on;
loglog(T2_new(~flag_new), Q_new(~flag_new), 's', 'MarkerFaceColor', [0.06 0.73 0.51], 'MarkerEdgeColor', 'k');
loglog(T2_new(flag_new), Q_new(flag_new), 'd', 'MarkerSize', 9, ...
       'MarkerFaceColor', [0.94 0.27 0.27], 'MarkerEdgeColor', 'k');
idx = find(flag_new);
text(T2_new(idx), Q_new(idx), arrayfun(@(k) sprintf('  %d', 55 + k), idx, 'UniformOutput', false));
xline(model.T2_lim_new, '--', 'Color', [0.5 0.5 0.5]);
yline(model.Q_lim, ':', 'Color', [0.5 0.5 0.5]);
xlabel("Hotelling's T^2"); ylabel('Q residual (SPE)');
title('T^2 against Q'); grid on;
legend('Calibration samples', 'New, within limits', 'New, flagged', 'Location', 'east');

% Local functions must be at the end of a MATLAB script
function model = fit_pca_model(X_cal, A, alpha)
    % PCA model on trusted calibration samples, with T² and Q limits
    n = size(X_cal, 1);
    model.mean = mean(X_cal, 1);
    [~, S, V] = svd(X_cal - model.mean, 'econ');
    eigenvalues = diag(S).^2 / (n - 1);
    model.loadings = V(:, 1:A);
    model.eigenvalues = eigenvalues(1:A)';

    % T² limits: calibration samples and new samples
    F = finv(1 - alpha, A, n - A);
    model.T2_lim_cal = A * (n - 1) / (n - A) * F;
    model.T2_lim_new = A * (n^2 - 1) / (n * (n - A)) * F;

    % Q limit (Jackson-Mudholkar), percentile of the calibration Q as a fallback
    eig_res = eigenvalues(A+1:end);
    theta1 = sum(eig_res); theta2 = sum(eig_res.^2); theta3 = sum(eig_res.^3);
    [~, Q_cal] = project_pca(model, X_cal);
    model.Q_lim = prctile(Q_cal, 100 * (1 - alpha));
    if theta1 > 0 && theta2 > 0
        h0 = 1 - 2 * theta1 * theta3 / (3 * theta2^2);
        if h0 > 0
            c_alpha = norminv(1 - alpha);
            model.Q_lim = theta1 * (c_alpha * sqrt(2 * theta2 * h0^2) / theta1 ...
                          + 1 + theta2 * h0 * (h0 - 1) / theta1^2)^(1 / h0);
        end
    end
end

function [T2, Q] = project_pca(model, X)
    % Hotelling's T² and Q residual of samples X under an existing model
    Xc = X - model.mean;
    P = model.loadings;
    T = Xc * P;
    T2 = sum(T.^2 ./ model.eigenvalues, 2);
    Q = sum((Xc - T * P').^2, 2);
end

project_pca <- function(model, X) {
  # Hotelling's T² and Q residual of samples X under an existing model
  Xc <- sweep(X, 2, model$mean)
  P <- model$loadings
  Tm <- Xc %*% P
  T2 <- rowSums(sweep(Tm^2, 2, model$eigenvalues, "/"))
  Q <- rowSums((Xc - Tm %*% t(P))^2)
  list(T2 = T2, Q = Q)
}

fit_pca_model <- function(X_cal, A = 2, alpha = 0.05) {
  # PCA model on trusted calibration samples, with T² and Q limits
  n <- nrow(X_cal)
  mu <- colMeans(X_cal)
  sv <- svd(sweep(X_cal, 2, mu))
  eigenvalues <- sv$d^2 / (n - 1)
  model <- list(mean = mu, loadings = sv$v[, 1:A, drop = FALSE],
                eigenvalues = eigenvalues[1:A])

  # T² limits: calibration samples and new samples
  Fq <- qf(1 - alpha, A, n - A)
  model$T2_lim_cal <- A * (n - 1) / (n - A) * Fq
  model$T2_lim_new <- A * (n^2 - 1) / (n * (n - A)) * Fq

  # Q limit (Jackson-Mudholkar), percentile of the calibration Q as a fallback
  eig_res <- eigenvalues[(A + 1):length(eigenvalues)]
  theta1 <- sum(eig_res); theta2 <- sum(eig_res^2); theta3 <- sum(eig_res^3)
  model$Q_lim <- unname(quantile(project_pca(model, X_cal)$Q, 1 - alpha))
  if (theta1 > 0 && theta2 > 0) {
    h0 <- 1 - 2 * theta1 * theta3 / (3 * theta2^2)
    if (h0 > 0) {
      c_alpha <- qnorm(1 - alpha)
      model$Q_lim <- theta1 * (c_alpha * sqrt(2 * theta2 * h0^2) / theta1 +
                      1 + theta2 * h0 * (h0 - 1) / theta1^2)^(1 / h0)
    }
  }
  model
}

# Simulated NIR-like spectra
set.seed(42)
n_samples <- 60
n_wavelengths <- 200
wavelengths <- seq(1000, 2500, length.out = n_wavelengths)

comp1 <- exp(-((wavelengths - 1400)^2) / 20000)
comp2 <- exp(-((wavelengths - 1900)^2) / 30000)
concentrations <- matrix(runif(n_samples * 2) * 2, ncol = 2)

X <- concentrations %*% rbind(comp1, comp2) +
     matrix(rnorm(n_samples * n_wavelengths, sd = 0.005), nrow = n_samples)

# Three suspicious samples (R indices start at 1)
comp3 <- exp(-((wavelengths - 2200)^2) / 10000)
X[56, ] <- 5 * comp1 + 4 * comp2 + rnorm(n_wavelengths, sd = 0.005)  # high T²
X[57, ] <- 0.8 * comp1 + 1.0 * comp2 + 1.5 * comp3                   # high Q
X[58, ] <- 4 * comp1 + 0.5 * comp2 + 2.0 * comp3                     # both

# Model on the first 55 (trusted) samples; the last five are new samples
X_cal <- X[1:55, ]
X_new <- X[56:60, ]
model <- fit_pca_model(X_cal, A = 2, alpha = 0.05)
cal <- project_pca(model, X_cal)
new <- project_pca(model, X_new)

flag_cal <- cal$T2 > model$T2_lim_cal | cal$Q > model$Q_lim
flag_new <- new$T2 > model$T2_lim_new | new$Q > model$Q_lim
cat(sprintf("T2 limit (new samples) = %.2f, Q limit = %.2e\n", model$T2_lim_new, model$Q_lim))
cat(sprintf("Calibration samples over a limit: %d of %d\n", sum(flag_cal), nrow(X_cal)))
cat(sprintf("Sample %d: T2 = %.2f, Q = %.2e%s\n", 55 + seq_along(new$T2), new$T2, new$Q,
            ifelse(flag_new, "  <- flagged", "")), sep = "")

plot(cal$T2, cal$Q, log = "xy", pch = 19, col = "#3b82f6",
     xlim = range(c(cal$T2, new$T2)), ylim = range(c(cal$Q, new$Q)),
     xlab = expression("Hotelling's " * T^2), ylab = "Q residual (SPE)",
     main = expression(T^2 * " against Q"))
points(new$T2, new$Q, pch = ifelse(flag_new, 18, 15), cex = 1.8,
       col = ifelse(flag_new, "#ef4444", "#10b981"))
text(new$T2[flag_new], new$Q[flag_new], labels = 55 + which(flag_new), pos = 4)
abline(v = model$T2_lim_new, lty = 2, col = "gray50")
abline(h = model$Q_lim, lty = 3, col = "gray50")
grid()

An outlier is not a verdict, it is a question. The useful questions are concrete: is the sample extreme within the patterns we already know (high $T^2$, high leverage), does it contain something new (high $Q$), and does it actually change the model (high Cook's distance)?

The squared-error methods at the heart of chemometrics are very sensitive to a few bad samples, which is why these diagnostics matter. But the limits are screening tools. The decision to keep or remove a sample belongs to the chemistry and the laboratory records, not to a threshold alone.

Judge suspects against a model built from samples you trust, or the outliers will mask themselves. When you do not yet trust any reference set, start with robust methods (MCD, ROBPCA), which are not fooled by the outliers they look for.