Linear Regression for Data Scientists: Fit, Interpret, and Diagnose

Linear regression looks simple, and that is why it fails quietly. A model with a high R² can still violate every assumption that makes its coefficients trustworthy. This guide covers the full workflow: fit the model, read the coefficients correctly, and run the diagnostics that separate a defensible model from a misleading one.

Quick Takeaways

• Fit with statsmodels when you need inference (p-values, confidence intervals). Use scikit-learn when you need prediction pipelines.
• Interpret each coefficient as the expected change in y for a one-unit change in xⱼ, holding all other predictors constant.
• Diagnose with residual plots, VIF, Q-Q plots, and heteroscedasticity tests before trusting any p-value.
• Fix violations with transformations, robust standard errors, or regularization (Ridge, Lasso).

What Linear Regression Models

Ordinary Least Squares (OLS) fits a line (or hyperplane) that minimizes the sum of squared residuals:

ŷ = β₀ + β₁x₁ + β₂x₂ + … + βₚxₚ

The estimator solves β̂ = (XᵀX)⁻¹Xᵀy. Each βⱼ is a partial effect: the change in the expected response per unit change in xⱼ, with other features held fixed.

The Five Assumptions (LINE + No Multicollinearity)

Normality of errors matters for inference in small samples. With large n, the central limit theorem protects your coefficient estimates. Linearity and independence matter at every sample size.

Step 1: Data Exploration and Preparation

Start with the target distribution, missingness, and pairwise relationships. This example uses the California housing dataset.

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split

# Load data as a DataFrame
data = fetch_california_housing(as_frame=True)
df = data.frame

# Quick audit: dtypes, nulls, summary statistics
print(df.info())
print(df.describe().T)

# Correlation of every feature with the target
print(df.corr()["MedHouseVal"].sort_values(ascending=False))

Check for three issues before modeling:

1. Skewed target. A right-skewed y often produces fanning residuals. Consider log(y).
2. Outliers and leverage points. A few extreme rows can dominate the fit.
3. Capped values. This dataset caps house values at 5.0, which creates a ceiling effect.

# Remove the censored ceiling rows and log-transform skewed predictors
df = df[df["MedHouseVal"] < 5.0].copy()
df["log_AveRooms"] = np.log(df["AveRooms"])
df["log_Population"] = np.log(df["Population"])

# Split BEFORE any fitting to prevent leakage
X = df[["MedInc", "HouseAge", "log_AveRooms", "AveBedrms", "log_Population", "AveOccup"]]
y = df["MedHouseVal"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

Step 2: Fit the Model

Option A: statsmodels for Inference

# add_constant appends the intercept column; OLS does not add one by default
X_train_sm = sm.add_constant(X_train)
model = sm.OLS(y_train, X_train_sm).fit()

print(model.summary())

Option B: scikit-learn for Prediction

from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score

pipe = make_pipeline(StandardScaler(), LinearRegression())

# 5-fold cross-validated RMSE (sklearn returns negative scores)
cv_rmse = -cross_val_score(
    pipe, X_train, y_train, cv=5, scoring="neg_root_mean_squared_error"
)
print(f"CV RMSE: {cv_rmse.mean():.3f} ± {cv_rmse.std():.3f}")

Step 3: Interpret the Output

Read the summary() table in this order.

Coefficients (coef)

A coefficient of 0.42 on MedInc means each one-unit rise in median income (in $10k) is associated with a 0.42 increase in median house value (in $100k), holding other predictors constant. Association is not causation unless your design supports it.

Log Transformations Change the Reading

Fit Statistics

Standard Errors, t-values, and Confidence Intervals

The t-statistic is β̂ / SE(β̂). A p-value < 0.05 means the data are inconsistent with a true slope of zero. It does not measure effect size. Always report the 95% confidence interval (model.conf_int()) alongside the point estimate.

# Tidy coefficient table with confidence intervals
results = pd.concat(
    [model.params, model.conf_int(), model.pvalues], axis=1
)
results.columns = ["coef", "ci_low", "ci_high", "p_value"]
print(results.round(4))

Step 4: Diagnose the Model

Diagnostics decide whether Step 3 means anything.

4.1 Linearity and Equal Variance: Residuals vs. Fitted

import matplotlib.pyplot as plt

fitted = model.fittedvalues
resid = model.resid

fig, ax = plt.subplots(figsize=(7, 4))
ax.scatter(fitted, resid, alpha=0.3, s=8)
ax.axhline(0, color="red", linestyle="--")
ax.set_xlabel("Fitted values")
ax.set_ylabel("Residuals")
ax.set_title("Residuals vs. Fitted")
plt.show()

4.2 Normality: Q-Q Plot and Jarque-Bera

from statsmodels.stats.stattools import jarque_bera

sm.qqplot(resid, line="45", fit=True)
plt.show()

jb_stat, jb_p, skew, kurt = jarque_bera(resid)
print(f"Jarque-Bera p-value: {jb_p:.4f}, skew={skew:.2f}, kurtosis={kurt:.2f}")

Points that peel away from the 45-degree line at the tails signal heavy tails. With thousands of rows, Jarque-Bera rejects almost always. Judge the Q-Q plot’s severity, not just the p-value.

4.3 Heteroscedasticity: Breusch-Pagan Test

from statsmodels.stats.diagnostic import het_breuschpagan

bp_stat, bp_p, _, _ = het_breuschpagan(resid, model.model.exog)
print(f"Breusch-Pagan p-value: {bp_p:.4g}")

# If p < 0.05, refit with heteroscedasticity-consistent standard errors
robust = sm.OLS(y_train, X_train_sm).fit(cov_type="HC3")
print(robust.summary().tables[1])

Robust errors leave the coefficients unchanged. They correct the standard errors, which fixes the p-values and confidence intervals.

4.4 Multicollinearity: Variance Inflation Factor

from statsmodels.stats.outliers_influence import variance_inflation_factor

vif = pd.DataFrame({
    "feature": X_train_sm.columns,
    "VIF": [
        variance_inflation_factor(X_train_sm.values, i)
        for i in range(X_train_sm.shape[1])
    ],
})
print(vif[vif["feature"] != "const"].round(2))

High VIF inflates standard errors and makes coefficients flip signs between samples. It does not hurt predictive accuracy much, so it matters most when you interpret coefficients.

4.5 Independence: Durbin-Watson

from statsmodels.stats.stattools import durbin_watson
print(f"Durbin-Watson: {durbin_watson(resid):.2f}")

Values near 2.0 indicate no first-order autocorrelation. Values below 1.5 or above 2.5 warrant action, especially for time-ordered data. Use statsmodels.tsa models or cov_type="HAC" for time series.

4.6 Influential Observations: Cook’s Distance and Leverage

influence = model.get_influence()
cooks_d = influence.cooks_distance[0]
leverage = influence.hat_matrix_diag

# Common rule of thumb: investigate points with D > 4/n
threshold = 4 / len(X_train)
flagged = np.where(cooks_d > threshold)[0]
print(f"{len(flagged)} points exceed Cook's D threshold of {threshold:.5f}")

Do not delete flagged points automatically. Inspect them. Data entry errors deserve removal. Legitimate extreme values deserve a robust method such as sm.RLM or Huber regression.

Step 5: Evaluate on Held-Out Data

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

X_test_sm = sm.add_constant(X_test)
pred = model.predict(X_test_sm)

rmse = np.sqrt(mean_squared_error(y_test, pred))
mae = mean_absolute_error(y_test, pred)
r2 = r2_score(y_test, pred)
print(f"Test RMSE={rmse:.3f}  MAE={mae:.3f}  R²={r2:.3f}")

Compare train and test metrics. A large gap points to overfitting, which is common when you add many interaction or polynomial terms.

Fixing Violations: Decision Guide

Regularization in Practice

from sklearn.linear_model import RidgeCV, LassoCV

# Alphas searched on a log grid; scaling is mandatory for penalized models
ridge = make_pipeline(
    StandardScaler(), RidgeCV(alphas=np.logspace(-3, 3, 50))
).fit(X_train, y_train)

lasso = make_pipeline(
    StandardScaler(), LassoCV(cv=5, random_state=42)
).fit(X_train, y_train)

print("Ridge alpha:", ridge[-1].alpha_)
print("Lasso nonzero coefs:", (lasso[-1].coef_ != 0).sum())

Real-World Use Cases

Use Case 1: Predicting Customer Lifetime Value (CLV)

A retailer regresses 12-month spend on recency, order count, average basket size, and acquisition channel. Spend is right-skewed, so the analyst models log(spend) and reads coefficients as percentage effects. A Breusch-Pagan test flags fanning residuals, so the team reports HC3 confidence intervals to marketing. The model ranks channels by incremental value with defensible uncertainty ranges.

Use Case 2: Pricing Elasticity

An e-commerce team fits a log-log model of units sold on price, controlling for seasonality and promotions. The price coefficient is the price elasticity. A value of -1.8 means a 1% price increase predicts a 1.8% drop in demand. VIF checks catch the overlap between promotion flags and price discounts before the coefficients reach a pricing committee.

Use Case 3: Sensor Calibration in Manufacturing

An engineer regresses a reference measurement on raw readings from a low-cost sensor. Durbin-Watson reveals strong autocorrelation because readings are time-ordered, so the engineer applies HAC standard errors. Cook’s distance isolates three readings from a known equipment fault, which are excluded with documentation.

Common Mistakes

1. Forgetting the intercept in statsmodels. sm.OLS fits through the origin unless you call add_constant.
2. Chasing R². It never decreases as you add features. Use adjusted R² and held-out error.
3. Interpreting standardized and raw coefficients interchangeably. Standardized coefficients compare feature importance. Raw coefficients carry units.
4. Extrapolating. The fitted line is valid only within the observed range of X.
5. Ignoring leakage. Fit scalers and imputers on training data only. Pipelines enforce this.
6. Treating p-values as importance. With large n, trivial effects become significant.

FAQ

What is a good R² value for linear regression?

No universal threshold exists. In physical sciences, R² > 0.9 is common. In social science and marketing, 0.2–0.5 can be useful. Judge the model by held-out RMSE, residual diagnostics, and whether predictions improve decisions.

How do you check linear regression assumptions in Python?

Plot residuals against fitted values for linearity and equal variance. Use a Q-Q plot for normality, the Durbin-Watson statistic for independence, and VIF for multicollinearity. The statsmodels.stats.diagnostic module provides formal tests such as Breusch-Pagan.

What is the difference between R² and adjusted R²?

R² measures the share of variance explained and rises with every added predictor. Adjusted R² penalizes extra predictors, so it increases only when a new feature improves the model more than chance would.

When should you use Ridge or Lasso instead of OLS?

Use Ridge when predictors are highly correlated (high VIF) and you want stable coefficients. Use Lasso when you suspect many irrelevant features and want automatic selection. Use plain OLS when you have few, well-behaved predictors and need classical inference.