TopolearnInference is only as good as the model's assumptions, and diagnostics tell you when they fail.
Every p-value and interval from OLS rests on assumptions: a correct functional form, independent errors, constant variance, roughly normal residuals, and no single point dominating. statsmodels has a test or plot for each.
Functional form and variance: plot res.resid against res.fittedvalues. A curve suggests a missing transformation or term; a funnel shape suggests heteroscedasticity. Test it with het_breuschpagan(res.resid, res.model.exog) from statsmodels.stats.diagnostic. linear_reset(res) tests for omitted non-linearity.
Normality: sm.qqplot(res.resid, line="s"), and the Jarque-Bera test in the summary. With large samples, mild non-normality barely matters for coefficient inference.
Influence: infl = res.get_influence() gives leverage (hat_matrix_diag), studentized residuals and Cook's distance (infl.cooks_distance[0]). A common rule of thumb flags points with Cook's distance above 4/n. sm.graphics.influence_plot(res) shows them. Investigate flagged points; don't delete them automatically.
Multicollinearity: variance_inflation_factor(exog, i) from statsmodels.stats.outliers_influence; values above roughly 5-10 mean a coefficient's variance is inflated by correlation with other predictors. The summary's condition number warns about the same thing.
Autocorrelation (time-ordered data): the Durbin-Watson statistic and acorr_ljungbox(res.resid).
Example
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.outliers_influence import variance_inflation_factor
rng = np.random.default_rng(7)
df = pd.DataFrame({"x1": rng.uniform(1, 10, 200)})
df["x2"] = df["x1"] * 0.9 + rng.normal(0, 0.5, 200) # nearly collinear with x1
df["y"] = 3 + 2 * df["x1"] + rng.normal(0, df["x1"]) # noise grows with x1
res = smf.ols("y ~ x1 + x2", data=df).fit()
lm_stat, lm_p, f_stat, f_p = het_breuschpagan(res.resid, res.model.exog)
print(f"{lm_p:.1e}") # 3.8e-07 (heteroscedasticity detected)
exog = res.model.exog
print([round(float(variance_inflation_factor(exog, i)), 1) for i in (1, 2)]) # [24.9, 24.9]
cooks = res.get_influence().cooks_distance[0]
print((cooks > 4 / len(df)).sum()) # 12 (points worth a look)Remember
- Plot residuals vs fitted values first
- Test heteroscedasticity with Breusch-Pagan; check VIF for collinearity
- Use Cook's distance to find influential observations
Watch out for
- Deleting every flagged influential point instead of understanding why it is unusual.