TopolearnWrong standard errors produce falsely confident conclusions even when the coefficients are fine.
When residual variance is not constant, or observations are grouped (students in schools, repeated measures on customers), OLS coefficients are still unbiased, but the default standard errors are wrong, usually too small.
The fix is to keep the estimates and change the covariance estimator at fit time with cov_type:
fit(cov_type="HC3"): heteroscedasticity-consistent ("robust", "White") standard errors. HC3 behaves best in small samples; HC0-HC2 also exist.fit(cov_type="cluster", cov_kwds={"groups": df["school_id"]}): cluster-robust standard errors, allowing arbitrary correlation within each group. You need a reasonable number of clusters (dozens, not five).fit(cov_type="HAC", cov_kwds={"maxlags": 4}): Newey-West standard errors for time series with autocorrelation.
Coefficients are identical across these options; only bse, pvalues, conf_int() and the tests change. The summary notes which covariance type was used.
cov_type also works for GLM and discrete models such as logit. For clustered data where you also want to model the within-group correlation, mixed models or GEE are the next step.
In practice, robust standard errors are a cheap default for cross-sectional data, and clustering at the level where treatment or sampling happened is essential in experiments and panel data.
Example
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
rng = np.random.default_rng(8)
schools = np.repeat(np.arange(40), 25) # 40 schools x 25 students
tutoring = rng.integers(0, 2, 40)[schools] # assigned per school
school_effect = rng.normal(0, 5, 40)[schools]
df = pd.DataFrame({"school": schools, "tutoring": tutoring,
"score": 60 + 3 * tutoring + school_effect + rng.normal(0, 5, 1000)})
naive = smf.ols("score ~ tutoring", data=df).fit()
robust = smf.ols("score ~ tutoring", data=df).fit(cov_type="HC3")
clustered = smf.ols("score ~ tutoring", data=df).fit(
cov_type="cluster", cov_kwds={"groups": df["school"]})
for r in (naive, robust, clustered):
print(round(r.params["tutoring"], 2), round(r.bse["tutoring"], 2), round(r.pvalues["tutoring"], 4))
# 2.8 0.49 0.0 (same coefficient everywhere...)
# 2.8 0.49 0.0
# 2.8 1.91 0.1424 (...but the honest uncertainty is 4x larger)Remember
cov_typechanges standard errors, never coefficients- HC3 for heteroscedasticity, cluster for grouped data, HAC for autocorrelation
- Cluster at the level where treatment or sampling happened
Watch out for
- Using default standard errors when the treatment was assigned to groups, not individuals.