Regression Diagnostics

Regression diagnostics are the systematic checks and analyses performed after fitting a regression model to identify and correct issues that compromise model accuracy, reliability, and validity. They matter when a model’s predictions deviate from expectations, when residuals exhibit patterns, when outliers skew coefficients, or when assumptions underlying the regression—linearity, independence, homoscedasticity, normality—fail to hold. Use them when a model is ready for deployment, when a new feature impacts performance, or when a stakeholder questions why the model “looks good but doesn’t work well.”

Model Fit Assessment

Evaluate how well the model explains the variation in the response variable.

Use R² and adjusted R² to quantify the proportion of variance explained, but do not rely solely on R². Compute it as:

import numpy as np

y_true = np.array([1.2, 3.4, 2.1, 5.6, 4.3])
y_pred = np.array([1.1, 3.6, 2.0, 5.8, 4.0])

ss_res = np.sum((y_true - y_pred) ** 2)
ss_tot = np.sum((y_true - y_true.mean()) ** 2)
r_squared = 1 - (ss_res / ss_tot)


n, p = len(y_true), 3  # 3 predictors
adj_r_squared = 1 - (1 - r_squared) * (n - 1) / (n - p - 1)

A high R² does not imply good fit if residuals are correlated or heteroscedastic. Always pair R² with residual plots.

Residual Analysis

Plot residuals against fitted values to diagnose common problems.

import matplotlib.pyplot as plt

residuals = y_true - y_pred
fitted_values = y_pred

plt.scatter(fitted_values, residuals)
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.title('Residual vs. Fitted Plot')
plt.show()

Look for:

• Non-constant variance (heteroscedasticity): spread of residuals increases with fitted values.
• Curvature: systematic non-linearity in the residual pattern.
• Outliers: points far from the horizontal line.

Use Cook’s distance to identify influential observations:

from statsmodels.regression.linear_model import OLS
from statsmodels.datasets.boston import load
import pandas as pd

data = load()
X, y = data.exog, data.endog
model = OLS(y, X).fit()

# Cook's distance
influence = model.get_influence()
cooks_d = influence.cooks_distance[0]

# Plot Cook's distance
plt.plot(cooks_d, 'bo')
plt.axhline(y=4 / len(cooks_d), color='red', linestyle='--', label='Threshold (4/n)')
plt.xlabel('Observation Index')
plt.ylabel("Cook's Distance")
plt.title("Influence of Observations")
plt.legend()
plt.show()

Observations with Cook’s distance > 4/n are influential. Values > 1 indicate strong influence.

Leverage and Influence

Assess how much each observation affects the regression coefficients.

Compute leverage from the hat matrix:

import numpy as np

X = np.column_stack([np.ones(len(y)), X])  # Add intercept column
h = X @ np.linalg.inv(X.T @ X) @ X.T  # Hat matrix
leverage = np.diag(h)  # Diagonal elements are leverages

High leverage points are extreme in predictor space. Combine leverage with residuals:

• High leverage + large residual → influential.
• High leverage + small residual → good fit, but pulls model toward them.

Use leverage-residual plots (also called Cook’s plots):

plt.scatter(leverage, residuals, alpha=0.6)
plt.xlabel('Leverage')
plt.ylabel('Residuals')
plt.title('Leverage vs. Residuals')
plt.axhline(y=0, color='red', linestyle='--')
plt.axvline(x=2 * X.shape[1] / len(y), color='green', linestyle='--', label='Mean leverage (2p/n)')
plt.legend()
plt.show()

Normality of Residuals

Check whether residuals follow a normal distribution.

Use Q–Q plots:

import scipy.stats as stats

stats.probplot(residuals, dist="norm", plot=plt)
plt.title("Q–Q Plot of Residuals")
plt.show()

Deviations from the diagonal line indicate:

• Heavy tails (points above line at ends): outliers or high kurtosis.
• Light tails: short tails, less extreme values.
• Skewness: S-shaped curve.

Use Shapiro-Wilk test for formal normality testing:

w, p = stats.shapiro(residuals)
print(f"Shapiro-Wilk test: W = {w:.4f}, p-value = {p:.4f}")

A p-value < 0.05 rejects normality. Do not rely on visual inspection alone.

Multicollinearity Diagnosis

Detect high correlation between predictors.

Compute Variance Inflation Factor (VIF) for each predictor:

from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = pd.DataFrame()
vif_data["Variable"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X, i) for i in range(X.shape[1])]

print(vif_data)

VIF > 5 suggests moderate multicollinearity; VIF > 10 indicates strong multicollinearity.

For a predictor with VIF = 12, the variance of its coefficient estimate is 12 times larger than it would be if the predictor were uncorrelated with others.

Use correlation matrix heatmaps:

import seaborn as sns

corr_matrix = pd.DataFrame(X).corr()
sns.heatmap(corr_matrix, annot=True, cmap="coolwarm", center=0)
plt.title("Correlation Matrix of Predictors")
plt.show()

High off-diagonal values (> 0.7 or < -0.7) indicate strong pairwise relationships.

Outlier Detection

Find observations that deviate significantly from the model.

Use studentized residuals (internally and externally Studentized):

# Internally studentized residuals
residuals_std = residuals / np.sqrt(model.mse_resid)
# Externally studentized (leave-one-out)
n = len(y)
studentized_residuals = residuals / np.sqrt(model.mse_resid * (1 - leverage))

# Identify outliers
outlier_threshold = 2.5  # or use 3 for strict criteria
outliers = np.where(np.abs(studentized_residuals) > outlier_threshold)[0]

Points with |studentized residual| > 2.5 are likely outliers. Use boxplots of residuals:

plt.boxplot(residuals)
plt.title("Boxplot of Residuals")
plt.ylabel("Residuals")
plt.show()

Outliers in residuals may be due to measurement error, data entry mistakes, or structural changes.

Model Misspecification

Check whether the assumed model form is correct.

Test for non-linearity using partial regression plots (also called added variable plots):

from statsmodels.graphics.regressionplots import plot_added_variable

fig, ax = plt.subplots()
plot_added_variable(model, exog_idx=1, ax=ax)
plt.title("Added Variable Plot for Predictor 1")
plt.show()

If the plot shows a non-linear pattern, consider adding polynomial terms or splines.

Use Ramsey RESET test to detect omitted variables or incorrect functional form:

from statsmodels.stats.diagnostic import linear_reset

reset_test = linear_reset(model, power=2, verbose=True)
print(reset_test)

A significant p-value (e.g., < 0.05) suggests the model is misspecified.

Heteroscedasticity Detection

Test for non-constant error variance.

Use Breusch-Pagan test:

from statsmodels.stats.diagnostic import het_breuschpagan

bp_test = het_breuschpagan(residuals, X)
lm_stat, lm_pval, f_stat, f_pval = bp_test

print(f"Breusch-Pagan test: LM = {lm_stat:.4f}, p-value = {lm_pval:.4f}")

A p-value < 0.05 indicates heteroscedasticity.

Use White test for a more general form:

from statsmodels.stats.diagnostic import het_white

white_test = het_white(residuals, X)
white_stat, white_pval = white_test

print(f"White test: Chi2 = {white_stat:.4f}, p-value = {white_pval:.4f}")

White test does not assume a specific form of heteroscedasticity.

Serial Correlation (Autocorrelation)

Check for correlation between residuals over time or sequence.

Use Durbin-Watson test:

from statsmodels.stats.stattools import durbin_watson

dw_stat = durbin_watson(residuals)

print(f"Durbin-Watson statistic: {dw_stat:.4f}")

# Interpretation:
if dw_stat < 1.5:
    print("Positive autocorrelation")
elif dw_stat > 2.5:
    print("Negative autocorrelation")
else:
    print("Little to no autocorrelation")

Values near 2 indicate no autocorrelation. Values < 1.5 suggest positive autocorrelation.

Use autocorrelation function (ACF) plots:

from statsmodels.tsa.stattools import acf

acf_values = acf(residuals, nlags=20)
plt.stem(acf_values, use_line_collection=True)
plt.axhline(y=0, color='red', linestyle='--')
plt.axhline(y=1.96/np.sqrt(len(residuals)), color='green', linestyle='--', label='95% CI')
plt.xlabel('Lag')
plt.ylabel('ACF')
plt.title("Autocorrelation Function of Residuals")
plt.legend()
plt.show()

Significant spikes at multiple lags indicate autocorrelation.

Diagnosing and Correcting Failure Modes

When a regression diagnostic fails, apply targeted corrections.

• Heteroscedasticity: use weighted least squares (WLS) with weights inverse to variance estimates.
• Non-normality: use robust standard errors (Huber-White) or bootstrapped confidence intervals.
• Multicollinearity: remove redundant predictors, apply ridge regression, or use principal component regression.
• Outliers: fit model with and without high-leverage points; use robust regression (e.g., MM-estimation).
• Autocorrelation: include lagged response variables or use generalized least squares (GLS).

Workflow for Production Regression Diagnostics

1. Fit model using training data.
2. Generate diagnostics: residual vs. fitted, Q–Q, leverage-residual, Cook’s distance.
3. Compute VIF, Breusch-Pagan, Durbin-Watson, Shapiro-Wilk.
4. Flag observations with Cook’s distance > 1, studentized residuals > 2.5, leverage > 2p/n.
5. Summarize diagnostics in a report or dashboard.
6. Investigate high-leverage and influential points.
7. Iterate: refine model, retrain, re-diagnose.

Common Pitfalls

• Ignoring leverage: high-leverage points affect model fit but are not always outliers.
• Confusing studentized residuals with standardized residuals: studentized residuals use leave-one-out variance estimates; standardized use overall variance.
• Overlooking serial correlation: time-series data often show autocorrelation that is ignored.
• Using R² alone: a model may have high R² but poor predictive performance due to heteroscedasticity or non-normality.
• Failing to validate diagnostics across data splits: diagnostics should be computed on training, validation, and test sets.

Final Output: Diagnostic Report

A complete regression diagnostics report includes:

• R² and adjusted R².
• Residual vs. fitted plot.
• Q–Q plot of residuals.
• Leverage-residual plot.
• Cook’s distance plot.
• Shapiro-Wilk test result.
• Breusch-Pagan test result.
• Durbin-Watson statistic.
• VIF table.
• List of influential and outlier observations.

This report is the foundation for model review, stakeholder communication, and production validation.

This page was rewritten on 17 October 2026. It replaced a templated version that was not specific to its own title. How these pages are written and checked. If anything here is wrong, tell us at [email protected] and we will correct it.