Linear Regression Assumptions: What Breaks and How to Check

“What are the assumptions of linear regression?” is asked in almost every data science interview, from new-grad screens to senior loops. Reciting “linearity, independence, normality, homoscedasticity” earns little. What interviewers listen for is the next layer: which violations make the coefficients wrong and which only make the p-values wrong, how you would detect each one, and what you would do about it. That is the knowledge that stops you from presenting a regression to stakeholders that says each advertising dollar returns four, when the real answer is two.

Before you start

You should know what a regression line, a coefficient and a residual (observed minus fitted value) are, and have met standard errors and confidence intervals. The code uses Python 3.14’s statistics.linear_regression, correlation and fmean, plus a short hand-written two-predictor fit; in practice you would use statsmodels or scikit-learn. Outputs in comments come from real runs with the seeds shown, with the blocks run in order.

The short answer

Ordinary least squares assumes: the mean of y is linear in the features; errors are independent; errors have constant variance; errors are roughly normal (only important for small-sample inference); there is no perfect multicollinearity; and, most importantly, errors are uncorrelated with the features (exogeneity, which omitted confounders break). Nonlinearity and omitted variables bias the coefficients themselves. Heteroscedasticity and dependence leave coefficients unbiased but make standard errors and p-values wrong. I check with residual plots, fix variance problems with robust or clustered standard errors, and treat causal claims from regression with care.

How it works

OLS picks the line that minimises the sum of squared residuals. For one feature, the slope is cov(x, y) / var(x) and the line passes through the means. Each assumption protects a different part of the output:

Interviewers love the distinction in the second column, because it decides whether a problem changes your conclusion or only your confidence in it.

Step-by-step walkthrough

Step 1: Fit, then look at residuals instead of the summary

Anscombe’s quartet is four small datasets built in 1973 to make this point.

import math, random
from statistics import linear_regression, correlation, fmean

x = [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5]
quartet = {
    "I": (x, [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68]),
    "II": (x, [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74]),
    "III": (x, [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73]),
    "IV": ([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8], [6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89]),
}
for name, (xs, ys) in quartet.items():
    fit = linear_regression(xs, ys)
    print(name, round(fit.slope, 3), round(fit.intercept, 2), round(correlation(xs, ys), 2))
# I 0.5 3.0 0.82
# II 0.5 3.0 0.82
# III 0.5 3.0 0.82
# IV 0.5 3.0 0.82

Identical slopes, intercepts and correlations. Plotted, set I is a reasonable linear cloud, II is a smooth curve, III is a perfect line with one outlier dragging it, and IV is a vertical stack of points plus one point that alone determines the slope. A regression summary cannot tell these apart; residuals can.

Step 2: Check linearity with the residual pattern

Sort set II by x and print its residuals:

xs, ys = quartet["II"]
fit = linear_regression(xs, ys)
print([round(y - (fit.intercept + fit.slope * xi), 2) for xi, y in sorted(zip(xs, ys))])
# [-1.9, -0.76, 0.13, 0.76, 1.14, 1.27, 1.14, 0.76, 0.13, -0.76, -1.9]

If the model were right, residuals would scatter randomly around zero. These rise and fall in a perfect arch: negative at both ends, positive in the middle. The line is systematically wrong at every x, so the slope describes nothing real. Adding an x² term, transforming a variable or using a flexible model fixes the structure. Note that linearity is about the parameters: y = a + b·x + c·x² is still linear regression.

Step 3: Check constant variance, and use robust standard errors

In many business datasets, spread grows with size: big customers vary more than small ones. Here the noise standard deviation grows with x. We simulate 4,000 studies and check how often each 95% interval for the slope contains the true slope of 2.

def slope_ses(xs, ys):
    fit = linear_regression(xs, ys)
    n, xbar = len(xs), fmean(xs)
    sxx = sum((v - xbar) ** 2 for v in xs)
    resid = [y - fit.intercept - fit.slope * v for v, y in zip(xs, ys)]
    classic = math.sqrt(sum(e * e for e in resid) / (n - 2) / sxx)
    robust = math.sqrt(sum((v - xbar) ** 2 * e * e for v, e in zip(xs, resid)) / sxx ** 2 * n / (n - 2))
    return fit.slope, classic, robust

rng = random.Random(12)
def sample(n=200):
    xs = [rng.uniform(0, 10) for _ in range(n)]
    ys = [5 + 2.0 * v + rng.gauss(0, 0.3 + v ** 1.5) for v in xs]   # noise grows with x
    return xs, ys

slopes, hits_classic, hits_robust = [], 0, 0
for _ in range(4000):
    b, se_c, se_r = slope_ses(*sample())
    slopes.append(b)
    hits_classic += abs(b - 2.0) <= 1.96 * se_c
    hits_robust += abs(b - 2.0) <= 1.96 * se_r
print(round(fmean(slopes), 3), hits_classic / 4000, hits_robust / 4000)   # 2.004 0.8965 0.946

The slope is unbiased (average 2.004), but the classic standard error is too small: its “95%” interval covers only 90% of the time. The heteroscedasticity-robust (HC1, or White) standard error uses each observation’s own squared residual, weighted by how far its x sits from the mean, and restores 95%. Robust standard errors are cheap insurance, which is why many analysts use them by default.

Step 4: Check independence, especially with repeated measures

Event data usually has many rows per user. If a feature varies by user (plan, country, signup cohort) and each user has their own baseline, rows from the same user are correlated, and treating 1,000 rows as 1,000 independent observations overstates the evidence. Here the feature has no effect:

def slope_and_se(xs, ys):
    fit = linear_regression(xs, ys)
    n, xbar = len(xs), fmean(xs)
    sxx = sum((v - xbar) ** 2 for v in xs)
    rss = sum((y - fit.intercept - fit.slope * v) ** 2 for v, y in zip(xs, ys))
    return fit.slope, math.sqrt(rss / (n - 2) / sxx)

crng = random.Random(17)
def one_study(users=50, per_user=20):
    xs, ys, ux, uy = [], [], [], []
    for _ in range(users):
        feature = crng.uniform(0, 1)            # user-level feature, no true effect
        user_effect = crng.gauss(0, 1)          # shared by all of this user's rows
        rows = [user_effect + crng.gauss(0, 1) for _ in range(per_user)]
        xs += [feature] * per_user
        ys += rows
        ux.append(feature)
        uy.append(fmean(rows))
    b, se = slope_and_se(xs, ys)
    bu, seu = slope_and_se(ux, uy)
    return abs(b / se) > 1.96, abs(bu / seu) > 2.011   # t critical values for 998 and 48 df

results = [one_study() for _ in range(10_000)]
print(fmean(r for r, _ in results), fmean(u for _, u in results))   # 0.5535 0.052

The row-level regression finds a “significant” effect in 55% of studies where none exists. Aggregating to one row per user (the unit that actually varies) brings it back to 5%. Cluster-robust standard errors or mixed models do the same job without discarding row-level detail. Time series have the same problem through autocorrelation.

Step 5: Check multicollinearity with VIF

When two features carry nearly the same information, OLS cannot decide how to split the credit between them. Predictions stay fine; individual coefficients become unstable.

def two_predictor_ols(x1, x2, y):
    m1, m2, my = fmean(x1), fmean(x2), fmean(y)
    a = [v - m1 for v in x1]
    b = [v - m2 for v in x2]
    c = [v - my for v in y]
    s11, s22 = sum(u * u for u in a), sum(v * v for v in b)
    s12 = sum(u * v for u, v in zip(a, b))
    s1y, s2y = sum(u * w for u, w in zip(a, c)), sum(v * w for v, w in zip(b, c))
    det = s11 * s22 - s12 ** 2
    return (s22 * s1y - s12 * s2y) / det, (s11 * s2y - s12 * s1y) / det

mrng = random.Random(30)
for _ in range(4):
    visits = [mrng.gauss(10, 2) for _ in range(100)]
    sessions = [v + mrng.gauss(0, 0.2) for v in visits]      # nearly a copy of visits
    revenue = [1.0 * v + 1.0 * s + mrng.gauss(0, 2) for v, s in zip(visits, sessions)]
    b1, b2 = two_predictor_ols(visits, sessions, revenue)
    r = correlation(visits, sessions)
    print(round(b1, 2), round(b2, 2), round(b1 + b2, 2), round(1 / (1 - r * r)))
# -0.84 2.82 1.97 150
# 2.65 -0.77 1.88 70
# 3.34 -1.58 1.77 83
# 0.52 1.43 1.95 56

Both true coefficients are 1. Across four samples the estimates swing from -1.58 to 3.34 and even flip sign, while their sum stays near 2. The variance inflation factor, 1 / (1 - R²) where R² comes from regressing one feature on the others (here simply r²), measures how much collinearity inflates a coefficient’s variance; values above 5 to 10 are a warning, and these are 56 to 150. Combine the features, drop one, or use ridge regression if you need stable coefficients.

Worked scenario

A marketing analyst regresses weekly sales on ad spend and reports that every extra unit of spend returns 4.37 units of sales. But spend is planned around demand: the team spends more in high season, and high season sells more on its own.

rng = random.Random(3)
season = [rng.gauss(0, 1) for _ in range(5000)]           # seasonal demand
ads = [0.8 * s + rng.gauss(0, 0.6) for s in season]        # spend follows demand
sales = [2.0 * a + 3.0 * s + rng.gauss(0, 1) for a, s in zip(ads, season)]
print(round(linear_regression(ads, sales).slope, 3))                   # 4.369
print([round(v, 3) for v in two_predictor_ols(ads, season, sales)])    # [2.034, 2.964]
r = correlation(ads, season)
print(round(r, 3), round(1 / (1 - r * r), 2))                          # 0.796 2.73

Seasonality is an omitted confounder: it drives both spend and sales, so the simple regression credits ads with the season’s effect. The omitted variable bias formula predicts it exactly: the bias is the omitted coefficient (3) times the slope of season on ads (0.8 / 1.0 = 0.8), so 2 + 2.4 = 4.4. Adding seasonality recovers about 2.03. The VIF of 2.73 shows the two features are correlated but not dangerously so. The fixed report includes seasonality and says plainly that other unmeasured drivers (promotions, competitor activity) may still bias the estimate; a geo-holdout experiment would measure the ad effect directly.

Common mistake

• “y must be normally distributed.” The assumption is about the errors, conditional on the features, and it matters only for small-sample inference. A skewed y can have perfectly normal residuals, or the reverse.
• “Heteroscedasticity biases the coefficients.” It does not; it breaks the standard errors. Nonlinearity and omitted variables bias coefficients.
• Treating a high R² as proof the model is right. Anscombe’s sets II and IV have the same R² as set I.
• Reading a coefficient as a causal effect without arguing that confounders are controlled.
• Dropping a “non-significant” feature that is collinear with another and concluding it does not matter; the pair may matter a lot jointly.

Verify the behavior

The Anscombe check is a one-liner worth keeping in any regression notebook: assert the summaries agree, then look at residuals. For the robust standard error, rerun the Step 3 simulation with constant noise (replace 0.3 + v ** 1.5 with 3.0) and confirm both intervals cover close to 95%: robust errors cost little when they are not needed.

for name, (xs, ys) in quartet.items():
    fit = linear_regression(xs, ys)
    assert abs(fit.slope - 0.5) < 0.01 and abs(fit.intercept - 3.0) < 0.01, name
    worst = max(abs(y - fit.intercept - fit.slope * v) for v, y in zip(xs, ys))
    print(name, round(worst, 2))
# I 1.92
# II 1.9
# III 3.24
# IV 1.84

The largest residual singles out set III’s outlier; set IV’s influential point has a small residual precisely because it pulls the line through itself, which is why leverage (Cook’s distance) is a separate check.

Follow-up questions

• How do you interpret a coefficient after a log transform of y? A one-unit increase in x multiplies y by exp(b), roughly a 100 × b percent change for small b.
• What is the difference between R² and adjusted R²? Adjusted R² penalises extra features, so it can fall when a useless feature is added; neither measures causal validity.
• What are leverage and influence? Leverage is how unusual a point’s x values are; influence (Cook’s distance) combines leverage with residual size to measure how much the fit moves if the point is dropped.
• When would you use ridge or lasso instead? When collinearity or many features make OLS coefficients unstable; they add a little bias to cut variance, and lasso can zero out features.

Interview exercise

A house price model includes square footage, bedrooms and bathrooms. The bathroom coefficient is negative, even though bathrooms and price are positively correlated. A PM asks whether adding a bathroom lowers a home’s value. What do you tell them?

Answer and reasoning

A regression coefficient is a conditional effect: the change in predicted price for one more bathroom holding square footage and bedrooms fixed. At fixed size, another bathroom means less space elsewhere, so a negative conditional coefficient can be sensible even though bathrooms and price rise together overall. Collinearity makes it worse: bathrooms, bedrooms and footage are strongly correlated, so I would check the VIFs and expect the bathroom coefficient’s interval to be wide and possibly include zero. Finally, the model describes associations in existing homes, not the effect of renovating; houses with extra bathrooms may differ in age or location. My answer: the data does not say a bathroom lowers value, and to estimate a renovation effect I would compare similar homes before and after renovations or use a model built for that question.

Continue learning

• Practise the theory in the Statistics & Data Science chapter and test yourself with the data science MCQs.
• See how omitted confounders mislead pooled comparisons in Simpson’s paradox and correlation vs causation.
• Connect regularisation to unstable coefficients in overfitting and L1 vs L2 regularization.
• Read the NIST/SEMATECH handbook on linear least squares regression and how to tell if a model fits.