10 - Simple and Multiple Linear Regression

You have measured the response time of a service as the load varies and you want to predict how long it will take at a load you have not tried. You draw a line through the points, add a 95% band of uncertainty, and that band is wrong: if you use the confidence interval that your software calls “95%”, a future measurement falls inside it only 44.5% of the time, as the simulation below shows. There are two different intervals, one for the mean response and one for a single future measurement, and confusing them is the most frequent mistake in regression.

This article, the tenth in the Probability and Statistics series, builds the least-squares line starting from five points, tells you what its assumptions guarantee and where they break. You need the theory of estimators and the t-test from the previous articles; from the chi-square you carry over the idea of comparing what you observe with what the model predicts, which here is applied to the residuals.

The least-squares line

The simple linear regression model says that the response y_i is a line plus an error: y_i = β₀ + β₁x_i + ε_i, with independent errors ε_i, of mean zero and constant variance σ². The parameters β₀ (intercept) and β₁ (slope: by how much y changes on average for one unit of x) are unknown.

The least squares method chooses the line that minimizes the sum of the squared residuals, SSE = Σ(y_i − b₀ − b₁x_i)², where the residual is the vertical distance of each point from the line. Setting the derivatives to zero gives b₁ = S_xy/S_xx with S_xy = Σ(x_i − x̄)(y_i − ȳ) and S_xx = Σ(x_i − x̄)², and b₀ = ȳ − b₁x̄: the line always passes through the point of the means.

The example is deliberately small, so that it can be redone by hand: x = 1, 2, 3, 4, 5 and y = 2.1, 3.9, 6.2, 7.8, 10.1. The means are 3 and 6.02, S_xx = 10 and S_xy = 19.90, so b₁ = 1.99 and b₀ = 6.02 − 1.99·3 = 0.05. The estimate of the error variance is s² = SSE/(n − 2): you divide by n − 2 and not by n because two parameters were estimated from the data, as for the sample variance in the article on estimators.

The standard error of the slope is s/√S_xx. The program repeats the calculations, writes the same problem in matrix form (normal equations X′Xb = X′y, also solved with the QR decomposition and with np.linalg.lstsq) and compares with scipy.stats.linregress and with statsmodels.

# Simple linear regression: worked exercise with 5 points, hand formulas, matrix form and check with scipy and statsmodels
import numpy as np
from scipy import stats
import statsmodels.api as sm

x = np.array([1, 2, 3, 4, 5.0])
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1])
n = len(x)

Sxx = ((x - x.mean()) ** 2).sum()
Sxy = ((x - x.mean()) * (y - y.mean())).sum()
b1 = Sxy / Sxx                               # slope
b0 = y.mean() - b1 * x.mean()                # intercept: the line passes through (mean x, mean y)
e = y - (b0 + b1 * x)                        # residuals
sse = (e ** 2).sum()
sst = ((y - y.mean()) ** 2).sum()
s2 = sse / (n - 2)                           # estimate of sigma^2: n - 2 degrees of freedom (two estimated parameters)
se1 = np.sqrt(s2 / Sxx)
t1 = b1 / se1
p1 = 2 * stats.t.sf(abs(t1), n - 2)
q = stats.t.ppf(0.975, n - 2)
print(f"Sxx = {Sxx:.1f}, Sxy = {Sxy:.2f}, slope = {b1:.3f}, intercept = {b0:.3f}")
print(f"residuals = {np.round(e, 3).tolist()}, sum of residuals zero: {abs(e.sum()) < 1e-12}, SSE = {sse:.4f}, SST = {sst:.4f}")
print(f"R^2 = {1 - sse / sst:.5f}, estimated sigma^2 = {s2:.5f}")
print(f"standard error of the slope = {se1:.4f}, t = {t1:.2f} with {n - 2} d.f., p = {p1:.2e}")
print(f"95% CI of the slope = [{b1 - q * se1:.3f}, {b1 + q * se1:.3f}]")

# Matrix form: normal equations X'X b = X'y, QR decomposition, lstsq
X = np.column_stack([np.ones(n), x])
b_normali = np.linalg.solve(X.T @ X, X.T @ y)
Q, R = np.linalg.qr(X)
b_qr = np.linalg.solve(R, Q.T @ y)
b_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"coefficients: normal {np.round(b_normali, 6).tolist()}, QR {np.round(b_qr, 6).tolist()}, lstsq {np.round(b_lstsq, 6).tolist()}")

r = stats.linregress(x, y)
m = sm.OLS(y, X).fit()
print(f"scipy.linregress: slope {r.slope:.3f}, standard error {r.stderr:.4f}, p = {r.pvalue:.2e}, R^2 = {r.rvalue ** 2:.5f}")
print(f"statsmodels:      slope {m.params[1]:.3f}, standard error {m.bse[1]:.4f}, p = {m.pvalues[1]:.2e}, R^2 = {m.rsquared:.5f}")

# Output:
# Sxx = 10.0, Sxy = 19.90, slope = 1.990, intercept = 0.050
# residuals = [0.06, -0.13, 0.18, -0.21, 0.1], sum of residuals zero: True, SSE = 0.1070, SST = 39.7080
# R^2 = 0.99731, estimated sigma^2 = 0.03567
# standard error of the slope = 0.0597, t = 33.32 with 3 d.f., p = 5.94e-05
# 95% CI of the slope = [1.800, 2.180]
# coefficients: normal [0.05, 1.99], QR [0.05, 1.99], lstsq [0.05, 1.99]
# scipy.linregress: slope 1.990, standard error 0.0597, p = 5.94e-05, R^2 = 0.99731
# statsmodels:      slope 1.990, standard error 0.0597, p = 5.94e-05, R^2 = 0.99731

The three matrix methods give the same coefficients (0.05 and 1.99), and SciPy and statsmodels also coincide with the hand calculations: standard error 0.0597, t = 33.32 with 3 degrees of freedom, p = 5.94 · 10⁻⁵, 95% interval [1.800; 2.180]. The residuals sum to zero, as they must when there is an intercept. The index R² = 1 − SSE/SST = 0.99731 says that the line explains 99.7% of the total variability of the y.

For the normal equations it is better to use the QR decomposition and not to invert X′X: the condition number of X′X is the square of that of X, and with nearly collinear variables digits are lost. The formulas for the general case with p coefficients, and the orthogonal projection property, are in the regression lecture of the MIT 18.650 course.

What the assumptions guarantee

With errors of mean zero, constant variance and uncorrelated, the least-squares estimators are the best among the linear unbiased ones (Gauss-Markov theorem). To do exact inference you need more: Normal errors.

Then the estimates are Normal, (n − p)s²/σ² follows a chi-square with n − p degrees of freedom, and is independent of the estimates; hence the t test on a single coefficient with n − p degrees of freedom, stated in the same MIT lecture. For the simple line p = 2. To test all the coefficients together the F test is used (which for a single explanatory variable is the square of the t).

The practical value lies in the intervals. For a point x₀ three 95% intervals can be built, all with the t quantile at n − 2 degrees of freedom. The interval for the slope is b₁ ± t·s/√S_xx. The interval for the mean response at x₀ is ŷ₀ ± t·s√h₀, with h₀ = 1/n + (x₀ − x̄)²/S_xx, the leverage of the point: it grows as you move away from the mean of the x. The prediction interval for a new measurement is ŷ₀ ± t·s√(1 + h₀): the extra “1” is the error ε of the new measurement, which no amount of data eliminates.

The program simulates 100,000 experiments with 30 fixed points between 0 and 10, true slope 2, intercept 1 and Normal errors with standard deviation 1.5, and for each one builds the three intervals at x₀ = 9 with a new future measurement. It counts how many times each interval contains its target. It adds a fourth line, the error described in the opening: the interval for the mean used as if it were a prediction interval. At the end it compares with get_prediction from statsmodels.

# Confidence interval for the mean, prediction interval and interval for the slope: actual coverage
import numpy as np
from scipy import stats
import statsmodels.api as sm

n, R = 30, 100_000
beta0, beta1, sigma = 1.0, 2.0, 1.5
rng = np.random.default_rng(11)
x = np.linspace(0, 10, n)                             # fixed design
x0 = 9.0                                              # point at which we want to predict
Sxx = ((x - x.mean()) ** 2).sum()
q = stats.t.ppf(0.975, n - 2)

Y = beta0 + beta1 * x + sigma * rng.standard_normal((R, n))
xm = x.mean()
b1 = ((x - xm) * (Y - Y.mean(axis=1, keepdims=True))).sum(axis=1) / Sxx
b0 = Y.mean(axis=1) - b1 * xm
e = Y - b0[:, None] - b1[:, None] * x
s = np.sqrt((e ** 2).sum(axis=1) / (n - 2))
yhat0 = b0 + b1 * x0
h0 = 1 / n + (x0 - xm) ** 2 / Sxx                     # leverage of the point x0
se_media = s * np.sqrt(h0)
se_prev = s * np.sqrt(1 + h0)
media_vera = beta0 + beta1 * x0
y_nuovo = media_vera + sigma * rng.standard_normal(R)  # one future observation for each experiment

print(f"n = {n}, x0 = {x0}, experiments = {R}; standard error of a coverage: {np.sqrt(0.95 * 0.05 / R):.4f}")
print(f"CI of the slope covers the true value:        {(abs(b1 - beta1) <= q * s / np.sqrt(Sxx)).mean():.4f}")
print(f"CI of the mean at x0 covers the true mean:    {(abs(yhat0 - media_vera) <= q * se_media).mean():.4f}")
print(f"CI of the mean used as a prediction:          {(abs(yhat0 - y_nuovo) <= q * se_media).mean():.4f}")
print(f"prediction interval covers the new value:     {(abs(yhat0 - y_nuovo) <= q * se_prev).mean():.4f}")
print(f"ratio of the half-widths, prediction/mean: exact {np.sqrt((1 + h0) / h0):.3f}, simulated {(se_prev / se_media).mean():.3f}")

# Check with statsmodels on one sample
m = sm.OLS(Y[0], sm.add_constant(x)).fit()
pr = m.get_prediction(np.array([[1.0, x0]])).summary_frame(alpha=0.05)
print(f"statsmodels check (first experiment): mean at x0 {pr['mean'].iloc[0]:.4f} versus {yhat0[0]:.4f}; "
      f"CI half-width {(pr['mean_ci_upper'].iloc[0] - pr['mean'].iloc[0]):.4f} versus {q * se_media[0]:.4f}; "
      f"prediction half-width {(pr['obs_ci_upper'].iloc[0] - pr['mean'].iloc[0]):.4f} versus {q * se_prev[0]:.4f}")

# Output:
# n = 30, x0 = 9.0, experiments = 100000; standard error of a coverage: 0.0007
# CI of the slope covers the true value:        0.9504
# CI of the mean at x0 covers the true mean:    0.9506
# CI of the mean used as a prediction:          0.4453
# prediction interval covers the new value:     0.9497
# ratio of the half-widths, prediction/mean: exact 3.425, simulated 3.425
# statsmodels check (first experiment): mean at x0 17.9370 versus 17.9370; CI half-width 0.7687 versus 0.7687; prediction half-width 2.6326 versus 2.6326

The three correct intervals cover as promised: 95.04% for the slope, 95.06% for the mean, 94.97% for the prediction (standard error 0.0007). The interval for the mean, used as a prediction, covers only 44.53% of future measurements: it is too narrow, because it ignores the variability of the single measurement. The half-widths stand in the ratio √((1 + h₀)/h₀) = 3.425, both by the formula and in the simulation: the prediction interval is more than three times wider. The comparison with statsmodels at the end (half-widths 0.7687 and 2.6326) gives the same values, so the implemented formula is that of the library.

More variables: three things that change

In multiple regression the response depends on p − 1 variables: y = β₀ + β₁x₁ + … + ε, and each coefficient β_j is the average change in y for one unit of x_j holding the other variables fixed. Estimation, tests and intervals are the same, with n − p degrees of freedom. Three phenomena arise only with more variables, and the program measures them.

First: adding variables raises R² even if they are pure chance. With n = 30 and 10 variables independent of y, the expected value of R² is p/(n − 1) = 10/29, an exact result (the law of R² is a Beta). The remedy is the adjusted R², 1 − (1 − R²)(n − 1)/(n − p − 1), which compensates for the degrees of freedom lost.

Second: omitting a variable correlated with the included ones shifts the coefficient. If the true model has x₁ and x₂ (coefficients 2 and 3, correlation 0.6 between the two standardized variables) and you estimate only x₁, the expected slope is not 2 but 2 + 3 · 0.6 = 3.8: it absorbs the effect of the omitted variable. Third: collinearity. When two explanatory variables are strongly correlated (correlation ρ), the standard deviation of each coefficient grows by a factor √VIF, with VIF = 1/(1 − ρ²) (variance inflation factor): the model still predicts well, but it cannot attribute the effect to one variable or the other.

# Multiple regression: R^2 with pure-noise variables, omitted-variable bias, collinearity
import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

def ols(X, y):
    """Coefficients via QR; also returns the classical standard errors."""
    Q, R = np.linalg.qr(X)
    b = np.linalg.solve(R, Q.T @ y)
    e = y - X @ b
    s2 = e @ e / (len(y) - X.shape[1])
    Rinv = np.linalg.inv(R)
    return b, np.sqrt(s2 * np.diag(Rinv @ Rinv.T)), e

# 1) Adding p pure-noise variables raises R^2 even if there is no relationship
n, p, Rep = 30, 10, 20_000
r2 = np.empty(Rep); r2_agg = np.empty(Rep)
for i in range(Rep):
    X = np.column_stack([np.ones(n), rng.standard_normal((n, p))])
    y = rng.standard_normal(n)                              # y does not depend on any X
    _, _, e = ols(X, y)
    r2[i] = 1 - e @ e / ((y - y.mean()) @ (y - y.mean()))
    r2_agg[i] = 1 - (1 - r2[i]) * (n - 1) / (n - p - 1)       # R^2 adjusted for degrees of freedom
print(f"1) n = {n}, p = {p} noise variables: simulated mean R^2 {r2.mean():.4f}, exact p/(n-1) = {p / (n - 1):.4f}; "
      f"mean adjusted R^2 {r2_agg.mean():.4f}")
print(f"   R^2 exceeds 0.5 in {(r2 > 0.5).mean():.4f} of the cases (exact, Beta law: {stats.beta.sf(0.5, p / 2, (n - p - 1) / 2):.4f})")

# 2) Omitted variable: y = 1 + 2 x1 + 3 x2 + e, with corr(x1, x2) = 0.6; y is regressed on x1 alone
n, rho, Rep = 200, 0.6, 20_000
b_corto = np.empty(Rep); b_pieno = np.empty(Rep)
for i in range(Rep):
    z = rng.standard_normal((n, 2))
    x1 = z[:, 0]; x2 = rho * z[:, 0] + np.sqrt(1 - rho ** 2) * z[:, 1]
    y = 1 + 2 * x1 + 3 * x2 + rng.standard_normal(n)
    b_corto[i] = ols(np.column_stack([np.ones(n), x1]), y)[0][1]
    b_pieno[i] = ols(np.column_stack([np.ones(n), x1, x2]), y)[0][1]
print(f"2) slope of x1 (true 2): without x2 mean {b_corto.mean():.4f} (expected 2 + 3*0.6 = {2 + 3 * rho:.4f}); with x2 mean {b_pieno.mean():.4f}")

# 3) Collinearity: standard deviation of the slope of x1 as the correlation with x2 grows
print("3) corr(x1, x2)  VIF = 1/(1-r^2)  simulated std. dev.  expected 1/sqrt(n(1-r^2))")
for rho in (0.0, 0.9, 0.99):
    b = np.empty(5_000)
    for i in range(b.size):
        z = rng.standard_normal((n, 2))
        x1 = z[:, 0]; x2 = rho * z[:, 0] + np.sqrt(1 - rho ** 2) * z[:, 1]
        y = 1 + 2 * x1 + 3 * x2 + rng.standard_normal(n)
        b[i] = ols(np.column_stack([np.ones(n), x1, x2]), y)[0][1]
    print(f"   {rho:13.2f}  {1 / (1 - rho ** 2):12.2f}  {b.std():19.4f}  {1 / np.sqrt(n * (1 - rho ** 2)):.4f}")

# Output:
# 1) n = 30, p = 10 noise variables: simulated mean R^2 0.3442, exact p/(n-1) = 0.3448; mean adjusted R^2 -0.0010
#    R^2 exceeds 0.5 in 0.1080 of the cases (exact, Beta law: 0.1097)
# 2) slope of x1 (true 2): without x2 mean 3.8013 (expected 2 + 3*0.6 = 3.8000); with x2 mean 1.9992
# 3) corr(x1, x2)  VIF = 1/(1-r^2)  simulated std. dev.  expected 1/sqrt(n(1-r^2))
#             0.00          1.00               0.0717  0.0707
#             0.90          5.26               0.1656  0.1622
#             0.99         50.25               0.5109  0.5013

With 10 noise variables, R² is on average 0.3442 (exact 0.3448) and exceeds 0.5 in 10.8% of cases (exact 0.1097): a model with no relationship at all to the response appears to explain a third of the variance. The adjusted R² has mean −0.0010, that is, zero. Omitting x₂, the slope of x₁ has mean 3.8013, against 3.8 expected; including it, 1.9992, close to the true 2.

With correlation 0.90 the VIF is 5.26 and the standard deviation of the coefficient goes from 0.0717 to 0.1656; with 0.99 the VIF is 50.25 and the standard deviation is 0.5109, 7 times larger. A non-significant coefficient in the presence of strong collinearity does not mean that the variable does not matter: it means that the data are not enough to separate it from the other.

When the error variance is not constant

The second assumption, constant variance (homoscedasticity), is the one that breaks most often: the error of a measurement of time, price or income tends to grow with the level. The estimated coefficients remain unbiased, but the classical standard error s√(X′X)⁻¹_jj is wrong, and with it intervals and p-values. White’s heteroscedasticity-robust standard errors, in the HC3 variant, estimate the variance with (X′X)⁻¹ [Σ e_i²/(1 − h_i)² x_ix_i′] (X′X)⁻¹, where h_i is the leverage of observation i; statsmodels offers them with cov_type="HC3", described in the documentation of RegressionResults.

The program simulates 20,000 experiments with 100 points in [0, 1], true slope 2, and an error standard deviation that grows from 0.2 to 3.2 (σ_i = 0.2 + 3x_i²). With a fixed design the true standard deviation of the slope is computed exactly with the “sandwich” formula with the true variances; the program compares the formula with the simulation, then measures how much the 95% interval built with the classical standard error and with HC3 covers. It checks its own implementation of HC3 against statsmodels.

# When the error variance grows with x: classical standard errors versus robust HC3 errors
import numpy as np
import statsmodels.api as sm
from scipy import stats

n, R = 100, 20_000
rng = np.random.default_rng(99)
x = np.linspace(0, 1, n)                              # fixed design
X = np.column_stack([np.ones(n), x])
sigma_i = 0.2 + 3.0 * x ** 2                          # standard deviation of the error: from 0.2 to 3.2
XtXi = np.linalg.inv(X.T @ X)
H = X @ XtXi @ X.T
h = np.diag(H)

# True standard deviation of the slope estimator (exact formula, design fixed)
cov_vera = XtXi @ (X.T * sigma_i ** 2) @ X @ XtXi
sd_vera = np.sqrt(cov_vera[1, 1])

Y = 1.0 + 2.0 * x + sigma_i * rng.standard_normal((R, n))
B = Y @ X @ XtXi                                      # coefficients of all experiments, shape (R, 2)
E = Y - B @ X.T
s2 = (E ** 2).sum(axis=1) / (n - 2)
se_classico = np.sqrt(s2 * XtXi[1, 1])
# HC3: (X'X)^-1 [sum e_i^2 / (1-h_i)^2 x_i x_i'] (X'X)^-1
w = E ** 2 / (1 - h) ** 2
meat = np.einsum('ri,ij,ik->rjk', w, X, X)
se_hc3 = np.sqrt(np.einsum('a,rab,b->r', XtXi[1], meat, XtXi[:, 1]))

q = stats.t.ppf(0.975, n - 2)
cop = lambda se: (abs(B[:, 1] - 2.0) <= q * se).mean()
print(f"true std. dev. of the slope (formula) = {sd_vera:.4f}; simulated = {B[:, 1].std():.4f}")
print(f"classical standard error: mean {se_classico.mean():.4f}  -> 95% CI covers {cop(se_classico):.4f}")
print(f"HC3 standard error:       mean {se_hc3.mean():.4f}  -> 95% CI covers {cop(se_hc3):.4f}")

# Check of HC3 with statsmodels on the first experiment
m = sm.OLS(Y[0], X).fit(cov_type="HC3")
print(f"statsmodels check (first experiment): HC3 {m.bse[1]:.6f} versus {se_hc3[0]:.6f}")

# Output:
# true std. dev. of the slope (formula) = 0.6326; simulated = 0.6378
# classical standard error: mean 0.5116  -> 95% CI covers 0.8852
# HC3 standard error:       mean 0.6316  -> 95% CI covers 0.9432
# statsmodels check (first experiment): HC3 0.469747 versus 0.469747

The true standard deviation of the slope is 0.6326 (simulated 0.6378). The classical standard error is on average 0.5116, that is, it underestimates by 19%, and the resulting interval covers 88.52% of the time instead of 95%. HC3 is on average 0.6316, close to the true value, and the interval covers 94.32%. The check with statsmodels coincides to six digits (0.469747) on the first experiment. With a line and 100 points the cost of HC3 is nil: if you have doubts about homoscedasticity, using the robust error does not go wrong when the variance is constant and protects when it is not, as in this case.

Four mistakes that are really made

• Presenting the confidence interval of the mean as a prediction band: it covers future measurements 44.5% of the time.
• Adding variables until R² goes up: with 10 noise variables on 30 points it goes up to 0.34 on average.
• Reading a coefficient without asking which correlated variables are missing: the slope of 2 becomes 3.8.
• Trusting the classical standard errors when the error grows with x: the “95%” interval covers 88.5%.

The four mistakes are mistakes of reading: the model does what it must, and the trouble is that the printed number (95%, R², coefficient) is read as if it held under conditions that do not hold. For each one you have a measured number that says by how much: if your problem has a plausible omitted variable, a growing variance or two strongly correlated variables, you already have an idea of the order of magnitude of the error.

You can now build by hand and with QR the least-squares line, compute standard errors, t tests and R², distinguish the interval for the mean from the prediction interval and verify their coverage with a simulation, recognise the effects of useless, omitted and collinear variables, and switch to HC3 robust errors when the variance is not constant. The linear model is also the destination of ANOVA, which is a regression with group variables, and in the article on Bayesian estimation you will see how a prior distribution changes the way a parameter is estimated.

The same series also includes Exploratory Data Analysis (EDA) with Python.