Multiple Linear Regression

The idea of fitting a straight line to data goes back centuries, but chemistry rarely gives us the luxury of a single predictor. When a spectrophotometer records absorbance at 200 wavelengths and we want to predict the concentration of an analyte, we need a model that can handle all 200 variables at once. That model is multiple linear regression (MLR).

The mathematical machinery behind MLR — the matrix formulation of least squares — was assembled piece by piece during the first half of the twentieth century. Carl Friedrich Gauss had already described the method of least squares around 1795, and Adrien-Marie Legendre published the first formal account in 1805. But it was the statisticians of the early 1900s who generalized the idea to multiple predictors. R.A. Fisher, working at the Rothamsted Experimental Station in the 1920s and 1930s, developed the analysis of variance and the general linear model framework that placed multiple regression on rigorous statistical footing. His 1925 book Statistical Methods for Research Workers gave scientists a practical toolbox for handling several variables simultaneously [1].

Fisher’s framework was not purely theoretical. He applied multiple regression to agricultural field trials at Rothamsted, where crop yield depended on rainfall, temperature, soil nitrogen, and other factors simultaneously. His insight was that you could estimate the effect of each factor while accounting for all the others — a problem that single-variable analysis could not solve. The mathematical tools he developed — matrix algebra applied to linear models, the F-distribution for testing groups of coefficients, the concept of degrees of freedom — became the foundation that chemometricians would later build upon [3].

In chemistry, the need for multiple regression became acute with the rise of instrumental analysis in the mid-twentieth century. A UV-Vis spectrum contains absorbance readings at dozens or hundreds of wavelengths, and each wavelength carries partial information about the analyte concentration. Using just one wavelength — simple linear regression — throws away most of the available data. Multivariate calibration, the practice of building models that use many wavelengths to predict concentration, became a cornerstone of chemometrics by the 1970s and 1980s [2]. MLR was the natural first attempt, and understanding where it succeeds and where it fails is essential for appreciating the more advanced methods (PCR, PLS) that came later.

Why one variable is not enough

In the simple linear regression article, we modeled a response $y$ as a function of a single predictor $x$ . That works when one measurement clearly determines the outcome — for example, absorbance at a single wavelength predicting a single analyte via Beer’s Law.

But consider a more realistic scenario. You have a mixture of three dyes in solution, and you measure the UV-Vis spectrum at 50 wavelengths. Each wavelength responds to all three dyes, with different sensitivities. No single wavelength can disentangle the contributions. You need a model that combines information from multiple wavelengths to predict each dye’s concentration.

This is where multiple linear regression enters. Instead of one slope and one intercept, we have many slopes — one for each predictor variable — plus an intercept:

$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdots +{\beta }_p{x}_p+\epsilon$$

Each coefficient ${\beta }_j$ tells us how much the response changes when predictor $x_j$ increases by one unit, holding all other predictors constant. That last phrase is critical. In simple regression, the slope has an unambiguous interpretation. In multiple regression, the interpretation of each coefficient depends on what other variables are in the model.

Spectroscopy example

In a calibration model using 50 wavelengths, ${\beta }_{17}$ tells you how much the predicted concentration changes when absorbance at wavelength 17 increases by one unit, assuming absorbance at all other wavelengths stays the same. In practice, neighboring wavelengths are highly correlated, so this “holding constant” interpretation becomes slippery — a preview of the multicollinearity problem we will address shortly.

A concrete example: multi-analyte determination

Suppose you are analyzing pharmaceutical tablets that contain two active ingredients, A and B. You measure the UV-Vis spectrum at five wavelengths. Beer’s Law tells us that absorbance is additive:

$$a_j={\epsilon }_{Aj}\cdot c_A\cdot l+{\epsilon }_{Bj}\cdot c_B\cdot l+\text{noise}$$

where $a_j$ is the absorbance at wavelength $j$ , ${\epsilon }_{Aj}$ and ${\epsilon }_{Bj}$ are the molar absorptivities of A and B at that wavelength, $c_A$ and $c_B$ are the concentrations, and $l$ is the path length.

To determine $c_A$ from the spectrum, simple regression at one wavelength would fail if B also absorbs there. But with multiple wavelengths, we can set up an MLR model where the absorbances at all five wavelengths are the predictors and the known concentration of A (from reference analyses) is the response. The model learns a set of coefficients that, when applied to a new spectrum, disentangle the contributions of A and B and return the concentration of A. This is the essence of multivariate calibration.

The model in matrix form

Writing out $y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdots$ for every observation and every predictor quickly becomes tedious. Matrix notation compresses everything into a single equation.

Suppose we have $n$ observations and $p$ predictor variables. We arrange the data into:

• A response vector $\mathbf{y}$ of size $n\times 1$ , containing the measured values (e.g., concentrations)
• A design matrix $\mathbf{X}$ of size $n\times (p+1)$ , where each row is one observation. The first column is all ones (for the intercept), and the remaining $p$ columns hold the predictor values (e.g., absorbances at different wavelengths)
• A coefficient vector $\mathbf{\beta }$ of size $(p+1)\times 1$ , containing the intercept and all slopes
• An error vector $\mathbf{\epsilon }$ of size $n\times 1$ , containing the noise

The model is:

$$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{\epsilon }$$

Written out, the design matrix looks like:

$$\mathbf{X}=\left(\begin{array}{ccccc}1& x_{11}& x_{12}& \cdots & x_{1p}\\ 1& x_{21}& x_{22}& \cdots & x_{2p}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n1}& x_{n2}& \cdots & x_{np}\end{array}\right)$$

where $x_{ij}$ is the value of predictor $j$ for observation $i$ . The column of ones ensures that when we multiply $\mathbf{X}\mathbf{\beta }$ , the first element of $\mathbf{\beta }$ acts as the intercept ${\beta }_0$ .

Why matrices?

Beyond notational convenience, the matrix form reveals the geometry of regression. The columns of $\mathbf{X}$ span a subspace, and the fitted values $\hat{\mathbf{y}}=\mathbf{X}\hat{\mathbf{\beta }}$ are the orthogonal projection of $\mathbf{y}$ onto that subspace. The residual vector $\mathbf{e}=\mathbf{y}-\hat{\mathbf{y}}$ is perpendicular to every column of $\mathbf{X}$ . This geometric view makes many results (like $R^2$ being a cosine squared) intuitive.

The normal equations

Just as in simple regression, we want to minimize the sum of squared residuals:

$$S(\mathbf{\beta })=\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2=(\mathbf{y}-\mathbf{X}\mathbf{\beta }{)}^T(\mathbf{y}-\mathbf{X}\mathbf{\beta })$$

To find the minimum, we take the derivative with respect to $\mathbf{\beta }$ and set it to zero. Expanding the expression:

$$S(\mathbf{\beta })={\mathbf{y}}^T\mathbf{y}-2{\mathbf{\beta }}^T{\mathbf{X}}^T\mathbf{y}+{\mathbf{\beta }}^T{\mathbf{X}}^T\mathbf{X}\mathbf{\beta }$$

Differentiating with respect to $\mathbf{\beta }$ and setting the gradient to zero:

$$\frac{\partial S}{\partial \mathbf{\beta }}=-2{\mathbf{X}}^T\mathbf{y}+2{\mathbf{X}}^T\mathbf{X}\mathbf{\beta }=0$$

This gives us the normal equations:

$${\mathbf{X}}^T\mathbf{X}\hat{\mathbf{\beta }}={\mathbf{X}}^T\mathbf{y}$$

If ${\mathbf{X}}^T\mathbf{X}$ is invertible (we will see shortly when it is not), the solution is:

$$\hat{\mathbf{\beta }}=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$$

This is the ordinary least squares (OLS) estimator. The matrix $({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T$ is sometimes called the pseudo-inverse of $\mathbf{X}$ (denoted ${\mathbf{X}}^{+}$ ), and it appears throughout chemometrics.

Never actually invert the matrix

In practice, you should not compute $({\mathbf{X}}^T\mathbf{X}{)}^{-1}$ explicitly. It is numerically unstable, especially when the predictors are correlated. Instead, use a QR decomposition or singular value decomposition (SVD) to solve the normal equations. Every major numerical library does this automatically when you call its least-squares solver.

What the solution means

The OLS estimator has a clean interpretation: $\hat{\mathbf{\beta }}$ is the vector of coefficients that makes the predicted values $\hat{\mathbf{y}}=\mathbf{X}\hat{\mathbf{\beta }}$ as close as possible to the observed values $\mathbf{y}$ , measured by Euclidean distance (sum of squared differences).

Under the standard assumptions (errors are independent, identically distributed with mean zero and constant variance ${\sigma }^2$ ), OLS has a remarkable property: the Gauss-Markov theorem guarantees that among all linear unbiased estimators, OLS has the smallest variance. In other words, no other unbiased linear method can do better. This result is what makes OLS the default starting point for regression.

Assumptions of MLR

The Gauss-Markov optimality result depends on several assumptions. When they hold, OLS is hard to beat. When they are violated, the consequences range from mild inefficiency to completely misleading results.

1. Linearity: The true relationship between $\mathbf{y}$ and $\mathbf{X}$ is linear. If the relationship is curved, the residuals will show systematic patterns. In spectroscopy, Beer’s Law is linear only within a concentration range — at high concentrations, deviations from linearity appear.

2. Full rank: The design matrix $\mathbf{X}$ has rank $p+1$ (no predictor is a perfect linear combination of the others). If this fails, ${\mathbf{X}}^T\mathbf{X}$ is singular and cannot be inverted. Near-violations (multicollinearity) produce nearly singular matrices with all the associated instabilities.

3. Homoscedasticity: The error variance is constant across all observations: $\text{Var}({\epsilon }_i)={\sigma }^2$ for all $i$ . In analytical chemistry, measurement precision often decreases at low concentrations (higher relative noise), violating this assumption. Weighted least squares is the standard remedy.

4. Independence: The errors ${\epsilon }_i$ are mutually independent. Time-series data (e.g., process monitoring) often violates this through autocorrelation. For calibration with independent samples, this assumption is usually reasonable.

5. Normal errors (for inference): The Gauss-Markov theorem requires only assumptions 1-4 for OLS to be the best linear unbiased estimator. But to construct confidence intervals and perform hypothesis tests (F-test, t-tests), we additionally need ${\epsilon }_i\sim \mathcal{N}(0,{\sigma }^2)$ . For large sample sizes, the central limit theorem provides approximate normality even when the errors are not exactly Gaussian.

Assumptions in practice

In chemometric calibration with spectroscopic data, assumptions 1, 3, and 4 are usually reasonable, but assumption 2 is almost always violated (the wavelengths are highly correlated). This is why MLR is often a stepping stone rather than the final model in spectroscopy.

The multicollinearity problem

Everything described so far works beautifully — as long as ${\mathbf{X}}^T\mathbf{X}$ is invertible and well-conditioned. In many chemistry applications, it is neither.

Multicollinearity occurs when two or more predictor variables are highly correlated with each other. In spectroscopy, this is the rule rather than the exception: absorbance at 500 nm and absorbance at 501 nm are measuring nearly the same thing. They are almost perfectly correlated.

When predictors are highly correlated, ${\mathbf{X}}^T\mathbf{X}$ becomes nearly singular. Its eigenvalues span many orders of magnitude, and the inverse $({\mathbf{X}}^T\mathbf{X}{)}^{-1}$ amplifies small perturbations into enormous changes in the coefficients. The practical consequences are severe:

• Wildly unstable coefficients: Adding or removing a single observation can flip the sign of a coefficient. One coefficient might be +5000 and its neighbor -4998, with the two nearly canceling out. The individual values are meaningless even though their combination might still predict reasonably.

• Inflated standard errors: The variance of each coefficient estimate is proportional to the corresponding diagonal element of $({\mathbf{X}}^T\mathbf{X}{)}^{-1}$ . When the matrix is nearly singular, these diagonal elements are enormous, making every coefficient statistically insignificant even though the model as a whole may fit well.

• Poor prediction on new data: A model with wild coefficients is essentially memorizing the noise in the training data. It will generalize poorly to new samples.

Diagnosing multicollinearity

The condition number of ${\mathbf{X}}^T\mathbf{X}$ is the primary diagnostic:

$$\kappa ({\mathbf{X}}^T\mathbf{X})=\frac{{\lambda }_{\max }}{{\lambda }_{\min }}$$

where ${\lambda }_{\max }$ and ${\lambda }_{\min }$ are the largest and smallest eigenvalues. Rules of thumb:

• $\kappa <30$ : Multicollinearity is not a practical problem
• $30<\kappa <300$ : Moderate multicollinearity, coefficients becoming unreliable
• $\kappa >300$ : Severe multicollinearity, OLS should not be trusted

For a typical NIR spectrum with 1000 wavelengths and 50 samples, the condition number can easily exceed $10^{10}$ . In such cases, OLS is hopeless and we must turn to regularized methods.

Another useful diagnostic is the Variance Inflation Factor (VIF) for each predictor:

$${\text{VIF}}_j=\frac{1}{1-R_j^2}$$

where $R_j^2$ is the $R^2$ obtained by regressing predictor $j$ on all other predictors. A VIF above 10 is traditionally considered problematic.

The chemometrics perspective

Multicollinearity is not a bug in spectroscopic data — it is a fundamental property of the data. Neighboring wavelengths are correlated because they respond to overlapping molecular transitions. The solution is not to remove “redundant” wavelengths (you would lose information), but to use methods designed for correlated predictors: ridge regression, PCR, or PLS.

Model evaluation

Once you have fitted an MLR model, you need to assess its quality. Several complementary metrics exist.

Adjusted R-squared

The ordinary $R^2$ measures the proportion of variance explained:

$$R^2=1-\frac{S{S}_{\text{res}}}{S{S}_{\text{tot}}}=1-\frac{{\sum }_{i=1}^n(y_i-{\hat{y}}_i{)}^2}{{\sum }_{i=1}^n(y_i-\bar{y}{)}^2}$$

But $R^2$ always increases (or stays the same) when you add more predictors, even useless ones. The adjusted $R^2$ penalizes for model complexity:

$$R_{\text{adj}}^2=1-\frac{S{S}_{\text{res}}/(n-p-1)}{S{S}_{\text{tot}}/(n-1)}=1-(1-R^2)\frac{n-1}{n-p-1}$$

where $n$ is the number of observations and $p$ is the number of predictors. Adding a useless variable increases $p$ without reducing $S{S}_{\text{res}}$ enough, so $R_{\text{adj}}^2$ decreases — signaling that the variable did not help.

The F-test

The overall F-test asks: “Does this model explain significantly more variance than a model with no predictors (just the mean)?”

$$F=\frac{(S{S}_{\text{tot}}-S{S}_{\text{res}})/p}{S{S}_{\text{res}}/(n-p-1)}=\frac{M{S}_{\text{reg}}}{M{S}_{\text{res}}}$$

A large F-value (relative to the $F_{p,n-p-1}$ distribution) indicates that at least some of the predictors are useful. In chemometrics, where we often have many correlated predictors, the overall F-test is almost always significant — it does not tell you which predictors matter.

Individual t-tests

To test whether a specific coefficient ${\beta }_j$ is significantly different from zero:

$$t_j=\frac{{\hat{\beta }}_j}{\text{SE}({\hat{\beta }}_j)}$$

where $\text{SE}({\hat{\beta }}_j)=\sqrt{{\hat{\sigma }}^2[({\mathbf{X}}^T\mathbf{X}{)}^{-1}{]}_{jj}}$ . Under multicollinearity, the standard errors are inflated (as discussed above), so coefficients that are genuinely important may appear insignificant. This is another symptom of the disease, not a failure of the test itself.

RMSEP and cross-validation

For predictive models, the most honest assessment is the root mean square error of prediction (RMSEP) computed on data that was not used to build the model:

$$\text{RMSEP}=\sqrt{\frac{1}{n_{\text{test}}}\sum _{i=1}^{n_{\text{test}}}(y_i-{\hat{y}}_i{)}^2}$$

Leave-one-out cross-validation (LOOCV) or k-fold cross-validation provides a practical estimate of RMSEP when a separate test set is not available. In LOOCV, each sample is left out once, the model is built on the remaining $n-1$ samples, and the left-out sample is predicted. The process repeats for all $n$ samples, yielding $n$ prediction errors. For MLR, there is an elegant shortcut: the LOOCV RMSE can be computed from a single model fit using the hat matrix $\mathbf{H}=\mathbf{X}({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T$ :

$$\text{RMSECV}=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(\frac{y_i-{\hat{y}}_i}{1-h_{ii}}\right)}^2}$$

where $h_{ii}$ is the $i$ -th diagonal element of $\mathbf{H}$ , known as the leverage of observation $i$ . Samples with high leverage (close to 1) have a disproportionate influence on the model and produce large cross-validated residuals.

Code implementation

Python

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

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

# True concentrations
y = np.random.uniform(0.1, 5.0, n_samples)

# Generate spectra: each wavelength responds to concentration
# with some overlap (correlated predictors)
wavelengths = np.linspace(400, 700, n_wavelengths)
X = np.zeros((n_samples, n_wavelengths))
for j in range(n_wavelengths):
    # Each wavelength has a sensitivity that depends on position
    sensitivity = np.exp(-((wavelengths[j] - 550) ** 2) / 5000)
    X[:, j] = sensitivity * y + np.random.normal(0, 0.05, n_samples)

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

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

print(f"R²:          {r2:.4f}")
print(f"Adjusted R²: {r2_adj:.4f}")
print(f"RMSE (train): {rmse_train:.4f}")

# Cross-validated RMSE
cv_scores = cross_val_score(model, X, y,
                            cv=5, scoring='neg_mean_squared_error')
rmsecv = np.sqrt(-cv_scores.mean())
print(f"RMSECV:      {rmsecv:.4f}")

# --- Diagnose multicollinearity ---
XtX = X.T @ X
eigenvalues = np.linalg.eigvalsh(XtX)
condition_number = eigenvalues[-1] / eigenvalues[0]
print(f"\nCondition number: {condition_number:.1f}")

# Coefficient plot
plt.figure(figsize=(10, 4))
plt.bar(range(n_wavelengths), model.coef_)
plt.xlabel("Wavelength index")
plt.ylabel("Coefficient value")
plt.title("MLR Coefficients — note erratic values from multicollinearity")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

MATLAB

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

% True concentrations
y = 0.1 + (5.0 - 0.1) * rand(n_samples, 1);

% Generate spectra with correlated predictors
wavelengths = linspace(400, 700, n_wavelengths);
X = zeros(n_samples, n_wavelengths);
for j = 1:n_wavelengths
    sensitivity = exp(-((wavelengths(j) - 550)^2) / 5000);
    X(:, j) = sensitivity * y + 0.05 * randn(n_samples, 1);
end

% --- Fit multiple linear regression ---
% Add intercept column
X_aug = [ones(n_samples, 1), X];
beta = X_aug \ y;              % Solve normal equations via QR
y_pred = X_aug * beta;

% Training metrics
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);
RMSE_train = sqrt(SS_res / n);

fprintf('R²:           %.4f\n', R2);
fprintf('Adjusted R²:  %.4f\n', R2_adj);
fprintf('RMSE (train): %.4f\n', RMSE_train);

% --- Diagnose multicollinearity ---
XtX = X' * X;
eigenvalues = eig(XtX);
condition_number = max(eigenvalues) / min(eigenvalues);
fprintf('\nCondition number: %.1f\n', condition_number);

% Coefficient plot (excluding intercept)
figure;
bar(beta(2:end));
xlabel('Wavelength index');
ylabel('Coefficient value');
title('MLR Coefficients');
grid on;

R

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

# True concentrations
y <- runif(n_samples, 0.1, 5.0)

# Generate spectra with correlated predictors
wavelengths <- seq(400, 700, length.out = n_wavelengths)
X <- matrix(0, nrow = n_samples, ncol = n_wavelengths)
for (j in 1:n_wavelengths) {
  sensitivity <- exp(-((wavelengths[j] - 550)^2) / 5000)
  X[, j] <- sensitivity * y + rnorm(n_samples, 0, 0.05)
}
colnames(X) <- paste0("wl_", 1:n_wavelengths)

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

cat("R²:          ", summary_model$r.squared, "\n")
cat("Adjusted R²: ", summary_model$adj.r.squared, "\n")
cat("RMSE (train):", sqrt(mean(model$residuals^2)), "\n")

# --- Diagnose multicollinearity ---
XtX <- t(X) %*% X
eigenvalues <- eigen(XtX)$values
condition_number <- max(eigenvalues) / min(eigenvalues)
cat("\nCondition number:", condition_number, "\n")

# VIF (requires car package)
# library(car)
# vif(model)

# Coefficient plot (excluding intercept)
coefs <- coef(model)[-1]
barplot(coefs, names.arg = 1:n_wavelengths,
        xlab = "Wavelength index", ylab = "Coefficient value",
        main = "MLR Coefficients",
        col = ifelse(coefs > 0, "steelblue", "tomato"))
grid()

Variable selection

One common strategy for making MLR work with spectroscopic data is to select a small subset of wavelengths before fitting the model. If we can identify, say, 5 wavelengths out of 200 that carry most of the information about the analyte, then MLR on those 5 wavelengths avoids multicollinearity entirely.

Several approaches exist:

Forward selection starts with no predictors and adds them one at a time, choosing at each step the variable that most improves the model (highest partial F-statistic or lowest AIC). The process stops when no remaining variable provides a statistically significant improvement.

Backward elimination starts with all predictors and removes them one at a time, dropping the least useful variable at each step.

Stepwise selection combines both: at each step, variables can be added or removed. This is the most common approach in practice, though it has well-known problems — the selected variables depend heavily on the order of addition and removal, and the resulting p-values are not valid because of the multiple testing inherent in the search.

Best subset selection evaluates all $2^p$ possible subsets of predictors and picks the one that optimizes some criterion (AIC, BIC, or cross-validated RMSE). This is optimal in principle but computationally infeasible when $p$ is large (with 200 wavelengths, there are more subsets than atoms in the universe).

Why variable selection is tricky for spectra

Variable selection for spectroscopic data is inherently unstable. Because neighboring wavelengths are highly correlated, the selection algorithm might pick wavelength 502 nm in one run and 503 nm in another — both carry nearly the same information. This instability does not mean the resulting models are bad, but it makes interpretation difficult and highlights why regularization methods (ridge, LASSO) or latent variable methods (PCR, PLS) are generally preferred for spectroscopic calibration.

When to use and limitations

MLR works well when

Few, uncorrelated predictors When you have a small number of predictors that are not strongly correlated, OLS gives excellent, interpretable results.

More observations than predictors The system $n>p$ is required for OLS to have a unique solution, and $n\gg p$ is preferred for stable estimates.

Well-designed experiments Factorial or response surface designs produce orthogonal predictors, eliminating multicollinearity by construction.

Variable selection has been done If you have already selected a small set of informative wavelengths (e.g., via stepwise selection or expert knowledge), MLR can work well.

MLR struggles when

Spectroscopic data (many correlated wavelengths) This is the most common situation in chemometrics, and it is precisely where OLS fails. The condition number explodes and coefficients become meaningless.

More variables than observations (p > n) The system is underdetermined and ${\mathbf{X}}^T\mathbf{X}$ is singular. OLS has infinitely many solutions.

Noisy predictors When measurement noise in the predictors is comparable to the signal, OLS amplifies the noise through unstable coefficients.

No variable selection Throwing all available wavelengths into an MLR model almost always produces poor results for spectroscopic data.

Next steps

Multiple linear regression is the starting point for understanding every regression method in chemometrics. Its failure modes — particularly multicollinearity and the $p>n$ problem — motivate every technique that follows:

• Ridge regression: Adds a penalty term $\lambda \Vert \mathbf{\beta }{\Vert }^2$ to stabilize the coefficients. The simplest fix for multicollinearity.
• LASSO regression: Uses an L1 penalty $\lambda \Vert \mathbf{\beta }{\Vert }_1$ that drives some coefficients exactly to zero, performing automatic variable selection.
• Principal component regression (PCR): Replaces the original correlated predictors with a smaller set of uncorrelated principal components, then regresses on those.
• Partial least squares (PLS): Like PCR, but constructs components that maximize covariance with the response, not just variance in the predictors. The workhorse of chemometric calibration.

The key insight to carry forward: OLS is optimal under the Gauss-Markov conditions, but those conditions are violated in almost every spectroscopic dataset. Recognizing this is what separates naive regression from chemometric modeling.

References

[1] Fisher, R.A. (1925). Statistical Methods for Research Workers. Oliver and Boyd, Edinburgh.

[2] Martens, H., & Naes, T. (1989). Multivariate Calibration. John Wiley & Sons.

[3] Draper, N.R., & Smith, H. (1998). Applied Regression Analysis (3rd ed.). John Wiley & Sons.

[4] Wold, S., Sjostrom, M., & Eriksson, L. (2001). “PLS-regression: a basic tool of chemometrics.” Chemometrics and Intelligent Laboratory Systems, 58(2), 109-130.

[5] Brereton, R.G. (2003). Chemometrics: Data Analysis for the Laboratory and Chemical Plant. John Wiley & Sons.

[6] Montgomery, D.C., Peck, E.A., & Vining, G.G. (2012). Introduction to Linear Regression Analysis (5th ed.). John Wiley & Sons.