Linear regression

Much of engineering is about relationships between variables: concrete strength and water–cement ratio, sensor voltage and applied load, fuel consumption and speed, pavement deterioration and traffic. Regression finds the equation that best describes such a relationship from data, measures how well it fits, and, crucially, says how uncertain predictions made from it are. It is the statistical tool engineers use most, from calibrating instruments to building the empirical equations found throughout design codes.

This chapter covers scatter plots and correlation, least-squares fitting, residuals and R², inference about the slope, confidence and prediction intervals, and the traps: extrapolation, outliers, and mistaking correlation for causation.

Scatter plots and correlation

Always start with a scatter plot. It shows the form of the relationship (straight or curved), its direction, its strength, and any unusual points, none of which a single number can capture. The correlation coefficient r summarizes the strength of a straight-line relationship:

r lies between −1 and +1. Values near ±1 mean the points cluster tightly about a line; values near 0 mean no linear relationship. The sign gives the direction. r measures only linear association: a perfect U-shaped relationship can have r = 0.

Fitting a line by least squares

The simple linear regression model is

where ε is random scatter about the line, assumed independent with mean 0 and constant standard deviation σ. The least-squares line chooses the intercept \(b_0\) and slope \(b_1\) that minimize the sum of squared vertical distances (residuals) from the points to the line. Setting the partial derivatives to zero, the same optimization as MATH-209, gives

The line always passes through the point \((\bar{x}, \bar{y})\). The slope is the estimated change in y per unit change in x.

Residuals and goodness of fit

A residual is the vertical distance from a point to the line, \(e_i = y_i - \hat{y}_i\). The residual standard deviation estimates σ, the typical scatter about the line:

Two parameters were estimated, so the divisor is n − 2. The coefficient of determination \(R^2\) is the fraction of the variation in y explained by the line:

For the concrete data, the residual sum of squares is 17.13, so \(s_e = \sqrt{17.13/8} = 1.46\) MPa and \(R^2 = 1 - 17.13/477.54 = 0.964\): the water–cement ratio explains 96% of the variation in strength among these batches, and the remaining scatter is about ±1.5 MPa.

A high R² does not by itself mean the model is right. The residual plot (residuals against x or against fitted values) is the key diagnostic. Residuals should look like random scatter in a horizontal band. A curved pattern means a straight line is the wrong form; a funnel shape means the scatter grows with x (non-constant variance); a single huge residual flags an outlier to investigate.

Inference about the slope

The slope estimate varies from sample to sample, with standard error

A confidence interval is \(b_1 \pm t_{\alpha/2,\,n-2}\,SE(b_1)\), and the test of H₀: β₁ = 0 (no linear relationship) uses \(t = b_1/SE(b_1)\) with n − 2 degrees of freedom.

Predicting: confidence and prediction intervals

At a chosen value \(x_0\), the line gives \(\hat{y}_0 = b_0 + b_1x_0\). Two different intervals answer two different questions:

• Confidence interval for the mean response: where is the average y of all batches at \(x_0\)?

• Prediction interval for one new observation: where will the next single batch fall?

The extra 1 under the root is the scatter of an individual about the mean, so prediction intervals are always wider, and they do not shrink to zero however much data you collect. Both widen as \(x_0\) moves away from \(\bar{x}\): the model is best known near the middle of the data.

Extrapolation, outliers and causation

A regression is only trustworthy within the range of the data. The concrete line predicts 60 MPa at w/c = 0.25 and 11 MPa at 0.80, but the real relationship is curved (strength falls roughly exponentially with w/c, Abrams' law), very low ratios cannot be compacted properly, and no data were taken there. Extrapolation is the most common misuse of regression in practice.

Influential points, observations far from the others in x, can pull the line strongly toward themselves. Check them: a data-entry error should be corrected; a genuine unusual point should be reported, and the model fitted with and without it.

Finally, correlation is not causation. Across Canadian cities, the number of fire trucks dispatched correlates with fire damage, but sending trucks does not cause damage; bigger fires cause both. A third variable (a lurking or confounding variable) can create a correlation between two variables that do not affect each other. Only a designed experiment, in which x is deliberately set and other factors are controlled or randomized, as in the trial batches above, supports a causal conclusion directly.

Calibration and multiple regression

Calibrating an instrument is regression in reverse. Known inputs (reference loads) are applied, the outputs (voltages) recorded and a line fitted; in use, an output reading is converted back to the input by solving the line for x. A load cell giving 0.02, 1.03, 2.05, 2.98, 4.01 and 5.03 V at 0, 10, 20, 30, 40 and 50 kN fits \(V = 0.0257 + 0.09977L\), so a reading of 3.50 V means \(L = (3.50 - 0.0257)/0.09977 = 34.8\) kN.

Real responses usually depend on several variables. Multiple regression fits \(y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \cdots\) by the same least-squares principle, solved with matrices (MATH-125): for example, strength as a function of w/c, cement content and curing temperature, or energy use of a building as a function of outdoor temperature and occupancy. Each coefficient estimates the effect of its variable with the others held fixed. Software reports each coefficient with its standard error, and the same cautions about residual plots, extrapolation and causation apply.

import numpy as np
from scipy import stats

wc = np.array([0.35, 0.38, 0.40, 0.42, 0.45, 0.48, 0.50, 0.53, 0.55, 0.60])
fc = np.array([51.2, 50.3, 46.1, 45.9, 41.0, 40.6, 36.2, 36.9, 31.8, 30.4])

fit = stats.linregress(wc, fc)
print(f"strength = {fit.intercept:.2f} + ({fit.slope:.2f}) w/c")
print(f"r = {fit.rvalue:.3f}, R^2 = {fit.rvalue**2:.3f}, SE(slope) = {fit.stderr:.2f}")

# 95% prediction interval for one batch at w/c = 0.45
n = len(wc)
resid = fc - (fit.intercept + fit.slope * wc)
se = np.sqrt(np.sum(resid**2) / (n - 2))
x0 = 0.45
y0 = fit.intercept + fit.slope * x0
half = stats.t.ppf(0.975, n - 2) * se * np.sqrt(1 + 1/n + (x0 - wc.mean())**2 / np.sum((wc - wc.mean())**2))
print(f"At w/c = {x0}: {y0:.2f} MPa, 95% PI {y0 - half:.1f} to {y0 + half:.1f} MPa")

Output: strength = 82.54 + (-89.07) w/c; r = -0.982, R^2 = 0.964, SE(slope) = 6.07; At w/c = 0.45: 42.47 MPa, 95% PI 38.9 to 46.0 MPa.