Linear Regression

If you have ever prepared a set of standards, measured them, and drawn a straight line through the points, you have already built a linear regression model. It is probably the most used model in analytical chemistry, and it is also the seed from which almost every regression method in chemometrics grows.

In this entry we will look at it carefully: where the name comes from, what the model really says, how the best line is found, how to judge whether the line is any good, and where it quietly fails. Along the way you will be able to drag calibration points around and see the least squares machinery react in real time.

The word regression entered statistics through an unlikely door: the study of human height. In 1886, the English polymath Francis Galton published Regression towards Mediocrity in Hereditary Stature [1], based on the heights of 928 adult children from 205 families and their parents. He noticed that very tall parents tended to have tall children, but not quite as tall as themselves, and very short parents tended to have short children, but not quite as short. The children's heights "regressed" towards the population average.

Galton drew a straight line through that cloud of parent and child heights, and the name stuck to the line, even though today it has nothing to do with going backwards. Galton himself introduced the idea of correlation in 1888 [11], and his colleague Karl Pearson gave it its mathematical form in 1896, including the product-moment formula for the correlation coefficient $r$ that we still use [2].

The engine underneath, the method of least squares, is even older and comes from astronomy. Adrien-Marie Legendre published it in 1805 while computing the orbits of comets [3], and Carl Friedrich Gauss gave it a probabilistic justification in 1809 [4]. Gauss also claimed he had been using the method since 1795, which started one of the most famous priority disputes in the history of science [5].

For chemists, the decisive step came from the laws of light absorption. The Beer-Lambert law says that absorbance is proportional to concentration. Bouguer (1729) and Lambert (1760) had described how absorption depends on the thickness of the sample; Beer added the concentration in 1852 [6]. That proportionality turned a straight line into the natural way of converting an instrument reading into an amount of substance: the calibration curve.

Suppose you need to determine iron in drinking water. You complex Fe(II) with 1,10-phenanthroline, which gives an orange-red complex that absorbs strongly near 510 nm, and you measure absorbance in a UV-Vis spectrophotometer. Beer-Lambert tells you that absorbance should be proportional to concentration, but you do not know the exact proportionality constant for your instrument, cuvette and reagent batch. You need a calibration.

Step 3 is the linear regression. The question is what "best" means, and the answer is the same criterion we derived in your first step into chemometrics: the line that minimizes the sum of the squared vertical distances between the points and the line.

The simple linear regression model relates a response $y$ to a single predictor $x$:

Every symbol has a concrete meaning in the lab:

• $y_i$ is the measured response of the $i$-th standard (the absorbance).
• $x_i$ is the known value of the predictor (the iron concentration).
• ${\beta }_0$ is the intercept, the response expected when $x=0$. Ideally it is zero for a blank, but it is often slightly off because of reagent absorption, stray light or a detector offset.
• ${\beta }_1$ is the slope, how much the response changes per unit of concentration. Analysts call it the sensitivity of the method.
• ${\epsilon }_i$ is the error term: the part of the measurement that the straight line cannot explain, mostly random noise.

Notice a subtle but important detail: the model puts the noise in $y$ only. It assumes that the concentrations of the standards are known much more precisely than the absorbances are measured. In a calibration with carefully prepared standards this is usually true, and it is the reason why the residuals are measured vertically.

For any candidate line, the residual of the $i$-th point is the vertical gap between what we measured and what the line predicts:

Least squares chooses the intercept and slope that make the sum of squared residuals as small as possible:

Try it below with the six iron standards used in the examples later in this entry. Drag the standards and watch the fitted line, the residuals, and the quality numbers update. Switch on the squared residuals to see the quantity that is being minimized, and then try to beat least squares with your own line.

You probably noticed two things. First, you can never get a smaller sum of squares than the least squares line (that is the whole point of it). Second, one badly placed standard produces a huge square, so it pulls the line strongly towards itself. Squares punish large errors much more than small ones, which is both the strength and the weakness of the method.

To find the minimum, we take the partial derivatives of $S$ with respect to both parameters and set them to zero, exactly as in the normal equations entry:

These two equations are the normal equations for a straight line. Solving them gives:

where $\bar{x}$ and $\bar{y}$ are the averages of the concentrations and the responses. The hat on $\hat{\beta }$ reminds us that these are estimates computed from data, not the true (unknown) values.

The slope formula has a friendly reading. The numerator measures how $x$ and $y$ vary together (it is proportional to their covariance); the denominator measures how much $x$ varies on its own (proportional to its variance). If both go up together, the slope is positive. The second equation tells us something nice too: the least squares line always passes through the centre of the data, the point $(\bar{x},\bar{y})$.

A fitted line always exists, even through a cloud of random points. So after fitting, we need numbers and, above all, pictures that tell us whether the model describes the data.

The coefficient of determination $R^2$ measures the fraction of the variation in $y$ that the line explains:

Here $S{S}_{\text{res}}$ (the residual sum of squares) is the variation left over after fitting, and $S{S}_{\text{tot}}$ (the total sum of squares) is the variation of the responses around their mean, the variation you would have if you ignored $x$ entirely. For a straight line with an intercept, $R^2$ lies between 0 and 1, and it equals the square of Pearson's correlation coefficient $r$.

A well-behaved Beer-Lambert calibration easily gives $R^2>0.999$. That sounds reassuring, but it says surprisingly little, as the next demo shows.

$R^2$ has no units. The root mean square error (RMSE) gives the typical size of a residual in the units of the response:

An RMSE of 0.002 AU means that the standards typically sit about 0.002 absorbance units away from the line. Dividing by the slope converts it into concentration units, which is the number your clients actually care about.

The number a client receives is not the absorbance but the concentration read backwards through the line, and it inherits the uncertainty of both the line and the new measurement. If the unknown is measured $m$ times with mean response $y_0$, the standard deviation of its predicted concentration ${\hat{x}}_0$ is, to a good approximation [8]:

The three terms under the root have a plain reading: the noise of the new measurement (smaller with more replicates), the uncertainty of the line's height (smaller with more standards), and the uncertainty of its slope, which grows as $y_0$ moves away from the centre of the calibration. The line is best known at $(\bar{x},\bar{y})$ and less and less well towards the ends of the range, which is one more argument against extrapolation.

For the iron calibration in the examples below ($s_{y/x}=0.0026$ AU) and a single measurement of $y_0=0.370$, this gives $s_{x_0}=0.014$ mg/L. Multiplying by Student's $t$ for $n-2=4$ degrees of freedom (2.78 at 95 %) gives the confidence interval: $c=1.85\pm 0.04$ mg/L. A related figure is the limit of detection (LOD), the smallest concentration that can be told apart from a blank; a common estimate from the calibration is $3.3s_{y/x}/{\hat{\beta }}_1$, here about 0.044 mg/L.

Numbers alone are not enough. The single most informative diagnostic is the residual plot: the residuals $e_i=y_i-{\hat{y}}_i$ plotted against the fitted values (or against $x$). If the model is adequate, the residuals look like random noise: no trend, no curvature, no funnel.

The demo below shows five calibrations. All of them have $R^2$ above 0.99, and all of them would pass a quick glance at the line. Only the residual plot tells them apart, and for the last one only if the residuals are plotted in the order the standards were measured.

This is the same lesson as Anscombe's quartet [7], four small datasets with practically identical means, variances, correlation and regression line, but completely different shapes. Always plot your data and your residuals.

What the patterns usually mean:

• A curve (a smile or a frown) means the relationship is not linear over the range. Narrow the calibration range, add a quadratic term, or use a nonlinear model.
• A funnel means the noise grows with the signal (heteroscedasticity). The line can still be fine, but its uncertainties are wrong. Weighted least squares is the usual fix.
• Runs of positive and negative residuals in measurement order suggest drift or a missing variable.
• One isolated point far from zero flags a possible outlier. Investigate it before deleting it.

Simple regression handles a single predictor. But the same model can be written for any number of predictors, and the compact way to do it is with matrices. Stack all the measurements into a vector and the concentrations, together with a column of ones for the intercept, into a matrix:

$\mathbf{X}$ is called the design matrix: one row per sample, one column per term of the model. The least squares solution is then the familiar normal-equation formula:

This single expression reproduces the two summation formulas above, and it works unchanged when $\mathbf{X}$ has more columns. That is the door to multiple linear regression, and from there to ridge, PCR and PLS. The derivation, and the numerical cautions that come with it, are in the normal equations entry.

Nonlinearity. The Beer-Lambert law is linear only over a limited range. At high concentrations, chemical interactions, stray light or detector saturation bend the curve. Check the residuals, and restrict the range or change the model.

Outliers. One wrongly diluted standard can drag the line noticeably, especially with few standards, as you saw in the first demo. Identify it, find out what happened, and only then decide (see outliers). Robust regression methods are an option when outliers are frequent.

Heteroscedasticity. In many techniques the noise grows with the signal. The fitted line is then still unbiased, but the standard errors and confidence intervals from the usual formulas are not trustworthy. Weighted least squares, which minimises ${\sum }_i{w}_i(y_i-{\beta }_0-{\beta }_1{x}_i{)}^2$ with weights $w_i=1/s_i^2$ (where $s_i$ is the standard deviation of the response at that standard), is the standard remedy (see the homoscedasticity discussion).

Extrapolation. A calibration is only supported inside the range of the standards. A curve built from 0 to 4 mg/L says nothing reliable about a sample at 20 mg/L: dilute the sample instead.

Forcing through zero. It is tempting to drop the intercept because "a blank should read zero". Unless you have checked that the intercept is statistically indistinguishable from zero (its confidence interval, ${\hat{\beta }}_0\pm t{s}_{{\beta }_0}$, contains zero; for the iron data it is $0.0026\pm 0.0047$ AU), forcing the line through the origin can bias every prediction.

Here are complete, working implementations of simple linear regression for the iron calibration in three languages.

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

def linear_regression(x, y):
    """
    Fit y = beta_0 + beta_1 * x by ordinary least squares.

    Returns the intercept, slope, R-squared, RMSE and the residual
    standard deviation s_yx (which divides by n - 2).
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    n = len(x)

    x_mean, y_mean = x.mean(), y.mean()

    # Slope and intercept from the normal equations
    beta_1 = np.sum((x - x_mean) * (y - y_mean)) / np.sum((x - x_mean) ** 2)
    beta_0 = y_mean - beta_1 * x_mean

    # Predictions and residuals
    y_pred = beta_0 + beta_1 * x
    residuals = y - y_pred

    # Quality metrics
    ss_res = np.sum(residuals ** 2)
    ss_tot = np.sum((y - y_mean) ** 2)
    r_squared = 1 - ss_res / ss_tot
    rmse = np.sqrt(ss_res / n)
    s_yx = np.sqrt(ss_res / (n - 2))

    return beta_0, beta_1, r_squared, rmse, s_yx


# Iron calibration (Fe-phenanthroline, 510 nm)
conc = np.array([0.0, 0.5, 1.0, 2.0, 3.0, 4.0])              # mg/L
absorbance = np.array([0.004, 0.101, 0.199, 0.402, 0.594, 0.797])

b0, b1, r2, rmse, s_yx = linear_regression(conc, absorbance)
print(f"Intercept: {b0:.4f} AU")
print(f"Slope:     {b1:.4f} AU per mg/L")
print(f"R-squared: {r2:.5f}")
print(f"RMSE:      {rmse:.4f} AU   (s_yx = {s_yx:.4f} AU)")

# Predict an unknown sample
unknown_abs = 0.370
unknown_conc = (unknown_abs - b0) / b1
print(f"Unknown: A = {unknown_abs} -> c = {unknown_conc:.2f} mg/L")

# Uncertainty of the predicted concentration (m = 1 measurement of the unknown)
n, m = len(conc), 1
s_xx = np.sum((conc - conc.mean()) ** 2)
s_x0 = s_yx / b1 * np.sqrt(1 / m + 1 / n + (unknown_abs - absorbance.mean()) ** 2 / (b1 ** 2 * s_xx))
t = stats.t.ppf(0.975, n - 2)                                  # 2.776 for 4 degrees of freedom
print(f"95% interval: c = {unknown_conc:.2f} +/- {t * s_x0:.2f} mg/L")
print(f"LOD (3.3 s_yx / slope): {3.3 * s_yx / b1:.3f} mg/L")

# Calibration plot and residual plot
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

x_line = np.linspace(0, 4.2, 100)
axes[0].scatter(conc, absorbance, s=80, zorder=5, label="Standards")
axes[0].plot(x_line, b0 + b1 * x_line, "r-", linewidth=2, label=f"Fit (R² = {r2:.5f})")
axes[0].scatter(unknown_conc, unknown_abs, s=120, marker="*", color="green",
                zorder=5, label=f"Unknown ({unknown_conc:.2f} mg/L)")
axes[0].set_xlabel("Concentration (mg/L)")
axes[0].set_ylabel("Absorbance (AU)")
axes[0].set_title("Calibration")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

y_pred = b0 + b1 * conc
axes[1].scatter(y_pred, absorbance - y_pred, s=80, zorder=5)
axes[1].axhline(0, color="red", linestyle="--", alpha=0.7)
axes[1].set_xlabel("Fitted values (AU)")
axes[1].set_ylabel("Residuals (AU)")
axes[1].set_title("Residual plot")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

% Iron calibration (Fe-phenanthroline, 510 nm)
conc = [0.0; 0.5; 1.0; 2.0; 3.0; 4.0];                 % mg/L
absorbance = [0.004; 0.101; 0.199; 0.402; 0.594; 0.797];

[b0, b1, r2, rmse_val, s_yx] = linear_regression(conc, absorbance);

fprintf('Intercept: %.4f AU\n', b0);
fprintf('Slope:     %.4f AU per mg/L\n', b1);
fprintf('R-squared: %.5f\n', r2);
fprintf('RMSE:      %.4f AU   (s_yx = %.4f AU)\n', rmse_val, s_yx);

% Predict an unknown sample
unknown_abs = 0.370;
unknown_conc = (unknown_abs - b0) / b1;
fprintf('Unknown: A = %.3f -> c = %.2f mg/L\n', unknown_abs, unknown_conc);

% Uncertainty of the predicted concentration (m = 1 measurement of the unknown)
n = numel(conc); m = 1;
s_xx = sum((conc - mean(conc)).^2);
s_x0 = s_yx / b1 * sqrt(1/m + 1/n + (unknown_abs - mean(absorbance))^2 / (b1^2 * s_xx));
t = 2.776;   % Student's t, 95%, n - 2 = 4 degrees of freedom (tinv(0.975, n - 2))
fprintf('95%% interval: c = %.2f +/- %.2f mg/L\n', unknown_conc, t * s_x0);
fprintf('LOD (3.3 s_yx / slope): %.3f mg/L\n', 3.3 * s_yx / b1);

% Calibration plot and residual plot
figure('Position', [100 100 1000 400]);

subplot(1, 2, 1);
x_line = linspace(0, 4.2, 100);
scatter(conc, absorbance, 80, 'b', 'filled'); hold on;
plot(x_line, b0 + b1 * x_line, 'r-', 'LineWidth', 2);
scatter(unknown_conc, unknown_abs, 120, 'g', 'p', 'filled');
xlabel('Concentration (mg/L)'); ylabel('Absorbance (AU)');
title(sprintf('Calibration (R^2 = %.5f)', r2));
legend('Standards', 'Fit', 'Unknown', 'Location', 'northwest');
grid on;

subplot(1, 2, 2);
y_pred = b0 + b1 * conc;
scatter(y_pred, absorbance - y_pred, 80, 'b', 'filled');
yline(0, '--r', 'LineWidth', 1.5);
xlabel('Fitted values (AU)'); ylabel('Residuals (AU)');
title('Residual plot');
grid on;

% Local functions must be placed at the end of a script file
function [beta0, beta1, r_squared, rmse, s_yx] = linear_regression(x, y)
    x = x(:); y = y(:);
    n = numel(x);
    beta1 = sum((x - mean(x)) .* (y - mean(y))) / sum((x - mean(x)).^2);
    beta0 = mean(y) - beta1 * mean(x);
    residuals = y - (beta0 + beta1 * x);
    ss_res = sum(residuals.^2);
    ss_tot = sum((y - mean(y)).^2);
    r_squared = 1 - ss_res / ss_tot;
    rmse = sqrt(ss_res / n);
    s_yx = sqrt(ss_res / (n - 2));
end

linear_regression <- function(x, y) {
  # Fit y = beta0 + beta1 * x by ordinary least squares
  n <- length(x)
  beta1 <- sum((x - mean(x)) * (y - mean(y))) / sum((x - mean(x))^2)
  beta0 <- mean(y) - beta1 * mean(x)

  y_pred <- beta0 + beta1 * x
  residuals <- y - y_pred
  ss_res <- sum(residuals^2)
  ss_tot <- sum((y - mean(y))^2)

  list(beta0 = beta0, beta1 = beta1,
       r_squared = 1 - ss_res / ss_tot,
       rmse = sqrt(ss_res / n),
       s_yx = sqrt(ss_res / (n - 2)),
       fitted = y_pred, residuals = residuals)
}

# Iron calibration (Fe-phenanthroline, 510 nm)
conc <- c(0.0, 0.5, 1.0, 2.0, 3.0, 4.0)                 # mg/L
absorbance <- c(0.004, 0.101, 0.199, 0.402, 0.594, 0.797)

fit <- linear_regression(conc, absorbance)
cat(sprintf("Intercept: %.4f AU\n", fit$beta0))
cat(sprintf("Slope:     %.4f AU per mg/L\n", fit$beta1))
cat(sprintf("R-squared: %.5f\n", fit$r_squared))
cat(sprintf("RMSE:      %.4f AU   (s_yx = %.4f AU)\n", fit$rmse, fit$s_yx))

# Predict an unknown sample
unknown_abs <- 0.370
unknown_conc <- (unknown_abs - fit$beta0) / fit$beta1
cat(sprintf("Unknown: A = %.3f -> c = %.2f mg/L\n", unknown_abs, unknown_conc))

# Uncertainty of the predicted concentration (m = 1 measurement of the unknown)
n <- length(conc); m <- 1
s_xx <- sum((conc - mean(conc))^2)
s_x0 <- fit$s_yx / fit$beta1 *
  sqrt(1 / m + 1 / n + (unknown_abs - mean(absorbance))^2 / (fit$beta1^2 * s_xx))
t <- qt(0.975, df = n - 2)                              # 2.776 for 4 degrees of freedom
cat(sprintf("95%% interval: c = %.2f +/- %.2f mg/L\n", unknown_conc, t * s_x0))
cat(sprintf("LOD (3.3 s_yx / slope): %.3f mg/L\n", 3.3 * fit$s_yx / fit$beta1))

# Calibration plot and residual plot
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

plot(conc, absorbance, pch = 16, cex = 1.5, col = "steelblue",
     xlab = "Concentration (mg/L)", ylab = "Absorbance (AU)",
     main = sprintf("Calibration (R² = %.5f)", fit$r_squared))
abline(a = fit$beta0, b = fit$beta1, col = "red", lwd = 2)
points(unknown_conc, unknown_abs, pch = 8, cex = 2, col = "forestgreen")
grid()

plot(fit$fitted, fit$residuals, pch = 16, cex = 1.5, col = "steelblue",
     xlab = "Fitted values (AU)", ylab = "Residuals (AU)",
     main = "Residual plot")
abline(h = 0, col = "red", lty = 2, lwd = 2)
grid()

Linear regression is a straight line chosen by a precise criterion: the smallest possible sum of squared vertical residuals. Its formulas come straight out of the normal equations, and its justification comes from the assumptions about measurement noise that we met in the least squares series.

Three habits will save you from most of its traps: plot the data, plot the residuals, and stay inside the calibrated range. A high $R^2$ is a necessary courtesy, never a proof.

Everything else in this section builds on this one idea. Multiple linear regression adds more predictors. Ridge and lasso regression add a penalty to keep the coefficients sensible when predictors overlap. Principal component regression and partial least squares compress whole spectra into a few latent variables before regressing [10]. For a full treatment of regression diagnostics, Draper and Smith remains a classic [9].