Multiple Linear Regression

A straight-line calibration works beautifully when one measurement determines the answer. Chemistry is rarely that polite. The absorbance at your favourite wavelength also contains a little of every other species that absorbs there; a sensor response drifts with temperature; a chromatographic peak sits on the shoulder of another one.

Multiple linear regression (MLR) is the natural next step: instead of one predictor, the model uses several at once, and it tries to separate their effects. It is still least squares, it is still solved with the normal equations, and it is still linear in the parameters. What changes is that the predictors can now talk to each other, and when they say almost the same thing, the model gets confused in a very specific and instructive way.

Fitting one variable from several others grew out of the same correlation work that gave us simple regression. George Udny Yule, a student and collaborator of Pearson, applied least squares to regression on several variables in 1897 [10] and developed the notation for multiple and partial regression coefficients in the following decade [1]. The idea of estimating the effect of one factor while accounting for the others is already in his partial coefficients. Ronald A. Fisher, working at the Rothamsted agricultural station, then turned regression and the analysis of variance into everyday tools for experimental scientists with his 1925 book Statistical Methods for Research Workers [2]. His key contributions were the analysis of variance and the design of experiments, which let you estimate the effect of each factor (rainfall, fertiliser, temperature) cleanly while accounting for the others.

In chemistry, the need became pressing with instrumental analysis. A spectrum contains dozens or thousands of readings, and each one carries part of the information about the analyte. Building models that use many wavelengths at once, multivariate calibration, became a cornerstone of chemometrics in the 1970s and 1980s [3]. MLR was the first tool tried, and understanding where it succeeds and where it fails is the best way to understand why ridge, PCR and PLS exist.

Imagine a solution containing an analyte A and an interferent B whose absorption bands overlap. You calibrate at 510 nm, the maximum of A. Whatever B absorbs at 510 nm is silently counted as A. If the standards contained no B, every sample with B is overestimated; if they contained some B, samples with more B than the standards are overestimated and samples with less are underestimated.

The fix is to read a second wavelength where B absorbs strongly and A barely does. A model that uses both readings can learn how much of the 510 nm signal belongs to B and subtract it. Try it: change the amount of interferent in the unknown sample and compare the two predictions.

The one-wavelength model is not "wrong" in its own terms: it found the best straight line through standards in which B varied randomly, so it learned the average amount of B and is right only for samples that happen to contain that much. It simply has no information to tell A and B apart. The two-wavelength model has that information, and least squares uses it.

With $p$ predictors, the model becomes:

Each coefficient ${\beta }_j$ tells you how much the response changes when predictor $x_j$ increases by one unit while all the other predictors stay constant. That last phrase is the heart of MLR. In simple regression the slope means one thing. In multiple regression, the meaning of each coefficient depends on which other variables are in the model, because each coefficient describes only the part of $x_j$ that the other predictors cannot explain.

There are two natural ways to set up a spectroscopic calibration, and it helps to know both [4].

The classical (or direct) approach follows the Beer-Lambert law literally: each absorbance is a sum of contributions from every absorbing species,

where $k_{Aj}$ and $k_{Bj}$ are the sensitivities (molar absorptivity times path length) of A and B at wavelength $j$. To use it you need to know, or estimate, the pure-component contributions of every species in the sample.

The inverse approach turns things around and regresses the concentration of the analyte directly on the absorbances:

This is MLR with absorbances as predictors, the model used in the demo above. Its great advantage is that you only need reference values for the analyte you care about; the interferents are handled implicitly, as long as they vary in the calibration samples the way they will vary in future samples. This inverse form is the ancestor of PCR and PLS.

The demo above is calibrated with eight mixtures in which both species vary (concentrations in arbitrary units, absorbances in AU, with a little simulated noise):

Regressing $c_A$ on $a_{510}$ alone gives ${\hat{c}}_A=-0.104+0.941a_{510}$, with $R^2=0.934$. Adding the second wavelength gives

with $R^2$ practically equal to 1. The negative coefficient at 600 nm is the chemistry: one unit of B adds 0.282 AU at 510 nm and 0.686 AU at 600 nm, so B contributes $1.021\times 0.282\approx 0.288$ through the first term and $-0.426\times 0.686\approx -0.292$ through the second. The two cancel, and B disappears from the prediction. Nobody told the model the spectrum of B; least squares found the subtraction because B varied in the calibration set.

Writing the equation for every sample quickly becomes tedious. Matrix notation compresses it. With $n$ samples and $p$ predictors we collect:

• the response vector $\mathbf{y}$ ($n\times 1$), for example the reference concentrations;
• the design matrix $\mathbf{X}$ ($n\times (p+1)$), one row per sample: a 1 for the intercept followed by the predictor values;
• the coefficient vector $\mathbf{\beta }$ ($(p+1)\times 1$);
• the error vector $\mathbf{\epsilon }$ ($n\times 1$).

Here $x_{ij}$ is the value of predictor $j$ in sample $i$.

The criterion is exactly the one we used for a straight line: minimize the sum of squared residuals,

Setting its gradient to zero gives the normal equations, derived step by step in the normal equations entry:

If ${\mathbf{X}}^T\mathbf{X}$ is invertible, the solution is the ordinary least squares (OLS) estimator:

Under the standard assumptions (the model is right, errors have zero mean, constant variance and are uncorrelated), the Gauss-Markov theorem guarantees that OLS has the smallest variance among all linear unbiased estimators. We met this result at the end of your first step into chemometrics. Keep the word unbiased in mind: ridge regression will deliberately give it up.

1. Linearity. The response is a linear function of the predictors (in the parameters). Curvature shows up as structure in the residuals.
2. Full column rank. No predictor is an exact linear combination of the others, and there are more samples than coefficients (in fact $n>p+1$ is needed, otherwise there are no degrees of freedom left to estimate the noise). If this fails, ${\mathbf{X}}^T\mathbf{X}$ is singular and the solution is not unique.
3. Constant variance of the errors (homoscedasticity). If it fails, weighted least squares is the usual remedy.
4. Independent errors. Measurements in a time series often violate this through drift.
5. Normal errors, needed only for exact small-sample confidence intervals and tests, not for the estimates themselves.

Collinearity (or multicollinearity) means that some predictors are strongly correlated with each other, so that one of them can be almost perfectly predicted from the others. In spectroscopy it is the rule, not the exception: absorbances at 500 and 501 nm measure practically the same thing.

What does that do to MLR? The cleanest way to see it is to repeat the same experiment many times with fresh noise and look at where the coefficients land. In the demo below the true model is $y=x_1+x_2$ plus noise, and you control how correlated $x_1$ and $x_2$ are.

With uncorrelated predictors, the estimates gather tightly around the true values. As the correlation approaches 1, the cloud stretches into a long thin ellipse along the line $b_1+b_2=2$. The data know very well what the sum of the coefficients is, because that is what determines the predictions, but they hardly know how to split it between the two predictors. Any split along the ellipse fits the data almost equally well, so noise decides.

The practical consequences:

• Unstable coefficients. Adding or removing one sample can change a coefficient's size or even its sign. Values like +5000 and -4998 that nearly cancel are a classic symptom.
• Inflated standard errors. The variance of each coefficient is equal to ${\sigma }^2$ (the noise variance) times a diagonal element of $({\mathbf{X}}^T\mathbf{X}{)}^{-1}$, which becomes enormous. Coefficients look "not significant" even when the model as a whole predicts well.
• Fragile predictions outside the calibration pattern. Predictions are fine for samples that follow the same correlation structure as the calibration set, but a new sample that breaks it (a new interferent, a baseline shift) can be predicted wildly wrong.

The variance inflation factor of predictor $j$ measures how much the variance of its coefficient is inflated by its correlation with the others:

where $R_j^2$ is the $R^2$ of regressing predictor $j$ on all the other predictors. A VIF of 1 means no inflation; a common rule of thumb flags values above 10 [8]. With two predictors, $R_j^2$ is simply the squared correlation between them, which is the VIF shown in the demo.

The condition number looks at the whole matrix at once. It is the ratio of the largest to the smallest singular value of the predictor matrix,

Its value depends on how the columns are scaled. Belsley, Kuh and Welsch scale each column of the design matrix, intercept column included and without centring, to unit length, and suggested that condition indices around 30 or more indicate harmful collinearity [7]. Software often reports the condition number of the autoscaled predictors instead (as the code below does), which gives different numbers but tells the same story. Note the square root: the condition number of ${\mathbf{X}}^T\mathbf{X}$ is the square of that of $\mathbf{X}$, which is exactly why solving the normal equations explicitly is risky. For a near-infrared spectrum with hundreds of wavelengths and a few dozen samples, $\kappa (\mathbf{X})$ can be astronomically large.

The ordinary $R^2$ never decreases when you add a predictor, even a useless one. The adjusted $R^2$ charges a price for every extra coefficient:

where $n$ is the number of samples and $p$ the number of predictors. If a new variable does not reduce the residuals enough to pay for itself, $R_{\text{adj}}^2$ goes down.

The F-test asks whether the model explains significantly more than just the mean of $y$:

A large F (compared with the $F_{p,n-p-1}$ distribution) says that at least some predictors are useful. With spectra it is almost always significant, and it does not say which predictors matter.

To test whether one coefficient differs from zero:

where ${\hat{\sigma }}^2=S{S}_{\text{res}}/(n-p-1)$ is the estimated noise variance. Under collinearity the standard errors are inflated, so important variables can look insignificant. That is a symptom of the disease, not a failure of the test.

For a calibration, what matters most is how well the model predicts samples it has never seen. The root mean square error of prediction on an independent test set is the most honest number:

When there are not enough samples for a separate test set, cross-validation imitates one: leave a sample (or a group) out, fit the model on the rest, predict the left-out sample, and repeat until every sample has been predicted once. For MLR, leave-one-out cross-validation does not even require refitting. With the hat matrix $\mathbf{H}=\mathbf{X}({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T$, which maps the measured responses onto the fitted ones ($\hat{\mathbf{y}}=\mathbf{Hy}$),

The diagonal element $h_{ii}$ is the leverage of sample $i$: how strongly its own measurement pulls its own fitted value. Samples with leverage close to 1 sit far from the others in predictor space, have a large influence on the model, and receive large cross-validated residuals. The leverages average $(p+1)/n$, and a common rule of thumb flags samples with $h_{ii}>2(p+1)/n$ (or $3(p+1)/n$) for a closer look (see outliers).

One classic way to make MLR work with spectra is to keep only a few wavelengths. If five well-chosen wavelengths carry the information, MLR on those five avoids most of the collinearity.

• Forward selection starts from an empty model and adds, one at a time, the variable that improves it most.
• Backward elimination starts from all variables and removes the least useful one at each step.
• Stepwise selection allows both additions and removals. It is popular but has well-known problems: the result depends on the path taken, and the p-values printed at the end ignore the search that produced them.
• Best subset selection evaluates every possible subset. With $p$ variables there are $2^p$ subsets, which becomes impossible very quickly (200 wavelengths give about $1.6\times 10^{60}$ subsets). Search heuristics such as genetic algorithms, which evolve a population of candidate subsets, are a classic chemometric way around this [11].

The examples simulate a small spectroscopic calibration in which every wavelength responds to the analyte, so the predictors are strongly correlated. They fit MLR, report the usual statistics (with the leave-one-out RMSECV computed from the leverages, as above), and diagnose the collinearity. The final plot fits the same model to the two halves of the data: the coefficients change a lot, some of them even their sign, although both halves come from the same experiment.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error
import matplotlib.pyplot as plt

# --- Simulate spectroscopic calibration data ---
rng = np.random.default_rng(42)
n_samples, n_wavelengths = 50, 20

y = rng.uniform(0.1, 5.0, n_samples)                 # analyte concentration
wavelengths = np.linspace(400, 700, n_wavelengths)
sensitivity = np.exp(-((wavelengths - 550) ** 2) / 5000)
X = np.outer(y, sensitivity) + rng.normal(0, 0.05, (n_samples, n_wavelengths))

# --- Fit multiple linear regression ---
model = LinearRegression().fit(X, y)
y_pred = model.predict(X)

r2 = r2_score(y, y_pred)
n, p = X.shape
r2_adj = 1 - (1 - r2) * (n - 1) / (n - p - 1)
rmsec = np.sqrt(mean_squared_error(y, y_pred))

# Leave-one-out RMSECV from the leverages (no refitting needed)
X_aug = np.column_stack([np.ones(n), X])
h = np.diag(X_aug @ np.linalg.pinv(X_aug))
rmsecv = np.sqrt(np.mean(((y - y_pred) / (1 - h)) ** 2))

print(f"R²:           {r2:.4f}")
print(f"Adjusted R²:  {r2_adj:.4f}")
print(f"RMSEC:        {rmsec:.4f}")
print(f"RMSECV (LOO): {rmsecv:.4f}")

# --- Diagnose collinearity ---
# Condition number of the autoscaled X (Belsley's own convention scales the
# uncentred columns, intercept included, to unit length; the message is the same)
X_scaled = (X - X.mean(axis=0)) / X.std(axis=0, ddof=1)
print(f"\nCondition number of scaled X: {np.linalg.cond(X_scaled):.1f}")

# VIF of each predictor: diagonal of the inverse correlation matrix
corr = np.corrcoef(X, rowvar=False)
vif = np.diag(np.linalg.inv(corr))
print(f"Largest VIF: {vif.max():.0f}")

# Two halves of the same experiment give very different coefficients:
# that instability is the typical symptom of collinearity
coef_a = LinearRegression().fit(X[:25], y[:25]).coef_
coef_b = LinearRegression().fit(X[25:], y[25:]).coef_
print(f"Coefficients that change sign between halves: {np.sum(np.sign(coef_a) != np.sign(coef_b))} of {p}")

plt.figure(figsize=(10, 4))
plt.bar(wavelengths - 4, coef_a, width=8, label="first half")
plt.bar(wavelengths + 4, coef_b, width=8, label="second half")
plt.xlabel("Wavelength (nm)")
plt.ylabel("Coefficient")
plt.title("MLR coefficients from two halves of the data")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

% --- Simulate spectroscopic calibration data ---
rng(42);
n_samples = 50;
n_wavelengths = 20;

y = 0.1 + (5.0 - 0.1) * rand(n_samples, 1);
wavelengths = linspace(400, 700, n_wavelengths);
sensitivity = exp(-((wavelengths - 550).^2) / 5000);
X = y * sensitivity + 0.05 * randn(n_samples, n_wavelengths);

% --- Fit multiple linear regression (QR-based solve) ---
X_aug = [ones(n_samples, 1), X];
beta = X_aug \ y;
y_pred = X_aug * beta;

SS_res = sum((y - y_pred).^2);
SS_tot = sum((y - mean(y)).^2);
R2 = 1 - SS_res / SS_tot;
n = n_samples; p = n_wavelengths;
R2_adj = 1 - (1 - R2) * (n - 1) / (n - p - 1);

% Leave-one-out RMSECV from the leverages (no refitting needed)
H = X_aug * pinv(X_aug);
h = diag(H);
RMSECV = sqrt(mean(((y - y_pred) ./ (1 - h)).^2));

fprintf('R²:          %.4f\n', R2);
fprintf('Adjusted R²: %.4f\n', R2_adj);
fprintf('RMSEC:       %.4f\n', sqrt(SS_res / n));
fprintf('RMSECV (LOO): %.4f\n', RMSECV);

% --- Diagnose collinearity ---
% Condition number of the autoscaled X (Belsley's own convention scales the
% uncentred columns, intercept included, to unit length; the message is the same)
X_scaled = (X - mean(X)) ./ std(X);
fprintf('\nCondition number of scaled X: %.1f\n', cond(X_scaled));
vif = diag(inv(corrcoef(X)));
fprintf('Largest VIF: %.0f\n', max(vif));

% Two halves of the same experiment give very different coefficients:
% that instability is the typical symptom of collinearity
half = 25;
beta_a = [ones(half, 1), X(1:half, :)] \ y(1:half);
beta_b = [ones(n - half, 1), X(half+1:end, :)] \ y(half+1:end);
fprintf('Coefficients that change sign between halves: %d of %d\n', ...
        sum(sign(beta_a(2:end)) ~= sign(beta_b(2:end))), p);

figure;
bar(wavelengths, [beta_a(2:end), beta_b(2:end)]);
legend('first half', 'second half');
xlabel('Wavelength (nm)'); ylabel('Coefficient');
title('MLR coefficients from two halves of the data');
grid on;

# --- Simulate spectroscopic calibration data ---
set.seed(42)
n_samples <- 50
n_wavelengths <- 20

y <- runif(n_samples, 0.1, 5.0)
wavelengths <- seq(400, 700, length.out = n_wavelengths)
sensitivity <- exp(-((wavelengths - 550)^2) / 5000)
X <- outer(y, sensitivity) + matrix(rnorm(n_samples * n_wavelengths, 0, 0.05),
                                    n_samples, n_wavelengths)
colnames(X) <- paste0("wl_", round(wavelengths))

# --- Fit multiple linear regression ---
df <- data.frame(y = y, X)
model <- lm(y ~ ., data = df)
s <- summary(model)

# Leave-one-out RMSECV from the leverages (hat values)
h <- hatvalues(model)
rmsecv <- sqrt(mean((residuals(model) / (1 - h))^2))

cat("R²:           ", s$r.squared, "\n")
cat("Adjusted R²:  ", s$adj.r.squared, "\n")
cat("RMSEC:        ", sqrt(mean(residuals(model)^2)), "\n")
cat("RMSECV (LOO): ", rmsecv, "\n")

# --- Diagnose collinearity ---
# Condition number of the autoscaled X (Belsley's own convention scales the
# uncentred columns, intercept included, to unit length; the message is the same)
cat("\nCondition number of scaled X:", kappa(scale(X), exact = TRUE), "\n")
vif <- diag(solve(cor(X)))
cat("Largest VIF:", round(max(vif)), "\n")

# Two halves of the same experiment give very different coefficients:
# that instability is the typical symptom of collinearity
coef_a <- coef(lm(y ~ ., data = df[1:25, ]))[-1]
coef_b <- coef(lm(y ~ ., data = df[26:50, ]))[-1]
cat("Coefficients that change sign between halves:",
    sum(sign(coef_a) != sign(coef_b)), "of", n_wavelengths, "\n")

barplot(rbind(coef_a, coef_b), beside = TRUE, names.arg = round(wavelengths),
        legend.text = c("first half", "second half"),
        xlab = "Wavelength (nm)", ylab = "Coefficient",
        main = "MLR coefficients from two halves of the data")

MLR works well when:

• you have a few predictors that are not strongly correlated, for example concentration and temperature, or a handful of well-separated wavelengths;
• there are clearly more samples than coefficients ($n\gg p$);
• the predictors come from a designed experiment: two-level factorial designs make the coded factors orthogonal by construction, and good response-surface designs keep the correlations small.

It struggles when:

• the predictors are whole spectra, with many strongly correlated wavelengths;
• there are more variables than samples ($p>n$): ${\mathbf{X}}^T\mathbf{X}$ is then singular and infinitely many coefficient vectors fit the data perfectly;
• the predictors themselves are noisy, so that unstable coefficients amplify that noise into the predictions.

MLR is the straight line generalized to several predictors: one coefficient per variable, all estimated together by least squares, each one describing the effect of its variable with the others held fixed. When predictors carry independent information, as in the interference example, this is exactly what lets the model separate overlapping signals.

When predictors carry nearly the same information, the same mechanism backfires. The data still determine the predictions, but not the individual coefficients, which become large, unstable and uninterpretable. Recognizing this is what separates naive regression from chemometric modelling, and every method that follows is a different answer to it:

• Ridge regression adds a penalty $\lambda \Vert \mathbf{\beta }{\Vert }^2$ (usually not applied to the intercept) that keeps the coefficients small and stable.
• Lasso regression uses an L1 penalty $\lambda \Vert \mathbf{\beta }{\Vert }_1$ that sets some coefficients exactly to zero.
• Principal component regression replaces the correlated predictors with a few uncorrelated principal components.
• Partial least squares builds components that are relevant for predicting $y$, and is the workhorse of chemometric calibration.

For a thorough treatment of regression diagnostics, see Draper and Smith [5] or Montgomery, Peck and Vining [6]; for the chemometric view, Brereton [9] and Martens and Næs [3].