import matplotlib
matplotlib.use("Agg") # headless: render to buffers, never to a screen
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyfixest as pf
rng = np.random.default_rng(42) # one generator, pinned, drives the whole notebook
plt.rcParams["figure.figsize"] = (9, 5)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3
pd.set_option("display.width", 100)
pd.set_option("display.max_columns", 20)
print("numpy ", np.__version__)
print("pandas ", pd.__version__)
print("pyfixest", pf.__version__)nb7.5 — First-Look Regressions: Frozen Until the PAP Is Filed
Week 7, Chapter 7.5 companion. Maya runs the first estimate of her primary specification — once.
This is the most disciplined notebook in the whole camp, and the discipline is the point. By now Maya has done everything Week 7 asked: she generated her idea (7.1), pulled and documented her HMDA-style data (7.2), filed a pre-analysis plan as a tagged commit (7.3), built a clean analysis dataset (7.4), and written the identification memo whose one assumption her headline number rests on (7.5). Only now may she look at the result.
The notebook embodies one rule, the rule that the whole week was built to enforce:
You run only your pre-registered primary specification, exactly as the §4 spec line states, once.
Not the version that gives the prettiest star. Not the clustering level that tightens the t-stat. The one specification she committed to before she could see a coefficient. Chapter 7.3 showed why: with the data already in hand she could peek, so the only thing standing between her and the garden of forking paths is a plan she froze in advance. This notebook is the technological enforcement of that freeze — a gate that refuses to run the regression until PAP_FILED is True.
We do this on a synthetic panel with a known true effect, exactly as the Week 4 DiD notebooks did, for one reason: when you build the world yourself, you know the right answer, so you can check that an honestly-run primary specification — clustered standard errors and all — recovers the planted effect. On real data you never get that luxury; here we use it to prove the machinery is honest.
1. The analysis panel (schema like the nb7.4 output)
We load the analysis dataset Maya built in nb7.4. Here we generate it in-notebook so the notebook is self-contained and reproducible, but the schema deliberately matches what a real nb7.4 build emits: one row per county–year, with a treatment cohort, the outcome, the named controls, and the cluster key.
The setting is the fair-lending DiD of Chapter 7.5’s worked memo. Some states adopt a fair-lending examination program in a known year; never-adopting states are the clean controls. The outcome is the county-level minority–white mortgage-denial gap in percentage points (denial rate for minority applicants minus the rate for white applicants). The program, where active, narrows the gap — so the planted true effect is negative.
We write the data-generating process (DGP) as the two-way fixed-effects (TWFE) model run forward, in the notation of CONVENTIONS §3:
\texttt{gap}_{it} = \underbrace{\alpha_i}_{\text{county level}} + \underbrace{\lambda_t}_{\text{year shock}} + \underbrace{\tau \cdot \texttt{post\_treat}_{it}}_{\text{the planted effect}} + \mathbf{x}_{it}'\boldsymbol{\gamma} + \varepsilon_{it},
where \texttt{post\_treat}_{it}=1 once county i’s state has an active program, \tau is the true ATT we plant, \alpha_i is the county fixed effect, \lambda_t the year fixed effect, and \mathbf{x}_{it} the applicant-composition controls. The county and year fixed effects absorb \alpha_i and \lambda_t; if our spec is honest, the coefficient on the treatment indicator should land on \tau.
# --- Planted truth and panel dimensions ---------------------------------
TAU_TRUE = -1.80 # the programs NARROW the denial gap by 1.8 pp (the effect we must recover)
YEARS = list(range(2014, 2024)) # 10 years
N_STATES = 30 # 30 states -> 30 clusters (treatment varies at the STATE level)
COUNTIES_PER_STATE = 8 # 240 counties, ~2,400 county-year rows
GAMMA_INCOME = -0.9 # higher applicant income -> smaller raw gap
GAMMA_LTV = 1.2 # higher loan-to-value -> larger raw gap
SIGMA = 0.8 # idiosyncratic county-year noise (sd)
# Half the states adopt; adoption YEAR is staggered across a couple of cohorts.
adopt_cohort = {} # state -> year program takes effect (np.nan = never)
for s in range(N_STATES):
if s < 8: adopt_cohort[s] = 2018
elif s < 15: adopt_cohort[s] = 2020
else: adopt_cohort[s] = np.nan # never-adopters = clean controls
# State-level mean shift in the baseline gap (some states just have wider gaps).
state_level = {s: rng.normal(0.0, 1.0) for s in range(N_STATES)}
# Year shocks common to all counties (national mortgage-market conditions).
year_shock = {t: 0.10 * (t - 2014) + rng.normal(0.0, 0.3) for t in YEARS}
rows = []
for s in range(N_STATES):
G = adopt_cohort[s]
for c in range(COUNTIES_PER_STATE):
county_id = s * COUNTIES_PER_STATE + c
alpha_i = state_level[s] + rng.normal(0.0, 0.7) # county permanent level
for t in YEARS:
# Named controls (applicant composition), drawn per county-year.
income = rng.normal(0.0, 1.0) # standardized log applicant income
ltv = rng.normal(0.0, 1.0) # standardized loan-to-value
post_treat = 1 if (not np.isnan(G) and t >= G) else 0
gap = (alpha_i
+ year_shock[t]
+ TAU_TRUE * post_treat
+ GAMMA_INCOME * income
+ GAMMA_LTV * ltv
+ rng.normal(0.0, SIGMA))
rows.append({
"county": county_id,
"state": s,
"year": t,
"cohort_G": G, # year the state's program took effect (NaN = never)
"post_treat": post_treat, # 1 once this county's state has an active program
"income": income,
"ltv": ltv,
"gap": gap, # OUTCOME: minority-white denial gap, pp
})
panel = pd.DataFrame(rows)
panel["adopter"] = panel["cohort_G"].notna().astype(int)
print(f"county-year rows : {len(panel):,}")
print(f"states (clusters): {panel['state'].nunique()}")
print(f"counties : {panel['county'].nunique()}")
print(f"adopting states : {panel.loc[panel['adopter']==1, 'state'].nunique()}"
f" never-adopters: {panel.loc[panel['adopter']==0, 'state'].nunique()}")
print(f"PLANTED true ATT : {TAU_TRUE:+.2f} pp (kept hidden from the spec; used only to grade it)")
panel.head()A quick, non-confirmatory sanity look at the schema — the kind of glance at structure that does not touch any coefficient. We confirm the panel is balanced and that the treatment switches on at the right cohort years. Looking at the design (who is treated, when) is always allowed; looking at the outcome’s relationship to treatment before the gate is exactly what the freeze forbids.
# Balance + treatment-timing check (design only -- no outcome regression here).
counts = panel.groupby("year").size()
print("rows per year (should be constant -> balanced panel):")
print(counts.to_string())
print("\ntreated county-year cells by cohort:")
timing = (panel[panel["adopter"] == 1]
.groupby("cohort_G")["post_treat"].agg(["mean", "size"])
.rename(columns={"mean": "share_post", "size": "county_years"}))
print(timing.to_string())
assert panel.groupby(["county", "year"]).size().max() == 1, "panel must be one row per county-year"
print("\nSchema OK: balanced county-year panel, staggered adoption, never-adopter controls present.")2. The PAP-FILED gate — frozen until the plan is a tagged commit> ## PAP-FILED GATE>> The primary specification below does not run until the pre-analysis plan is filed.>> In Lab 7 the PAP and the identification memo are committed and tagged (git tag pap-filed). That> tag is the time-stamp that proves the spec was fixed before any coefficient was seen (Ch 7.3.3). The> flag PAP_FILED is the notebook’s stand-in for that tag. While it is False, the analysis cell below> refuses to run and raises an error. You flip it to True only once the plan is genuinely filed —> and once it is filed, you may not change the spec to chase a nicer number.This is the freeze made mechanical. Chapter 7.3 argued that “don’t peek” is a request for superhumanwillpower; a gate that physically halts the regression turns it into a fact. The assertion is small,but it is the entire ethic of the week compressed into one line: commit first, then test.
# Flip to True ONLY after the PAP + identification memo are filed as a tagged commit (Lab 7).
# Setting it True here represents Maya having filed `git tag pap-filed`.
PAP_FILED = True
# The pre-registered primary spec, frozen in advance, transcribed verbatim from the PAP's section 4.
# Changing any field after filing would break the contract (and would be logged as a deviation).
PRIMARY_SPEC = {
"outcome": "gap", # county-year minority-white denial gap (pp)
"treatment": "post_treat", # 1 once the county's state has an active exam program
"controls": ["income", "ltv"], # applicant-composition covariates (named in the PAP)
"fixed_effects": ["county", "year"], # county and year FE
"cluster": "state", # CLUSTER at the level treatment varies -- chosen in advance
"test": "two-sided", # H1 effect direction is negative; report two-sided to be strict
}
print("PAP_FILED =", PAP_FILED)
for k, v in PRIMARY_SPEC.items():
print(f" {k:14s}: {v}")3. Run the pre-registered primary specification — once
Here is the spec, in CONVENTIONS §4 form, exactly as the PAP’s §4 line states it:
Outcome:
gap(county-year minority–white denial gap, pp). Treatment / key regressor:post_treat(= 1 once the county’s state has an active fair-lending examination program). Controls:income,ltv(applicant-composition covariates). Fixed effects: county and year. Clustering: by state (the level at which treatment varies). Sample: the full balanced county-year panel, 2014–2023. Identifying assumption (one sentence): absent the examination programs, adopting and never-adopting states’ county-level denial gaps would have moved in parallel.
In pyfixest syntax that is one line: feols("gap ~ post_treat + income + ltv | county + year", ...) with the standard errors pre-specified as clustered by state via vcov={"CRV1": "state"}. The clustering choice is not a knob we turn after seeing the t-stat — Chapter 7.3 called that “p-hacking on the standard error.” It was nailed down in the PAP, so we read off whatever it gives.
The analysis cell asserts PAP_FILED before doing anything. If the plan were not filed, the cell would raise and refuse to estimate.
# ---- THE GATE: refuse to estimate until the PAP is filed --------------------
assert PAP_FILED, (
"PAP NOT FILED. The primary specification is frozen until the pre-analysis plan and "
"identification memo are committed and tagged (git tag pap-filed). Flip PAP_FILED to True "
"ONLY after filing -- and do not change the spec to chase a result."
)
# ---- Build the formula straight from the frozen PRIMARY_SPEC (no hand-editing) ----
rhs_controls = " + ".join(PRIMARY_SPEC["controls"])
fe = " + ".join(PRIMARY_SPEC["fixed_effects"])
formula = f'{PRIMARY_SPEC["outcome"]} ~ {PRIMARY_SPEC["treatment"]} + {rhs_controls} | {fe}'
print("Estimating (verbatim from PAP S4):")
print(" ", formula)
print(" SEs: clustered by", PRIMARY_SPEC["cluster"], "(pre-specified)\n")
# ---- Run it ONCE, with the PRE-SPECIFIED clustered standard errors -------------
primary = pf.feols(formula, data=panel, vcov={"CRV1": PRIMARY_SPEC["cluster"]})
print(primary.tidy().round(4).to_string())4. Report the headline number — with the pre-specified clustered SEs
Pull the one coefficient H1 is about — the coefficient on post_treat — together with its clustered standard error, t-statistic, p-value, and 95% confidence interval. This is the headline number the identification memo defends, reported exactly as promised. Because we built the world, we can also check the honest run against the planted truth: the confidence interval should comfortably cover \tau.
key = PRIMARY_SPEC["treatment"]
beta = float(primary.coef()[key])
se = float(primary.se()[key])
tval = float(primary.tstat()[key])
pval = float(primary.pvalue()[key])
ci = primary.confint().loc[key]
lo, hi = float(ci.iloc[0]), float(ci.iloc[1])
print("=" * 60)
print(" PRIMARY SPECIFICATION -- CONFIRMATORY RESULT (H1)")
print("=" * 60)
print(f" coefficient on {key:<12s}: {beta:+.4f} pp")
print(f" clustered SE (by state) : {se:.4f}")
print(f" t-statistic : {tval:+.3f}")
print(f" p-value (two-sided) : {pval:.4g}")
print(f" 95% CI (clustered) : [{lo:+.4f}, {hi:+.4f}]")
print("-" * 60)
print(f" planted true ATT : {TAU_TRUE:+.4f} pp")
covered = lo <= TAU_TRUE <= hi
print(f" CI covers the truth? : {covered}")
print("=" * 60)
# Honest-machinery checks: the pre-registered spec, run once, recovers the planted effect.
assert covered, "clustered 95% CI failed to cover the planted ATT -- spec is not recovering the truth"
assert beta < 0, "expected a NEGATIVE effect (programs narrow the gap); sign is wrong"
assert abs(beta - TAU_TRUE) < 0.30, "point estimate drifted too far from the planted truth"
print("\nRecovered: the frozen primary spec lands on the planted effect, with clustered SEs.")A coefficient plot of the single confirmatory estimate — point estimate and clustered 95% interval — against the planted truth. One number, reported once. The plot is the whole confirmatory result; there is nothing else on this canvas, because there is nothing else in the primary analysis.
fig, ax = plt.subplots(figsize=(7.5, 3.2))
ax.errorbar([beta], [0], xerr=[[beta - lo], [hi - beta]],
fmt="o", color="#3b6ea5", capsize=6, markersize=9,
label="primary estimate (clustered 95% CI)")
ax.axvline(TAU_TRUE, color="green", ls="--", lw=1.5, label=f"planted truth ({TAU_TRUE:+.2f})")
ax.axvline(0.0, color="grey", ls=":", lw=1, label="no effect")
ax.set_yticks([])
ax.set_xlabel("effect on minority-white denial gap (pp)")
ax.set_title("nb7.5 -- the ONE confirmatory number (primary spec, clustered SEs)")
ax.legend(loc="upper left", fontsize=9)
fig.tight_layout()
fig.savefig("nb75_primary_coef.png", dpi=110)
plt.close(fig)
print("saved figure: nb75_primary_coef.png")
print("That single estimate -- and only that estimate -- is the confirmatory result of the project.")—## The confirmatory analysis is now COMPLETE.The headline number above is the entire pre-registered result. Everything below this line isEXPLORATORY — generated after seeing the data, valuable as a source of new hypotheses, but not aconfirmatory test of this project’s H1. It is fenced off and labeled so it can never be mistaken forthe registered result. This is the confirmatory/exploratory partition of Chapter 7.3.1, drawn on thepage where a reader can see exactly which results were tests and which were searches.—## 5. EXPLORATORY — not confirmatory> EXPLORATORY ZONE. Nothing below was pre-registered. These are post-hoc looks: they generate> hypotheses for the next study, they do not test this one. Reporting any of them as if it were the> primary result would be exactly the laundering of exploration as confirmation that the PAP exists to> prevent.To make the danger concrete, we walk a few forks Maya might be tempted by after seeing her headline —the same “garden of forking paths” from Chapter 7.3.1. Each is a perfectly defensible choice inisolation; each is a different number; and choosing among them by which gives the prettiest star isp-hacking. We run them only to show the forks, all clearly stamped exploratory.
# EXPLORATORY: a handful of alternative specifications Maya did NOT pre-register.
# Shown to illustrate the garden of forking paths -- NOT to be reported as the result.
explore = {}
# Fork A: drop the controls.
explore["no controls"] = pf.feols("gap ~ post_treat | county + year",
data=panel, vcov={"CRV1": "state"})
# Fork B: cluster by county instead of state (a DIFFERENT, non-registered SE choice).
explore["cluster by county"] = pf.feols(
"gap ~ post_treat + income + ltv | county + year",
data=panel, vcov={"CRV1": "county"})
# Fork C: drop 2020 (a post-hoc sample exclusion).
explore["drop 2020"] = pf.feols(
"gap ~ post_treat + income + ltv | county + year",
data=panel[panel["year"] != 2020].copy(), vcov={"CRV1": "state"})
# Fork D: add a control on the causal pathway (over-controlling -- a tempting but wrong fork).
panel_x = panel.copy()
panel_x["loan_size"] = 0.5 * panel_x["ltv"] + 0.4 * panel_x["post_treat"] + rng.normal(0, 1, len(panel_x))
explore["over-control (bad)"] = pf.feols(
"gap ~ post_treat + income + ltv + loan_size | county + year",
data=panel_x, vcov={"CRV1": "state"})
forks = pd.DataFrame({
name: {
"beta": float(m.coef()["post_treat"]),
"se": float(m.se()["post_treat"]),
"pval": float(m.pvalue()["post_treat"]),
} for name, m in explore.items()
}).T
# Put the registered result at the top for contrast.
forks.loc["** PRIMARY (registered) **"] = [beta, se, pval]
forks = forks.reindex(["** PRIMARY (registered) **"] + list(explore.keys()))
print("EXPLORATORY forks vs. the one registered spec (all clearly labeled):\n")
print(forks.round(4).to_string())
print("\nEvery row is a defensible choice in isolation; only the top row is the confirmatory test.")
print("Picking a row by which gives the smallest p-value is the garden of forking paths (Ch 7.3.1).")Plot the forks side by side against the planted truth. They scatter — that scatter is the researcher degrees of freedom Chapter 7.3 warned about. Note especially the over-control fork: by conditioning on a variable that sits on the treatment’s causal pathway, it absorbs part of the very effect we are hunting and biases the estimate toward zero (the Mentor 4 over-controlling trap). The point of showing it is not to choose among these; it is to make visible why the choice had to be frozen before this picture existed.
fig, ax = plt.subplots(figsize=(8.5, 4.2))
ypos = np.arange(len(forks))[::-1]
colors = ["#b5651d" if "PRIMARY" not in name else "#3b6ea5" for name in forks.index]
ax.errorbar(forks["beta"], ypos, xerr=1.96 * forks["se"],
fmt="o", capsize=5, markersize=8, ecolor="grey", linestyle="none",
mfc="none")
for y, (name, row), col in zip(ypos, forks.iterrows(), colors):
ax.plot(row["beta"], y, "o", color=col, markersize=9)
ax.axvline(TAU_TRUE, color="green", ls="--", lw=1.5, label=f"planted truth ({TAU_TRUE:+.2f})")
ax.axvline(0.0, color="grey", ls=":", lw=1)
ax.set_yticks(ypos)
ax.set_yticklabels(forks.index, fontsize=9)
ax.set_xlabel("effect on denial gap (pp), 95% CI")
ax.set_title("EXPLORATORY: the garden of forking paths (blue = registered; orange = post-hoc)")
ax.legend(loc="lower right", fontsize=9)
fig.tight_layout()
fig.savefig("nb75_forks.png", dpi=110)
plt.close(fig)
print("saved figure: nb75_forks.png")6. Proving the gate works
Before we leave the discipline, let us show the gate has teeth. If PAP_FILED were False, the analysis cell’s assert would fire and refuse to estimate. We demonstrate it here in a contained way — catching the error so the notebook still runs end-to-end — so you can see the freeze is real and not decorative.
def run_primary(pap_filed_flag):
"""The gated analysis, factored out so we can demonstrate the freeze both ways."""
assert pap_filed_flag, (
"PAP NOT FILED -- primary specification is frozen until `git tag pap-filed`."
)
return pf.feols(formula, data=panel, vcov={"CRV1": PRIMARY_SPEC["cluster"]})
# With the PAP NOT filed, the gate must block the run:
try:
run_primary(pap_filed_flag=False)
print("UNEXPECTED: the gate did not fire.")
except AssertionError as e:
print("GATE FIRED (as it should) when PAP_FILED is False:")
print(" ", str(e))
# With the PAP filed, it runs and reproduces the same headline coefficient:
m_again = run_primary(pap_filed_flag=True)
print(f"\nGate open (PAP filed) -> coefficient reproduces: {float(m_again.coef()['post_treat']):+.4f} pp")
assert abs(float(m_again.coef()["post_treat"]) - beta) < 1e-9, "gated run must match the headline exactly"
print("The freeze is mechanical: no plan, no regression.")7. The bridge to Week 8 — robustness comes next, on the frozen result
You now have your headline number, and you have it honestly: one pre-registered specification, run once, with the standard errors you chose in advance. What you do not do is start swapping specifications to see which looks best — that was the whole point of the freeze.
So what do you do next? Week 8 is robustness, run on this frozen result. The threats-and-responses table you built in the identification memo (Ch 7.5) becomes the skeleton of your robustness section: Chapter 8.2 is, almost literally, the instruction “go run column 3 for every row of your table” — the placebo tests, the alternative standard errors, the sample and bandwidth perturbations you promised in your PAP’s §5 to report whichever way they come out. The exploratory forks in Section 5 above are not that: they are unregistered searches, useful only for generating next-study hypotheses. A planned robustness check is a fork you committed to walking in advance and to reporting honestly; that is what makes it confirmatory-grade rather than fishing.
The thread of the whole week, one last time: the PAP says you did not fish, and the identification memo says the pond is the right pond. This notebook ran the first cast — once, under the freeze. Week 8 is where you stress-test the fish you caught, on the result you are no longer allowed to change.
Your Turn
Move the truth. Change
TAU_TRUEto0.0(a true null) and re-run. The primary spec should now produce a confidence interval that covers zero — the pre-registered honest null Chapter 7.3.4 said is a genuine contribution, not a failure. Notice you did not get to go fishing for a star.Strain the clusters. In
adopt_cohort, make only 2 or 3 states ever adopt. Re-run the primary spec. Watch the clustered standard error balloon — with few treated clusters, cluster-robust inference is strained (Ch 4.1’s few-treated-clusters warning). What would your PAP’s §5 placebo permutation say here, and why is it the more honest summary?Audit a fork. Pick one EXPLORATORY fork from Section 5 whose number you liked better than the registered one. Write the single sentence you would have to put in a deviation log (Ch 7.3.3) to use it instead — including the honest clause: was the change prompted by seeing the coefficient? If yes, it is exploratory, full stop.