Run a regression in scikit-learn, and you get exactly two things back: a set of coefficients and a .predict() method. Ask it for a p-value on any of those coefficients, and there’s no built-in way to get one. This isn’t a missing feature scikit-learn’s maintainers overlooked; it’s a deliberate design choice they’ve stated directly: scikit-learn is built for prediction, not statistical inference, and its models optimize for minimizing prediction error rather than answering whether a given input actually matters. A long-standing GitHub feature request asking sklearn to add p-values to its linear models has sat open for years, and the standard answer in that thread is the same one this article is built around: if you want inference, use statsmodels instead.
That’s the actual gap this article covers. Statsmodels answers the questions sklearn was never built to answer: is this coefficient real or noise, does this added complexity earn its keep, can I trust these p-values in the first place, and which specific data points are quietly running the whole show. Five functions, one running example, real output from a test I ran before writing a word of this article.
We will use a simulated dataset of 200 house sales with four predictors: square footage, bedroom count, age, and distance to downtown. This dataset is built so that a couple of these checks catch something genuinely interesting rather than coming back clean every time.
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
rng = np.random.default_rng(7)
n = 200
size_sqft = rng.normal(1800, 500, n).clip(500, 4000)
bedrooms = rng.integers(1, 6, n)
age_years = rng.normal(20, 12, n).clip(0, 80)
distance_km = (0.004 * size_sqft + rng.normal(5, 3, n)).clip(0.2, None)
# noise variance grows with home size, on purpose, to give function #4
# something real to catch later in this article
noise = rng.normal(0, 1, n) * (5000 + size_sqft * 15)
price = 50000 + 120*size_sqft + 2000*bedrooms - 1800*distance_km - 300*age_years + noise
df = pd.DataFrame({"price": price, "size_sqft": size_sqft, "bedrooms": bedrooms,
"age_years": age_years, "distance_km": distance_km})
Prerequisites:
- Python 3.9+
- pip install statsmodels pandas numpy matplotlib (statsmodels 0.14.6 is the current stable release as of this writing). With that in place, here are the five functions worth knowing.
1. .summary() on a Fitted OLS Model
This is the single biggest reason people reach for statsmodels over sklearn. If you fit an ordinary least squares model with the formula API and the .summary() function, it prints a full statistical report in one call, not just coefficients.
model = smf.ols("price ~ size_sqft + bedrooms + age_years + distance_km", data=df).fit()
print(model.summary())
Running this on the housing data returns an R-squared of 0.720, meaning the four predictors together explain about 72% of the variation in price, along with an F-statistic of 125.6 and a vanishingly small p-value confirming the model as a whole is meaningfully better than no model at all. But the coefficient table is where the real value shows up. size_sqft, age_years, and distance_km all come back with p-values under 0.001 — genuinely significant predictors. bedrooms, on the other hand, comes back with a p-value of 0.808, nowhere close to significant, meaning once square footage is already accounted for, bedroom count isn’t adding real predictive information of its own. That’s a finding sklearn’s LinearRegression would never surface, since it would happily hand back a nonzero coefficient for bedrooms and let you assume it mattered.
What this tells you: not just what the model predicts, but which inputs are actually doing work and which ones are along for the ride.
2. anova_lm for Formal Nested Model Comparison
A higher R-squared on a more complex model doesn’t automatically mean the added complexity was worth it. R-squared mechanically increases (or at worst stays flat) every time you add a variable, whether or not that variable means anything. anova_lm answers the actual question: does the added complexity produce a statistically meaningful improvement, or just a cosmetic one?
from statsmodels.stats.anova import anova_lm
reduced_model = smf.ols("price ~ size_sqft", data=df).fit()
full_model = smf.ols("price ~ size_sqft + bedrooms + age_years + distance_km", data=df).fit()
comparison = anova_lm(reduced_model, full_model)
print(comparison)
This runs an F-test directly comparing the simple, size-only model against the full four-predictor model. The result: F = 10.87, p < 0.001. That is real and statistically grounded. The extra predictors as a group earn their place, even though Section 1 already showed that one of them — bedrooms individually — isn’t pulling meaningful weight on its own. Both things are true at once, and anova_lm is what lets you test the group-level claim formally instead of eyeballing whether R-squared moved.
What this tells you: whether a more complex model is actually justified, tested directly rather than assumed from a marginally higher R-squared.
3. variance_inflation_factor for Multicollinearity
If two or more predictors are highly correlated with each other, a regression can’t cleanly separate their individual effects, and the coefficients (and their p-values) become unstable — sometimes flipping sign entirely if you add or remove one observation. The variance inflation factor measures exactly this: how much a coefficient’s variance is inflated because its predictor is correlated with the others in the model.
from statsmodels.stats.outliers_influence import variance_inflation_factor X = df[["size_sqft", "bedrooms", "age_years", "distance_km"]] X_const = sm.add_constant(X) vif_data = pd.DataFrame() vif_data["feature"] = X_const.columns vif_data["VIF"] = [variance_inflation_factor(X_const.values, i) for i in range(X_const.shape[1])] print(vif_data)
The common rule of thumb treats a VIF above 5 as worth investigating and above 10 as a real problem. Running this on the housing data, every predictor came back under 1.4. That’s a clean result, and it’s worth including here specifically because a red-flag check clearing a model is just as useful an example as one catching a problem. In this case, distance_km was deliberately built to correlate somewhat with size_sqft, and the VIF check confirms that correlation isn’t strong enough to actually threaten the model’s coefficient estimates.
What this tells you: whether your predictors are independent enough that each coefficient means what it claims to mean, rather than being tangled up with another variable in the model.
4. het_breuschpagan for Non-Constant Error Variance
Every p-value in the summary table rests on an assumption that the model’s errors have roughly constant variance across the full range of predictions — a property called homoscedasticity. Violate that assumption and the standard errors, and therefore the p-values, are no longer trustworthy, even if the coefficients themselves are still roughly right. The Breusch-Pagan test checks this assumption directly instead of leaving you to eyeball a residual plot and guess.
from statsmodels.stats.diagnostic import het_breuschpagan
bp_test = het_breuschpagan(model.resid, model.model.exog)
labels = ["LM Statistic", "LM p-value", "F-Statistic", "F p-value"]
for label, value in zip(labels, bp_test):
print(f"{label}: {value:.5f}")
The housing dataset was deliberately built with noise that grows alongside home size (bigger, more expensive homes have more volatile pricing, which is realistic), and the test caught it cleanly: an LM p-value of 0.00001. That’s a clear violation of the constant-variance assumption.
5. OLSInfluence and Cook’s Distance
Every check so far has operated at the level of predictors — which variables matter, whether the extra ones earn their place, whether they’re tangled together, and whether the errors behave. This last one operates at the level of individual data points: which specific rows in your dataset are disproportionately steering the result.
from statsmodels.stats.outliers_influence import OLSInfluence
influence = OLSInfluence(model)
cooks_d = influence.cooks_distance[0]
threshold = 4 / len(df)
n_influential = (cooks_d > threshold).sum()
print(f"Threshold (4/n): {threshold:.4f}")
print(f"Influential points flagged: {n_influential} out of {len(df)}")
Cook’s distance measures how much the model’s fitted values would shift if a given point were removed entirely. Using the common 4/n threshold, this flagged 13 of the 200 points in the housing dataset as meaningfully influential, and one point in particular stood out with a Cook’s distance more than five times higher than the next highest value, visible as the darkest, largest point. That’s a concrete lead: go look at that exact row before trusting the model as final, since a single unusual sale (an oddly cheap mansion, a data entry error, a genuine outlier) can quietly reshape coefficients that get reported as if they applied evenly across the whole dataset.
What this tells you: not just whether the model is trustworthy in aggregate, but whether that trust is being carried disproportionately by a handful of specific rows.
Putting the Five Together
Used in sequence, these aren’t five disconnected diagnostics; they’re one coherent workflow. Fit the model and read the summary to see which predictors are doing real work. Run anova_lm to confirm the added complexity was worth including in the first place. Check VIF to make sure those predictors aren’t secretly measuring the same thing. Run Breusch-Pagan to confirm the p-values you just relied on are actually valid. Then check Cook’s distance to make sure no single row is quietly writing the whole story. Skip any one of these, and you’re trusting a number — a p-value, an R-squared, a coefficient — without checking the assumption that the number depends on.
Let’s test all five statsmodels functions using one coherent synthetic dataset (simulated house prices), so results stay comparable across every section.
"""
Tests all 5 statsmodels functions for the article, using one coherent
synthetic dataset (simulated house prices) so results stay comparable
across every section.
"""
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.outliers_influence import variance_inflation_factor, OLSInfluence
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.anova import anova_lm
rng = np.random.default_rng(7)
n = 200
size_sqft = rng.normal(1800, 500, n).clip(500, 4000)
bedrooms = rng.integers(1, 6, n)
age_years = rng.normal(20, 12, n).clip(0, 80)
# distance_to_downtown correlated with size (older, bigger suburban homes are farther out)
distance_km = 0.004 * size_sqft + rng.normal(5, 3, n)
distance_km = distance_km.clip(0.2, None)
# true price model: size and distance matter, bedrooms barely matters once size is controlled,
# and noise variance grows with size (heteroskedasticity, on purpose, to demonstrate flag 4 later)
noise = rng.normal(0, 1, n) * (5000 + size_sqft * 15)
price = (
50000
+ 120 * size_sqft
+ 2000 * bedrooms
- 1800 * distance_km
- 300 * age_years
+ noise
)
df = pd.DataFrame({
"price": price,
"size_sqft": size_sqft,
"bedrooms": bedrooms,
"age_years": age_years,
"distance_km": distance_km,
})
print("=" * 70)
print("1. OLS SUMMARY")
print("=" * 70)
model = smf.ols("price ~ size_sqft + bedrooms + age_years + distance_km", data=df).fit()
print(model.summary())
print("\n" + "=" * 70)
print("2. ANOVA_LM (nested model comparison)")
print("=" * 70)
reduced_model = smf.ols("price ~ size_sqft", data=df).fit()
full_model = smf.ols("price ~ size_sqft + bedrooms + age_years + distance_km", data=df).fit()
comparison = anova_lm(reduced_model, full_model)
print(comparison)
print("\n" + "=" * 70)
print("3. VARIANCE INFLATION FACTOR")
print("=" * 70)
X = df[["size_sqft", "bedrooms", "age_years", "distance_km"]]
X_const = sm.add_constant(X)
vif_data = pd.DataFrame()
vif_data["feature"] = X_const.columns
vif_data["VIF"] = [variance_inflation_factor(X_const.values, i) for i in range(X_const.shape[1])]
print(vif_data)
print("\n" + "=" * 70)
print("4. BREUSCH-PAGAN TEST")
print("=" * 70)
bp_test = het_breuschpagan(model.resid, model.model.exog)
labels = ["LM Statistic", "LM p-value", "F-Statistic", "F p-value"]
for label, value in zip(labels, bp_test):
print(f"{label}: {value:.5f}")
print("\n" + "=" * 70)
print("5. COOK'S DISTANCE / INFLUENCE")
print("=" * 70)
influence = OLSInfluence(model)
cooks_d = influence.cooks_distance[0]
threshold = 4 / n
n_influential = (cooks_d > threshold).sum()
print(f"Threshold (4/n): {threshold:.4f}")
print(f"Number of influential points flagged: {n_influential} out of {n}")
print(f"Top 5 Cook's distance values:\n{pd.Series(cooks_d).sort_values(ascending=False).head()}")
Output
Residuals vs Fitted | Image by Author
Wrapping Up
The honest way to describe the difference between these two libraries isn’t that one is better than the other; it’s that they’re answering different questions. Sklearn’s .fit() and .predict() are the right tools when the goal is a working prediction and nothing else. The moment the real question becomes whether a result is trustworthy — and why it came out the way it did — these five functions are what actually get you there, and none of them requires leaving Python or reaching for a separate statistical package to do it.
