Image Author
You run a regression. The R² is 0.87. The coefficients are significant. You write it up, ship the model, and move on. Then someone asks whether you tested for homoscedasticity. You check your notes. You did not. You checked R². You checked p-values. You called it done.
That is how most people use linear regression, and it is how most people end up with models that are technically fitted but statistically unreliable. The five assumptions of linear regression are not bureaucratic checkboxes — they are the conditions under which the underlying math actually works. Violate them and your standard errors are wrong, your p-values mislead you, your confidence intervals lie, and your coefficients may not mean what you think they mean.
This article is a diagnostic checklist. For each assumption: what it means in plain language, what breaks when you violate it, how to test it in Python, and what to do when you catch a problem. Run this checklist before you interpret a single coefficient.
The five assumptions are linearity, independence of errors, homoscedasticity, normality of residuals, and no multicollinearity. Linearity and independence are foundational — if either breaks, residual-based tests for the other three can mislead you. Check them in order. Some violations also compound each other: heteroscedasticity and non-normality often appear together, and both can be symptoms of an underlying linearity problem that was never addressed.
Assumption 1: Linearity
Linear regression assumes that the relationship between each predictor and the outcome variable is linear. This is not an arbitrary restriction — it is the entire premise of the model. If the true relationship is a curve and you fit a straight line, the model will be systematically wrong in predictable ways, and no amount of tuning will fix a fundamentally misspecified structure.
Why Most People Miss It
R² is the culprit. A model can produce a decent R² on non-linear data because R² measures the proportion of variance explained, not whether the relationship is correctly specified. A curved relationship has variance, and a linear model can explain a meaningful chunk of it while still being structurally wrong. People see 0.80, feel satisfied, and never look at the residuals.
What Breaks
The coefficients are biased. Predictions are systematically too high in some ranges and too low in others. The residual plot will show a curve or arch pattern instead of random scatter — the clearest signal that the model is misspecified.
How to Detect It
Plot residuals against fitted values. If the relationship is truly linear, the residuals should form a flat, randomly scattered horizontal band around zero. Any visible curve, arch, or trend is a red flag. For multiple predictors, partial regression plots isolate each predictor’s relationship with the outcome while controlling for the others.
For a formal test, the Ramsey RESET test checks whether adding powers of the fitted values improves the model; if it does, the original specification was missing non-linear structure.
Python Code
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.diagnostic import linear_reset
np.random.seed(42)
n = 200
# Generate data with a true non-linear (quadratic) relationship
x = np.linspace(1, 10, n)
y_true = 3 * x**2 - 10 * x + 5 # True relationship is curved
y = y_true + np.random.normal(0, 15, n) # Add noise
X = sm.add_constant(x) # Add intercept column to predictor matrix
# Fit a simple linear model (deliberately misspecified -- true relationship is quadratic)
model = sm.OLS(y, X).fit()
residuals = model.resid
fitted = model.fittedvalues
# --- Plot 1: Residuals vs Fitted for the misspecified linear model ---
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(fitted, residuals, alpha=0.5, color="steelblue", edgecolors="none")
axes[0].axhline(0, color="red", linewidth=1.2, linestyle="--")
axes[0].set_title("Linear Model -- Non-linear Data\n(Notice the arch in residuals)")
axes[0].set_xlabel("Fitted Values")
axes[0].set_ylabel("Residuals")
# Fix: add a quadratic term to the model
X_quad = sm.add_constant(np.column_stack([x, x**2])) # Include x and x-squared
model_quad = sm.OLS(y, X_quad).fit()
residuals_quad = model_quad.resid
fitted_quad = model_quad.fittedvalues
axes[1].scatter(fitted_quad, residuals_quad, alpha=0.5, color="seagreen", edgecolors="none")
axes[1].axhline(0, color="red", linewidth=1.2, linestyle="--")
axes[1].set_title("Quadratic Model -- Same Data\n(Residuals are now flat and random)")
axes[1].set_xlabel("Fitted Values")
axes[1].set_ylabel("Residuals")
plt.tight_layout()
plt.savefig("linearity_residuals.png", dpi=150)
plt.show()
# --- RESET test: formally check whether non-linear terms improve the model ---
reset_result = linear_reset(model, power=2, use_f=True)
print(f"\nRESET Test on linear model:")
print(f" F-statistic: {reset_result.fvalue:.4f}")
print(f" p-value: {reset_result.pvalue:.4f}")
print(f" Verdict: {'Non-linearity detected -- model misspecified' if reset_result.pvalue < 0.05 else 'No non-linearity detected'}")
reset_result_quad = linear_reset(model_quad, power=2, use_f=True)
print(f"\nRESET Test on quadratic model:")
print(f" F-statistic: {reset_result_quad.fvalue:.4f}")
print(f" p-value: {reset_result_quad.pvalue:.4f}")
print(f" Verdict: {'Non-linearity detected' if reset_result_quad.pvalue < 0.05 else 'No non-linearity detected -- model correctly specified'}")
Assumption 2: Independence of Errors
The residuals must be independent of each other — the error in one observation should tell you nothing about the error in another. This sounds abstract until you consider how data is actually collected: students from the same classroom, patients from the same hospital, transactions from the same user, measurements taken on consecutive days. Observations that share context tend to share error structure, and that correlation between errors is exactly what this assumption prohibits.
Why Most People Miss It
Independence is taught almost exclusively in the context of time series, which leads practitioners working with cross-sectional data to assume they are safe. They are often not. Any dataset with natural clustering — surveys from the same household, employees within the same company, students within the same school — can produce correlated errors that are invisible until you specifically look for them.
As research on regression assumptions confirms, this is one of the more serious violations. When errors are correlated, standard errors are typically understated, p-values are too small, and predictors that appear significant may not be.
How to Detect It
For time-ordered data, the Durbin-Watson test checks for first-order autocorrelation in residuals. The statistic ranges from 0 to 4, with a value near 2 indicating no autocorrelation. Values closer to 0 indicate positive autocorrelation; values closer to 4 indicate negative autocorrelation. Values between 1.5 and 2.5 are generally considered acceptable.
For clustered cross-sectional data, plot residuals against the grouping variable and look for within-group patterns.
The Independence Checklist
Work through these four questions before you trust your standard errors:
- Were observations drawn independently from the population, or are there natural groupings in your sampling design?
- Is there a time ordering to your data? If yes, run the Durbin-Watson test.
- Do observations share context that could produce correlated errors — same location, same subject, same time period?
- Have you plotted residuals against any potential grouping variable to check for within-group patterns?
Python Code
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.stattools import durbin_watson
np.random.seed(42)
n = 100
# --- Example 1: Residuals with positive autocorrelation (violation) ---
# Simulate time-series data where each error is correlated with the previous one
x = np.arange(n)
errors_autocorrelated = np.zeros(n)
errors_autocorrelated[0] = np.random.normal(0, 1)
for t in range(1, n):
# Each error is 0.8 * previous error + new noise -- strong positive autocorrelation
errors_autocorrelated[t] = 0.8 * errors_autocorrelated[t - 1] + np.random.normal(0, 1)
y_autocorrelated = 2 * x + errors_autocorrelated
X = sm.add_constant(x)
model_auto = sm.OLS(y_autocorrelated, X).fit()
dw_autocorrelated = durbin_watson(model_auto.resid)
print("=== Autocorrelated Errors (Violation) ===")
print(f"Durbin-Watson statistic: {dw_autocorrelated:.4f}")
print(f"Interpretation: {'Near 2 -- no autocorrelation' if 1.5 < dw_autocorrelated < 2.5 else 'Outside 1.5-2.5 range -- autocorrelation likely present'}")
# --- Example 2: Residuals with no autocorrelation (assumption met) ---
errors_independent = np.random.normal(0, 1, n) # Each error is truly independent
y_independent = 2 * x + errors_independent
model_indep = sm.OLS(y_independent, X).fit()
dw_independent = durbin_watson(model_indep.resid)
print("\n=== Independent Errors (Assumption Met) ===")
print(f"Durbin-Watson statistic: {dw_independent:.4f}")
print(f"Interpretation: {'Near 2 -- no autocorrelation detected' if 1.5 < dw_independent < 2.5 else 'Outside 1.5-2.5 -- investigate further'}")
# --- Visual: plot residuals over time to make autocorrelation visible ---
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(model_auto.resid, color="firebrick", linewidth=0.8)
axes[0].axhline(0, color="black", linewidth=1, linestyle="--")
axes[0].set_title(f"Autocorrelated Residuals\n(DW = {dw_autocorrelated:.2f} -- violation)")
axes[0].set_xlabel("Observation order")
axes[0].set_ylabel("Residual")
axes[1].plot(model_indep.resid, color="steelblue", linewidth=0.8)
axes[1].axhline(0, color="black", linewidth=1, linestyle="--")
axes[1].set_title(f"Independent Residuals\n(DW = {dw_independent:.2f} -- assumption met)")
axes[1].set_xlabel("Observation order")
axes[1].set_ylabel("Residual")
plt.tight_layout()
plt.savefig("independence_residuals.png", dpi=150)
plt.show()
Assumption 3: Homoscedasticity
The variance of the residuals must be constant across all levels of the predictors. If you plot residuals against fitted values, they should scatter in a flat, even horizontal band. Heteroscedasticity appears as a cone or fan shape — residuals spread out as fitted values get larger, or contract, or follow some other non-constant pattern.
Why It Is Dangerous
Heteroscedasticity does not necessarily wreck your coefficient estimates — predictions can still be in the right ballpark. What it wrecks is your standard errors and, by extension, every piece of inference you draw from the model. Your p-values are wrong. Your confidence intervals are wrong. A predictor might look significant when it is not, or vice versa.
According to Statistics Solutions, a scatterplot of residuals versus predicted values is the primary diagnostic; a cone-shaped pattern confirms heteroscedasticity.
The Homoscedasticity Checklist
Before trusting any inference from your model, work through this:
- Does a plot of residuals vs. fitted values show a cone, fan, or any other non-uniform spread pattern?
- Does residual variance appear to increase or decrease with the magnitude of fitted values?
- Does the Breusch-Pagan test return a p-value below 0.05? (The null hypothesis is homoscedasticity; rejecting it means heteroscedasticity is present.)
- If any of the above are true: apply a log or square root transformation to the outcome, use weighted least squares, or use robust (HC3) standard errors.
Python Code
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan
import matplotlib.pyplot as plt
np.random.seed(42)
n = 200
x = np.linspace(1, 10, n)
# Generate heteroscedastic data: variance of errors grows with x
# This is the classic "cone" pattern -- larger x = more spread in residuals
errors_hetero = np.random.normal(0, 0.8 * x, n) # Standard deviation scales with x
y_hetero = 4 * x + 10 + errors_hetero
# Generate homoscedastic data: constant error variance regardless of x
errors_homo = np.random.normal(0, 3, n) # Constant standard deviation
y_homo = 4 * x + 10 + errors_homo
X = sm.add_constant(x)
model_hetero = sm.OLS(y_hetero, X).fit()
model_homo = sm.OLS(y_homo, X).fit()
# --- Breusch-Pagan test ---
# Returns: (LM statistic, LM p-value, F statistic, F p-value)
# We use the LM p-value (index 1) as the primary result
bp_hetero = het_breuschpagan(model_hetero.resid, model_hetero.model.exog)
bp_homo = het_breuschpagan(model_homo.resid, model_homo.model.exog)
print("=== Heteroscedastic Model ===")
print(f"Breusch-Pagan LM statistic: {bp_hetero[0]:.4f}")
print(f"p-value: {bp_hetero[1]:.4f}")
print(f"Verdict: {'Heteroscedasticity detected -- standard errors are unreliable' if bp_hetero[1] < 0.05 else 'No heteroscedasticity detected'}")
print("\n=== Homoscedastic Model ===")
print(f"Breusch-Pagan LM statistic: {bp_homo[0]:.4f}")
print(f"p-value: {bp_homo[1]:.4f}")
print(f"Verdict: {'Heteroscedasticity detected' if bp_homo[1] < 0.05 else 'Assumption satisfied -- residual variance is constant'}")
# --- Residual plots ---
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(model_hetero.fittedvalues, model_hetero.resid,
alpha=0.5, color="firebrick", edgecolors="none")
axes[0].axhline(0, color="black", linewidth=1.2, linestyle="--")
axes[0].set_title(f"Heteroscedastic Model\n(Breusch-Pagan p = {bp_hetero[1]:.4f})")
axes[0].set_xlabel("Fitted Values")
axes[0].set_ylabel("Residuals")
axes[1].scatter(model_homo.fittedvalues, model_homo.resid,
alpha=0.5, color="seagreen", edgecolors="none")
axes[1].axhline(0, color="black", linewidth=1.2, linestyle="--")
axes[1].set_title(f"Homoscedastic Model\n(Breusch-Pagan p = {bp_homo[1]:.4f})")
axes[1].set_xlabel("Fitted Values")
axes[1].set_ylabel("Residuals")
plt.tight_layout()
plt.savefig("homoscedasticity_residuals.png", dpi=150)
plt.show()
# --- Fix: robust standard errors (HC3) ---
# HC3 (MacKinnon-White) is the most commonly recommended heteroscedasticity-
# consistent covariance estimator for small to medium samples
model_robust = model_hetero.get_robustcov_results(cov_type="HC3")
print("\n=== Comparison: Regular vs Robust Standard Errors ===")
print(f"{'Predictor':<12} {'OLS SE':>10} {'HC3 Robust SE':>15}")
print("-" * 40)
for i, name in enumerate(model_hetero.model.exog_names):
print(f"{name:<12} {model_hetero.bse[i]:>10.4f} {model_robust.bse[i]:>15.4f}")
Assumption 4: Normality of Residuals
The residuals — not the predictors, not the outcome variable — must be normally distributed around zero. Testing whether your input variables are normally distributed is irrelevant to this assumption. What the model requires is that the errors in its predictions follow a normal distribution.
Why It Is Misunderstood
Two layers of confusion exist here. First, people test the wrong thing — they run normality tests on their predictors or outcome, not on the residuals from the fitted model. Second, people over-weight this assumption relative to the others. As the research confirms, violating normality is less catastrophic than violating linearity or independence. In large samples, the Central Limit Theorem provides a natural cushion. This assumption matters most in small samples where you cannot rely on that cushion.
What Breaks
In small samples: confidence intervals become unreliable, prediction intervals are off, and hypothesis tests lose their validity. In large samples, the consequences are typically mild. A mild skew in residuals with n = 500 is a non-issue; the same skew with n = 30 is a real problem.
How to Detect It
The Q-Q (quantile-quantile) plot is the primary visual diagnostic. Plot the quantiles of your residuals against those of a theoretical normal distribution — if the assumption holds, the points should follow the diagonal reference line closely. Systematic S-curves indicate skewness; heavy tails show up as points pulling away from the line at both ends.
For formal testing, the Shapiro-Wilk test is preferred for small samples (n < 50). The Kolmogorov-Smirnov test scales better with larger samples. In very large datasets, both tests will flag even trivial departures as significant, so use visual inspection alongside any formal test.
Python Code
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n = 100
x = np.linspace(0, 10, n)
# --- Case 1: Normally distributed errors (assumption met) ---
y_normal = 3 * x + 5 + np.random.normal(0, 2, n)
X = sm.add_constant(x)
model_normal = sm.OLS(y_normal, X).fit()
resid_normal = model_normal.resid
# Shapiro-Wilk test: null hypothesis is that residuals are normally distributed
# p > 0.05 = fail to reject normality (assumption plausibly met)
# p < 0.05 = reject normality (assumption violated) sw_stat, sw_p = stats.shapiro(resid_normal) print("=== Normal Residuals ===") print(f"Shapiro-Wilk statistic: {sw_stat:.4f}") print(f"p-value: {sw_p:.4f}") print(f"Verdict: {'Normality not rejected -- assumption plausibly met' if sw_p > 0.05 else 'Normality rejected -- investigate further'}")
# --- Case 2: Skewed errors (assumption violated) ---
# Exponential distribution produces right-skewed, non-normal errors
y_skewed = 3 * x + 5 + np.random.exponential(scale=3, size=n) - 3
model_skewed = sm.OLS(y_skewed, X).fit()
resid_skewed = model_skewed.resid
sw_stat_sk, sw_p_sk = stats.shapiro(resid_skewed)
print("\n=== Skewed Residuals (Violation) ===")
print(f"Shapiro-Wilk statistic: {sw_stat_sk:.4f}")
print(f"p-value: {sw_p_sk:.4f}")
print(f"Verdict: {'Normality not rejected' if sw_p_sk > 0.05 else 'Normality rejected -- residuals are non-normal'}")
# --- Q-Q plots: visual inspection ---
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# probplot draws the Q-Q plot and returns the theoretical vs sample quantiles
sm.qqplot(resid_normal, line="s", ax=axes[0], alpha=0.6,
markerfacecolor="steelblue", markeredgewidth=0)
axes[0].set_title(f"Normal Residuals -- Q-Q Plot\n(Shapiro-Wilk p = {sw_p:.4f})")
axes[0].set_xlabel("Theoretical Quantiles")
axes[0].set_ylabel("Sample Quantiles")
sm.qqplot(resid_skewed, line="s", ax=axes[1], alpha=0.6,
markerfacecolor="firebrick", markeredgewidth=0)
axes[1].set_title(f"Skewed Residuals -- Q-Q Plot\n(Shapiro-Wilk p = {sw_p_sk:.4f})")
axes[1].set_xlabel("Theoretical Quantiles")
axes[1].set_ylabel("Sample Quantiles")
plt.tight_layout()
plt.savefig("normality_qq.png", dpi=150)
plt.show()
Assumption 5: No Multicollinearity
The predictor variables must not be highly correlated with each other. Some correlation is normal and expected. High multicollinearity — where one predictor can be closely predicted from a linear combination of the others — is the problem. The model cannot reliably separate the individual effects of predictors that move together.
Why It Is Dangerous
This is the assumption whose violation is most deceptive. The model’s overall predictions can still be accurate, and R² is unaffected. The damage shows up in the coefficients: high multicollinearity inflates standard errors, drives p-values up, and depresses t-statistics. A predictor that genuinely affects the outcome can appear non-significant simply because it shares too much information with another predictor. You end up with a model that predicts well but whose individual coefficients are meaningless for interpretation or inference.
The VIF Threshold Guide
The Variance Inflation Factor measures how much the variance of a coefficient is inflated by its correlation with the other predictors. Every predictor in your model gets its own VIF.
| VIF Range | Interpretation |
|---|---|
| VIF = 1 | No correlation with other predictors |
| 1 to 5 | Moderate correlation, generally acceptable |
| 5 to 10 | Concerning — investigate further |
| Above 10 | Severe multicollinearity — urgent action required |
VIF above 10 is the widely cited threshold for severe multicollinearity, though some practitioners use 5 as a more conservative cutoff.
What to Do
Remove one of the correlated predictors. If both carry theoretical importance, combine them into a composite index or average. Apply PCA to reduce correlated predictors to uncorrelated components. Alternatively, use Ridge regression, which handles multicollinearity by design through L2 regularization.
Python Code
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
np.random.seed(42)
n = 200
# Create a dataset with two highly correlated predictors
x1 = np.random.normal(0, 1, n)
x2 = 0.95 * x1 + np.random.normal(0, 0.1, n) # x2 is almost identical to x1
x3 = np.random.normal(0, 1, n) # x3 is independent -- no issue
y = 3 * x1 + 2 * x2 + 1.5 * x3 + np.random.normal(0, 1, n)
# Assemble the predictor DataFrame (no intercept -- VIF is computed per predictor)
X = pd.DataFrame({"x1": x1, "x2": x2, "x3": x3})
def compute_vif(df):
"""
Compute the Variance Inflation Factor for each column in df.
VIF = 1 / (1 - R^2) for each predictor regressed on all others.
Higher VIF = more multicollinearity with other predictors.
Parameters:
df : DataFrame of predictors (no intercept column)
Returns:
DataFrame with predictor names and their VIF values
"""
vif_data = pd.DataFrame()
vif_data["Predictor"] = df.columns
vif_data["VIF"] = [
variance_inflation_factor(df.values, i) for i in range(df.shape[1])
]
vif_data["Status"] = vif_data["VIF"].apply(
lambda v: "Severe" if v > 10 else ("Concerning" if v > 5 else "Acceptable")
)
return vif_data.sort_values("VIF", ascending=False).reset_index(drop=True)
print("=== VIF Results (with multicollinearity) ===")
vif_results = compute_vif(X)
print(vif_results.to_string(index=False))
# --- Fix: drop one of the correlated predictors ---
X_fixed = X.drop(columns=["x2"]) # Remove x2 since it is nearly identical to x1
print("\n=== VIF Results (after dropping x2) ===")
vif_fixed = compute_vif(X_fixed)
print(vif_fixed.to_string(index=False))
# --- Show how multicollinearity inflates standard errors ---
X_with_const = sm.add_constant(X)
X_fixed_with_const = sm.add_constant(X_fixed)
model_multi = sm.OLS(y, X_with_const).fit()
model_fixed = sm.OLS(y, X_fixed_with_const).fit()
print("\n=== Standard Error Comparison: Multicollinear vs Fixed ===")
print(f"{'Predictor':<8} {'SE (with x2)':>14} {'SE (without x2)':>17}")
print("-" * 42)
for pred in ["x1", "x3"]:
se_multi = model_multi.bse[pred]
se_fixed = model_fixed.bse[pred]
print(f"{pred:<8} {se_multi:>14.4f} {se_fixed:>17.4f}")
The Full Diagnostic Checklist
Run this before interpreting any coefficient. Every row is a check. Every failure needs a fix before you trust the output.
| Assumption | Visual Check | Formal Test | Python Function | Pass Criterion | Fix If Violated |
|---|---|---|---|---|---|
| Linearity | Residuals vs. fitted — no curve | RESET test | linear_reset() | p > 0.05 | Add polynomial term, log transform, or switch model |
| Independence | Residuals vs. time — no runs | Durbin-Watson | durbin_watson() | DW between 1.5 and 2.5 | Use clustered SE, mixed models, or account for time structure |
| Homoscedasticity | Residuals vs. fitted — no cone | Breusch-Pagan | het_breuschpagan() | p > 0.05 | Log/sqrt transform, weighted LS, or HC3 robust SE |
| Normality | Q-Q plot — points on the diagonal | Shapiro-Wilk | stats.shapiro() | p > 0.05 (small n) | Log transform outcome, weighted LS, or HC3 robust SE |
| No Multicollinearity | Correlation matrix — no pairs above 0.80 | VIF | variance_inflation_factor() | All VIF below 5 (10 max) | Drop predictor, combine, PCA, or use Ridge regression |
Conclusion
Your model’s predictions are only as trustworthy as the assumptions underneath them. R² is not a validity check. Significant p-values are not a validity check. The five assumptions above are the validity check — and they are the ones that get skipped precisely because nothing in the standard model output will tell you they are broken.
Run the diagnostic checklist before you interpret a single coefficient. All five tests are in statsmodels and scipy, and the whole suite takes fewer than 20 lines of code. The difference between a model you can trust and one that looks fine until it fails is almost always in whether those 20 lines were ever run.
