Linear Regression: Coefficients, Assumptions, and Diagnostics

A fitted linear regression gives coefficients, predictions, and often standard errors, but these answer different questions. This article focuses on interpreting a coefficient, choosing assumptions for uncertainty estimates, and checking what a diagnostic establishes. Linear Algebra for Machine Learning explains the least-squares calculation used to fit the model.

What a coefficient says

In \(y=\beta_0+\beta_1x_1+\cdots+\beta_px_p+\varepsilon\), the intercept \(\beta_0\) is the model value when all inputs are zero, and \(\varepsilon\) is the difference between y and the full linear expression. If \(E[\varepsilon\mid X]=0\), the linear expression is the conditional mean. For these untransformed inputs without interactions, \(\beta_j\) gives its change per unit of \(x_j\) with the other included variables fixed. A fitted coefficient estimates that relationship; if the mean model is misspecified, it describes a linear approximation instead. Terms such as \(x^2\) or interactions require interpreting the combined expression.

For example, a coefficient of 3 dollars per square meter adds 30 dollars to the predicted value for ten additional square meters, holding the other inputs fixed. That comparison can be poorly supported if the data contain no such combinations. Coefficients can change when other variables enter the model, but need not always do so. With exact collinearity, separate coefficients are not uniquely identified; software may return one solution, such as a minimum-norm solution. Neither a reported coefficient nor the phrase “holding fixed” establishes the effect of an intervention.

The first exercise constructs a sign reversal by omitting a variable related to both the included predictor and the outcome. It provides one explanation to investigate, not a diagnosis for every unexpected coefficient.

The assumptions, and what each one is for

Normality is not required to compute an OLS fit or to obtain unbiased coefficients under the mean assumption. It supplies exact finite-sample t and F reference distributions under the additional classical assumptions. Large-sample approximations can work without normality, but heavy tails, influential observations, dependence, or a large parameter count can make them inaccurate. A residual QQ plot compares empirical residual quantiles with normal quantiles; it is one diagnostic, not a pass/fail certificate for the analysis.

Heteroscedasticity means that error variance changes with the inputs. Dependence means errors from different rows are related. Either can invalidate the usual uncertainty formula while leaving OLS conditionally unbiased if the zero-conditional-mean and rank conditions still hold. For example, a time-series predictor related to the current error can violate the zero-conditional-mean condition as well. Distinguish the coefficient estimate from its standard error before choosing a remedy.

Heteroscedasticity, measured

In the code, Xd has an intercept column of ones and the predictor x. b holds the fitted intercept and slope; r = y - Xd @ b contains residuals, the observed minus fitted values. Residuals estimate errors but are not the unobserved errors themselves. The covariance formulas produce estimated coefficient variances; taking the square root gives a standard error. The experiment repeats the entire fit 4,000 times for each noise pattern.

import numpy as np

rng = np.random.default_rng(1)
R, n = 4000, 100

for power, label in ((0, "constant sd"), (2, "sd ~ x^2"), (3, "sd ~ x^3")):
    hit_ols = hit_robust = hit_hc3 = 0
    slopes = []
    for _ in range(R):
        x = rng.uniform(1, 5, n)
        y = 2 * x + rng.normal(size=n) * (x ** power) * 0.3
        Xd = np.column_stack([np.ones(n), x])
        b = np.linalg.lstsq(Xd, y, rcond=None)[0]
        slopes.append(b[1])
        r = y - Xd @ b
        XtXi = np.linalg.inv(Xd.T @ Xd)
        se_ols = np.sqrt((r @ r / (n - 2)) * XtXi[1, 1])
        se_hc = np.sqrt((XtXi @ (Xd.T @ (r[:, None] ** 2 * Xd)) @ XtXi)[1, 1])
        leverage = np.sum((Xd @ XtXi) * Xd, axis=1)
        hc3_weights = r ** 2 / (1 - leverage) ** 2
        se_hc3 = np.sqrt((XtXi @ (Xd.T @ (hc3_weights[:, None] * Xd)) @ XtXi)[1, 1])
        hit_hc3 += abs(b[1] - 2) <= 1.96 * se_hc3
        hit_ols += abs(b[1] - 2) <= 1.96 * se_ols
        hit_robust += abs(b[1] - 2) <= 1.96 * se_hc
    mc_se = np.std(slopes, ddof=1) / np.sqrt(R)
    print(f"{label:12s} mean slope {np.mean(slopes):.4f}"
          f"  mean MC SE {mc_se:.4f}  usual coverage {hit_ols / R:.4f}"
          f"   HC0 {hit_robust / R:.4f}   HC3 {hit_hc3 / R:.4f}")
# constant sd  mean slope 1.9994  mean MC SE 0.0004  usual coverage 0.9497   HC0 0.9463   HC3 0.9517
# sd ~ x^2     mean slope 2.0055  mean MC SE 0.0059  usual coverage 0.9042   HC0 0.9435   HC3 0.9485
# sd ~ x^3     mean slope 2.0378  mean MC SE 0.0282  usual coverage 0.8552   HC0 0.9287   HC3 0.9370

The generating slope is 2, and the errors have conditional mean zero at every x. That is why OLS is unbiased in this experiment despite changing variance. The mean of 4,000 fitted slopes is itself an estimate; its printed Monte Carlo standard error describes how much that average can vary. Individual fitted slopes are not all near 2. The interval coverage measures how often each constructed interval contains the generating slope.

For noise standard deviation proportional to x³, coverage is 0.8552 with the usual interval, 0.9287 with HC0, and 0.9370 with HC3. HC3 improves coverage in this run but remains below 0.95. The “usual” interval uses a common residual variance. HC0 uses squared residuals in a sandwich covariance estimate; it does not estimate each individual observation’s variance accurately from one residual. Under suitable sampling and moment conditions, the aggregate covariance estimate supports asymptotic inference. HC3 replaces each squared residual by \(r_i^2/(1-h_i)^2\), where leverage \(h_i\) is the corresponding diagonal of the projection matrix: the weight of that observation’s outcome in its own fitted value. It compensates for residual shrinkage at high-leverage points, but does not guarantee exact coverage or universal improvement.

All three interval methods use the normal cutoff 1.96. Even in the constant-variance Gaussian case, the exact classical cutoff would come from a t distribution with 98 degrees of freedom. Read the observed coverage as a finite simulation result; near 95%, its Monte Carlo standard error is about 0.0034 for 4,000 repetitions. For independent observations, HC estimates address heteroscedasticity; clustered or serially dependent observations need methods such as cluster or HAC covariance estimates. A fan-shaped residual plot suggests changing spread, but can coexist with a misspecified mean. Neither robust standard errors nor a residual plot repairs omitted-variable bias.

Training \(R^2\) under nested least-squares fits

import numpy as np

rng = np.random.default_rng(0)
n = 60
X0 = rng.normal(size=(n, 2))
y = X0 @ [2.0, -1.0] + rng.normal(size=n)
noise = rng.normal(size=(n, 57))

for extra in (0, 5, 20, 50, 57):
    X = np.hstack([X0, noise[:, :extra]])
    Xd = np.column_stack([np.ones(n), X])
    p = X.shape[1]
    assert np.linalg.matrix_rank(Xd) == p + 1
    resid = y - Xd @ np.linalg.lstsq(Xd, y, rcond=None)[0]
    r2 = 1 - resid @ resid / ((y - y.mean()) @ (y - y.mean()))
    residual_df = n - p - 1
    adjusted = f"{1 - (1 - r2) * (n - 1) / residual_df:+.4f}" if residual_df > 0 else "undefined"
    print(f"p={p:3d}  R2 {r2:.4f}   adjusted R2 {adjusted}")
# p=  2  R2 0.8239   adjusted R2 +0.8177
# p=  7  R2 0.8305   adjusted R2 +0.8077
# p= 22  R2 0.9033   adjusted R2 +0.8457
# p= 52  R2 0.9829   adjusted R2 +0.8555
# p= 59  R2 1.0000   adjusted R2 undefined

Each larger design contains all columns from the smaller designs, with the same rows, target, and intercept. For exact unpenalized least squares, adding columns cannot increase the minimized residual sum of squares: setting the new coefficients to zero reproduces the old fit. With a fixed, nonconstant target, training \(R^2\) therefore cannot fall. This guarantee does not apply to held-out \(R^2\), changed rows, or penalized fits. The extra columns here are generated independently of the target, so improved training fit need not imply better prediction.

Adjusted \(R^2=1-(1-R^2)(n-1)/(n-p-1)\) uses residual degrees of freedom when an intercept and \(p\) full-rank predictor columns are fitted. It can decrease as features are added, but can also increase after a noise feature happens to improve sample fit. It is a selection heuristic, not a guarantee of finding the generating model. At \(p=59\) and \(n=60\), the residual degrees of freedom are zero, so the code reports “undefined” without dividing by zero.

Evaluate model choices with a suitable held-out or cross-validation procedure. Report prediction error in the target’s units, such as RMSE or MAE, alongside R² and a relevant baseline. R² measures error relative to target variation in the evaluated sample; the same R² can correspond to very different absolute errors. Whether those errors are acceptable depends on the intended decision.

A coefficient that changes sign, on real data

California housing contains eight features for 20,640 census block groups. The target is a block group’s median house value, expressed in units of 100,000 dollars. The loader may download data on first use. We compare a rooms-only fit with an eight-feature fit on the same training rows, then evaluate the full model on the reserved random test split. Both coefficient comparisons use the training standardization.

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler

X, y = fetch_california_housing(return_X_y=True, as_frame=True)
names = list(X.columns)
Xtr, Xte, ytr, yte = train_test_split(
    X.to_numpy(), y.to_numpy(), test_size=0.3, random_state=0)

sc = StandardScaler().fit(Xtr)
ols = LinearRegression().fit(sc.transform(Xtr), ytr)
for name, coefficient in zip(names, ols.coef_):
    print(f"  {name:12s} {coefficient:9.4f}")
rooms = names.index("AveRooms")
beds = names.index("AveBedrms")
Ztr, Zte = sc.transform(Xtr), sc.transform(Xte)
short = LinearRegression().fit(Ztr[:, [rooms]], ytr)
print(f"rooms-only coefficient {short.coef_[0]:+.4f}   full {ols.coef_[rooms]:+.4f}")
print(f"training corr(rooms, bedrooms) {np.corrcoef(Xtr[:, rooms], Xtr[:, beds])[0, 1]:.4f}")
pred = ols.predict(Zte)
print(f"test R2 {ols.score(Zte, yte):.4f}   RMSE {np.sqrt(np.mean((yte-pred)**2)):.4f}")
print(f"training-mean baseline RMSE {np.sqrt(np.mean((yte-ytr.mean())**2)):.4f}")
#   MedInc          0.8449
#   HouseAge        0.1157
#   AveRooms       -0.2702
#   AveBedrms       0.2908
#   Population     -0.0108
#   AveOccup       -0.0281
#   Latitude       -0.8753
#   Longitude      -0.8496
# rooms-only coefficient +0.1837   full -0.2702
# training corr(rooms, bedrooms) 0.8253
# test R2 0.5926   RMSE 0.7370
# training-mean baseline RMSE 1.1547

The rooms-only coefficient is +0.1837, while the full-model AveRooms coefficient is −0.2702. They refer to different comparisons. In the full fit, increasing average rooms by one training standard deviation changes the fitted value by about −27,020 dollars while holding the other seven standardized inputs fixed. This is a block-group association in the fitted model, not the effect of adding a room to an individual house. AveBedrms and AveRooms are household averages, and the latter includes bedrooms.

AveRooms and AveBedrms are strongly correlated, but that fact alone does not explain the sign reversal or imply that the effects are impossible to estimate separately. Correlation can increase coefficient uncertainty; the conditional associations also depend on the other included predictors. A correlation below 1 is not exact collinearity, and the data do not identify merely the sum of these two coefficients. Assess uncertainty and coefficient stability rather than assigning a cause from one correlation.

On this random test split, the full model has R² 0.5926 and RMSE 0.7370 (about 73,700 dollars), versus RMSE 1.1547 for the training-mean baseline. These results do not validate a causal interpretation or establish performance in a new geographic area. A question about the effect of changing a house requires a causal design and appropriate individual-level measurements; a block-group regression alone does not answer it.

Exercises

1. An omitted variable reverses a sign. Generate data where practice reduces a score but skilled people practice more, and fit the score on practice with and without skill in the model.

Compare the fitted slopes with their population targets. The short and full models estimate different relationships.

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
n = 5000
skill = rng.normal(size=n)
practice = 2.0 * skill + 0.5 * rng.normal(size=n)      # skilled people practice more
score = 3.0 * skill - 0.5 * practice + 0.5 * rng.normal(size=n)

short = LinearRegression().fit(practice.reshape(-1, 1), score)
full = LinearRegression().fit(np.column_stack([practice, skill]), score)
print(f"score ~ practice          {short.coef_[0]:+.4f}")
print(f"score ~ practice + skill  {full.coef_[0]:+.4f}  (true -0.5)")
# score ~ practice          +0.9015
# score ~ practice + skill  -0.4873  (true -0.5)

2. Diagnose collinearity with the variance inflation factor. Build three features where two have population correlation 0.95, and compute each feature’s VIF from a regression on the others.

You should get: two features with a VIF near 8 and one near 1, and a rule for reading the number.

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(1)
n = 500
a = rng.normal(size=n)
b = 0.95 * a + np.sqrt(1 - 0.95 ** 2) * rng.normal(size=n)
c = rng.normal(size=n)
X = np.column_stack([a, b, c])

for j, name in enumerate("abc"):
    others = np.delete(X, j, axis=1)
    r2 = LinearRegression().fit(others, X[:, j]).score(others, X[:, j])
    print(f"{name}: R2 from the others {r2:.4f}   VIF {1 / (1 - r2):8.2f}")
# a: R2 from the others 0.8747   VIF     7.98
# b: R2 from the others 0.8750   VIF     8.00
# c: R2 from the others 0.0088   VIF     1.01

3. Two failure modes of extrapolation. Fit a linear model and a random forest to a logarithmic relationship observed on \([0, 10]\), then predict at \(x = 20\) and \(x = 50\).

Compare predictions at 20 and 50, then explain how each fitted function continues beyond the training inputs.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor

rng = np.random.default_rng(2)
x = rng.uniform(0, 10, 300)
y = 3 * np.log1p(x) + 0.2 * rng.normal(size=300)

lin = LinearRegression().fit(x.reshape(-1, 1), y)
rf = RandomForestRegressor(n_estimators=100, random_state=0).fit(x.reshape(-1, 1), y)

for xq in (5.0, 10.0, 20.0, 50.0):
    print(f"x={xq:5.1f}  truth {3 * np.log1p(xq):7.3f}"
          f"  linear {lin.predict([[xq]])[0]:7.3f}  forest {rf.predict([[xq]])[0]:7.3f}")
# x=  5.0  truth   5.375  linear   4.898  forest   5.374
# x= 10.0  truth   7.194  linear   7.916  forest   7.036
# x= 20.0  truth   9.134  linear  13.952  forest   7.036
# x= 50.0  truth  11.795  linear  32.059  forest   7.036