From estimates to uncertainty
A coefficient is a single number from a single sample. Its standard error says how much it would move across repeated samples (Lesson 2.2). From it we build:
- t-statistic = estimate / SE: how many standard errors away from zero the estimate is;
- p-value: how often you would see a t-statistic this large if the true effect were zero;
- 95% confidence interval ≈ estimate ± 1.96 × SE.
!pip install -q numpy pandas statsmodels matplotlib seaborn
import numpy as np, pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt, seaborn as sns
rng = np.random.default_rng(42)
sns.set_theme(style="whitegrid")
What a p-value promises: a placebo simulation
If a programme has no effect, a test at the 5% level should wrongly declare significance in 5% of studies. We check that by running 1,000 placebo studies.
def placebo_study(n=200):
d = rng.binomial(1, 0.5, n) # a 'programme' with zero effect
y = 50 + rng.normal(0, 10, n)
p = smf.ols("y ~ d", data=pd.DataFrame({"y": y, "d": d})).fit().pvalues["d"]
return p < 0.05
false_positive_rate = np.mean([placebo_study() for _ in range(1_000)])
print(f"Share of placebo studies with p < 0.05: {false_positive_rate:.3f}")
About 5%, as promised, when the assumptions behind the standard errors hold. Next we break one.
Heteroskedasticity and robust standard errors
Classical standard errors assume the noise has the same variance for everyone (homoskedasticity). Wages are a counterexample: they vary far more among people with postgraduate degrees than among early school-leavers.
N = 1_000
education = rng.uniform(8, 20, N)
errors = rng.normal(0, 150 * (education - 7) ** 2) # noise SD grows fast with education
wages = 25_000 + 3_000 * education + errors
df = pd.DataFrame({"education": education, "wages": wages})
classic = smf.ols("wages ~ education", data=df).fit()
robust = smf.ols("wages ~ education", data=df).fit(cov_type="HC3")
print(f"Coefficient (identical in both): {classic.params['education']:,.0f}")
print(f"Classical SE: {classic.bse['education']:,.1f} | Robust HC3 SE: {robust.bse['education']:,.1f}")
print(f"Classical SEs understate the uncertainty by {100 * (1 - classic.bse['education'] / robust.bse['education']):.0f}%")
plt.figure(figsize=(9, 4))
plt.scatter(classic.fittedvalues, classic.resid, alpha=.4)
plt.axhline(0, color="red"); plt.title("Residuals fan out: heteroskedasticity")
plt.xlabel("Fitted wage"); plt.ylabel("Residual"); plt.show()
The coefficient is unchanged: heteroskedasticity does not bias OLS. It biases the standard errors: here the classical ones are too small, so confidence intervals are too narrow and false positives too frequent. The Breusch–Pagan test makes the diagnosis formal:
from statsmodels.stats.diagnostic import het_breuschpagan
lm, lm_p, _, _ = het_breuschpagan(classic.resid, classic.model.exog)
print(f"Breusch-Pagan p-value: {lm_p:.2g} (small -> heteroskedasticity)")
Default to heteroskedasticity-robust standard errors (cov_type="HC1" or "HC3") in applied work. They
are valid whether or not the variance is constant, and they cost almost nothing.
Clustered data: when observations are not independent
Policies are often assigned to groups: schools, villages, provinces. Students in the same school share teachers and shocks, so their outcomes are correlated. Treating them as independent overstates how much information you have, and standard errors come out far too small.
def clustered_study(n_schools=40, pupils=25, cluster_se=False):
school = np.repeat(np.arange(n_schools), pupils)
treated_school = rng.permutation(np.repeat([0, 1], n_schools // 2))[school] # half the schools
school_shock = rng.normal(0, 5, n_schools)[school] # shared within a school
y = 50 + school_shock + rng.normal(0, 5, school.size) # NO true effect
d = pd.DataFrame({"y": y, "d": treated_school, "school": school})
if cluster_se:
fit = smf.ols("y ~ d", data=d).fit(cov_type="cluster", cov_kwds={"groups": d["school"]}, use_t=True)
else:
fit = smf.ols("y ~ d", data=d).fit(cov_type="HC1")
return fit.pvalues["d"] < 0.05
for flag in (False, True):
rate = np.mean([clustered_study(cluster_se=flag) for _ in range(500)])
print(f"Cluster-robust SE = {flag!s:5}: false-positive rate = {rate:.3f}")
Ignoring the clustering, a programme with no effect looks significant in a large share of studies. Clustering at the level of assignment brings the rate back near 5%. Rule: cluster at the level at which treatment was assigned. You will apply it in Modules 4 and 11.
Exercise
With only 8 schools, compare cluster-robust standard errors that use normal critical values
(use_t=False) with those that use the small-sample t correction (use_t=True, as above). Which keeps
the false-positive rate near 5%? What does that tell you about evaluating a policy rolled out in only a
handful of provinces?
Show a solution
def few_clusters(use_t, n_schools=8, pupils=25):
school = np.repeat(np.arange(n_schools), pupils)
d = pd.DataFrame({
"y": 50 + rng.normal(0, 5, n_schools)[school] + rng.normal(0, 5, school.size),
"d": rng.permutation(np.repeat([0, 1], n_schools // 2))[school],
"school": school,
})
fit = smf.ols("y ~ d", data=d).fit(cov_type="cluster", cov_kwds={"groups": d["school"]}, use_t=use_t)
return fit.pvalues["d"] < 0.05
for use_t in (False, True):
print(f"use_t={use_t!s:5}: false-positive rate = {np.mean([few_clusters(use_t) for _ in range(500)]):.3f}")
# With few clusters, normal critical values over-reject. The t correction helps, but with so few
# clusters inference stays fragile: consider randomisation inference (Module 4) or synthetic control (Module 12).
Checklist and recap
- I can read SEs, t-statistics, p-values and confidence intervals.
- I can check by simulation that a test has the right false-positive rate.
- I use robust SEs by default and cluster at the level of treatment assignment.