import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg") # non-interactive backend: render to buffers, never pop a window
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(42) # one seed for the whole notebook -> reproducible
# Maya's design constants, in the UNITS OF HER OUTCOME.
# `denied` is a 0/1 indicator, so effects are in PERCENTAGE POINTS (pp) of denial probability.
ALPHA = 0.05 # significance level committed to in the PAP
TARGET_POWER = 0.80 # the conventional 80% power target
MDE_PP = 2.0 # smallest gap (in pp) Maya says would matter to a real borrower
# These two numbers come FROM THE EXPLORATION SET ONLY (the 30% she is allowed to peek at).
# For a 0/1 outcome with overall denial rate ~12%, the residual SD is sqrt(p(1-p)) ~ 0.325.
P_BAR = 0.12 # baseline denial rate in the exploration set
RESID_SD = np.sqrt(P_BAR * (1 - P_BAR)) # ~0.325 (in probability units)
MINORITY_SHARE = 0.30 # fraction of applications flagged `minority`
N_REALISTIC = 8000 # rows she expects in the 70% confirmation set
print(f"alpha={ALPHA}, target power={TARGET_POWER}")
print(f"baseline denial rate p_bar = {P_BAR} -> residual SD ~ {RESID_SD:.4f}")
print(f"minority share = {MINORITY_SHARE}, realistic confirmation N = {N_REALISTIC}")
print(f"effect Maya cares about (MDE target): {MDE_PP} percentage points = {MDE_PP/100:.4f} in prob units")Notebook 7.3 — PAP Companion: A Power Calculation for Your Design
Maya has written the prose of her pre-analysis plan in Chapter 7.3 — directional hypotheses, one primary specification, a multiple-testing plan, falsification conditions. Everything except the one number a referee will ask about first: is your study even capable of detecting the effect you care about? A pre-registered design that is doomed to be underpowered is not a virtue; it is a confirmatory run wasted on a question your data were never large enough to answer. The honest move is to find that out now, on the exploration set, before you spend your one shot at the confirmation data.
This notebook is the direct sequel to nb1.5 (power curves and the t-test). We reuse that same machinery — true effect sits some number of standard errors above the null; the test fires past a critical value; power is the tail area beyond z_crit - ncp — and aim it at the five jobs the PAP companion has to do:
- an analytic power calculator (effect, SD, \alpha, N \to power) and its inverse (N for 80% power), for a two-group mean comparison and for a regression coefficient (residual SD and regressor variation);
- a Monte Carlo power simulator for your actual planned spec — simulate under the alternative, run the spec, count rejections — and confirm analytic \approx simulated;
- a minimum detectable effect (MDE) calculator given the realistic N you actually have;
- a runnable Benjamini–Hochberg FDR routine (p-values \to which survive at level q), with the worked income-band family from §7.3.2;
- a PAP-template emitter that prints a filled short-PAP skeleton you can paste into your own plan.
Throughout, the running example is Maya’s HMDA fair-lending design: the conditional denial gap (denied on minority, plus controls and fixed effects), the FinTech interaction, and the subgroup family she quarantines under FDR.
1. Analytic power for a two-group mean comparison
Start with the simplest version of Maya’s question, stripped of fixed effects: split the applicants into minority (n_1) and non-minority (n_0), and compare mean denial rates. The difference in means \hat\Delta=\bar y_1-\bar y_0 has standard error
\operatorname{se}(\hat\Delta)=\sigma\sqrt{\tfrac{1}{n_1}+\tfrac{1}{n_0}},
where \sigma is the (pooled) residual SD of the outcome. This is exactly the nb1.5 logic with a two-sample standard error in place of the one-sample \sigma/\sqrt N. The true effect \Delta sits \text{ncp}=\Delta/\operatorname{se} standard errors away from zero, the test fires past the critical value z_{\text{crit}}, and
\text{power}=\Pr\big(Z > z_{\text{crit}}-\text{ncp}\big)\quad(\text{one-sided}),\qquad \text{power}=\Pr\big(Z > z_{\alpha/2}-\text{ncp}\big)+\Pr\big(Z < -z_{\alpha/2}-\text{ncp}\big)\quad(\text{two-sided}).
We write one function and own every step, mirroring simulate_rejection_rate from nb1.5 but on the analytic side.
def power_two_group(delta, sigma, n1, n0, alpha=ALPHA, alternative="greater"):
"""Analytic power for a two-group difference in means (normal approximation).
delta : true difference in means (group1 - group0), in outcome units.
sigma : residual SD of the outcome (pooled).
n1,n0 : group sizes.
alternative: 'greater'/'less' -> one-sided; 'two-sided' -> two-sided.
Returns power in [0,1].
"""
se = sigma * np.sqrt(1.0 / n1 + 1.0 / n0)
ncp = delta / se # how many SEs the true effect sits from 0
if alternative == "two-sided":
z = stats.norm.ppf(1 - alpha / 2)
return stats.norm.sf(z - ncp) + stats.norm.cdf(-z - ncp)
elif alternative == "greater":
z = stats.norm.ppf(1 - alpha)
return stats.norm.sf(z - ncp)
elif alternative == "less":
z = stats.norm.ppf(1 - alpha)
return stats.norm.sf(z + ncp) # mirror image
else:
raise ValueError("alternative must be 'greater', 'less', or 'two-sided'")
# Maya's planned confirmation N, split by her minority share, one-sided (H1: gap > 0).
n1 = int(round(N_REALISTIC * MINORITY_SHARE)) # minority applications
n0 = N_REALISTIC - n1 # non-minority applications
delta = MDE_PP / 100.0 # 2 pp = 0.02 in probability units
pw = power_two_group(delta, RESID_SD, n1, n0, alpha=ALPHA, alternative="greater")
print(f"Groups: n1(minority)={n1}, n0={n0}")
print(f"True gap = {MDE_PP} pp ({delta:.3f}). Residual SD = {RESID_SD:.4f}")
print(f"Analytic ONE-sided power at N={N_REALISTIC}: {pw:.4f}")
print(f"Sanity: a zero effect must give power == alpha:",
f"{power_two_group(0.0, RESID_SD, n1, n0):.4f} (target {ALPHA})")2. Inverting it: the N for 80% power
Power is a knob you read; N is the lever you pull. The PAP question is the inverse: how big must the sample be to reach 80% power for the effect I care about? For the one-sided case the power equation
\text{power}=\Phi\!\Big(\frac{\Delta}{\sigma\sqrt{1/n_1+1/n_0}}-z_\alpha\Big)
inverts cleanly. Require power =1-\beta, so the argument of \Phi must equal z_{1-\beta}:
\frac{\Delta}{\sigma\sqrt{1/n_1+1/n_0}}-z_\alpha=z_{1-\beta} \;\Longrightarrow\; \frac{1}{n_1}+\frac{1}{n_0}=\Big(\frac{\Delta}{\sigma\,(z_\alpha+z_{1-\beta})}\Big)^{2}.
Holding the group fractions fixed (n_1=\pi N, n_0=(1-\pi)N), 1/n_1+1/n_0
=\frac{1}{N}\big(\tfrac{1}{\pi}+\tfrac{1}{1-\pi}\big), so N falls right out. We solve it in closed form, then verify it by feeding the answer back into power_two_group and checking we land on 0.80 — the same self-check discipline nb1.5 used (assert np.isclose).
def n_for_power_two_group(delta, sigma, alpha=ALPHA, power=TARGET_POWER,
frac1=0.5, alternative="greater"):
"""Smallest total N reaching `power` for a two-group mean comparison (normal approx).
frac1 : fraction of N assigned to group 1 (the rest to group 0).
Returns total N (rounded up to an integer).
"""
if alternative == "two-sided":
z_a = stats.norm.ppf(1 - alpha / 2)
else:
z_a = stats.norm.ppf(1 - alpha)
z_b = stats.norm.ppf(power)
# (1/n1 + 1/n0) required:
inv_sum_required = (delta / (sigma * (z_a + z_b))) ** 2
# with n1 = frac1*N, n0 = (1-frac1)*N: 1/n1 + 1/n0 = (1/frac1 + 1/(1-frac1)) / N
shape = 1.0 / frac1 + 1.0 / (1.0 - frac1)
N = shape / inv_sum_required
return int(np.ceil(N))
N_needed = n_for_power_two_group(delta, RESID_SD, alpha=ALPHA, power=TARGET_POWER,
frac1=MINORITY_SHARE, alternative="greater")
print(f"N needed for {int(TARGET_POWER*100)}% power at a {MDE_PP} pp gap (one-sided): {N_needed}")
# VERIFY the inversion: plug N_needed back in -> power should be ~0.80.
n1_chk = int(round(N_needed * MINORITY_SHARE))
n0_chk = N_needed - n1_chk
pw_chk = power_two_group(delta, RESID_SD, n1_chk, n0_chk, alpha=ALPHA, alternative="greater")
print(f"Check: power at N_needed = {pw_chk:.4f} (target {TARGET_POWER})")
assert abs(pw_chk - TARGET_POWER) < 0.01, "Inversion did not land on target power!"
# Two-sided costs more (1.96 vs 1.645): the price of not having had a directional prior.
N_two = n_for_power_two_group(delta, RESID_SD, frac1=MINORITY_SHARE, alternative="two-sided")
print(f"Two-sided N for {int(TARGET_POWER*100)}% power: {N_two} (more than one-sided {N_needed})")
print(f"\nMaya has ~{N_REALISTIC} rows; she needs ~{N_needed} (1-sided). "
f"{'ADEQUATE.' if N_REALISTIC >= N_needed else 'UNDERPOWERED -- see MDE in section 6.'}")3. Power for a regression coefficient (the spec she actually runs)
Maya’s primary specification is not a two-group t-test; it is a regression \texttt{denied}_i=\beta_0+\beta_1\texttt{minority}_i+\mathbf x_i'\boldsymbol\gamma+\dots, and she reads \beta_1 off the table. The power machinery is the same idea with one substitution. For OLS, the standard error of a slope is
\operatorname{se}(\hat\beta_1)=\frac{\sigma_\varepsilon}{\sqrt{N}\;\operatorname{sd}(\tilde x_1)},
where \sigma_\varepsilon is the residual SD (the noise left after the regression explains what it can) and \operatorname{sd}(\tilde x_1) is the SD of the key regressor after partialling out the other controls (its residual variation — what nb on Frisch–Waugh–Lovell taught: a coefficient’s precision depends on the part of x_1 the controls cannot explain). Two forces shrink the gap between this and the two-group formula: controls reduce \sigma_\varepsilon (good for power) but, if they are correlated with minority, they also reduce \operatorname{sd}(\tilde x_1) (bad for power — this is variance inflation). The PAP power calc must use the residualized regressor SD, not the raw one, or it will lie to you.
def power_regression_coef(beta, resid_sd, n, sd_x_resid, alpha=ALPHA, alternative="greater"):
"""Analytic power for a single OLS coefficient (normal approx).
beta : true coefficient (effect size, in outcome units per unit of x).
resid_sd : SD of the regression residual (sigma_eps).
n : sample size.
sd_x_resid : SD of the key regressor AFTER partialling out other controls.
"""
se = resid_sd / (np.sqrt(n) * sd_x_resid)
ncp = beta / se
if alternative == "two-sided":
z = stats.norm.ppf(1 - alpha / 2)
return stats.norm.sf(z - ncp) + stats.norm.cdf(-z - ncp)
elif alternative == "greater":
z = stats.norm.ppf(1 - alpha)
return stats.norm.sf(z - ncp)
elif alternative == "less":
z = stats.norm.ppf(1 - alpha)
return stats.norm.sf(z + ncp)
raise ValueError("bad alternative")
# For a binary regressor `minority` with share pi, its raw SD is sqrt(pi(1-pi)).
# After fixed effects + controls absorb some variation, the RESIDUAL sd is smaller; model
# that with a partial-R^2 of the controls explaining `minority` (variance inflation).
sd_x_raw = np.sqrt(MINORITY_SHARE * (1 - MINORITY_SHARE))
r2_controls_on_x = 0.20 # controls/FE explain 20% of minority's variance (assumed, from exploration set)
sd_x_resid = sd_x_raw * np.sqrt(1 - r2_controls_on_x) # residualized regressor SD
# Controls also shave the residual SD of the OUTCOME a little (say partial-R^2 ~ 0.10).
resid_sd_y = RESID_SD * np.sqrt(1 - 0.10)
pw_reg = power_regression_coef(delta, resid_sd_y, N_REALISTIC, sd_x_resid,
alpha=ALPHA, alternative="greater")
print(f"Raw regressor SD = {sd_x_raw:.4f}; residualized (after controls/FE) = {sd_x_resid:.4f}")
print(f"Outcome residual SD after controls = {resid_sd_y:.4f}")
print(f"Analytic power for beta_1 = {MDE_PP} pp at N={N_REALISTIC} (one-sided): {pw_reg:.4f}")
print("Note: residualizing the regressor RAISES its se and LOWERS power -- variance inflation.")4. A Monte Carlo power simulator for Maya’s actual spec
The analytic formulas are normal-approximation shortcuts. The honest check — and the thing a formula can never do for you — is to simulate the world under the alternative, run the exact spec you will run, and count how often it rejects. This is the nb1.5 move (simulate_rejection_rate) upgraded from a one-sample mean to Maya’s regression: draw a binary minority regressor, draw a binary denied outcome whose probability shifts by the true gap, fit OLS (a linear probability model — fine for power at these rates), and check whether \hat\beta_1 clears its critical value. We then confirm analytic \approx simulated, the central reassurance of the whole notebook.
def simulate_power_regression(beta, p_bar, n, minority_share, n_sims=2000,
alpha=ALPHA, alternative="greater", r=None):
"""Monte Carlo power for the minority coefficient in a linear probability model.
Data-generating process under H1:
minority_i ~ Bernoulli(minority_share)
P(denied=1) = p_bar + beta * (minority_i - minority_share) [centered so mean stays p_bar]
denied_i ~ Bernoulli(that probability), clipped to [0,1]
Fit OLS denied ~ const + minority; test beta_1 with a classic t-test.
Returns the empirical rejection rate (= power if beta != 0, size if beta == 0).
"""
if r is None:
r = np.random.default_rng(0)
if alternative == "two-sided":
z = stats.norm.ppf(1 - alpha / 2)
else:
z = stats.norm.ppf(1 - alpha)
rejects = 0
for _ in range(n_sims):
minority = (r.random(n) < minority_share).astype(float)
prob = p_bar + beta * (minority - minority_share)
prob = np.clip(prob, 1e-6, 1 - 1e-6)
denied = (r.random(n) < prob).astype(float)
# OLS of denied on [1, minority] -> slope and its classical SE.
X = np.column_stack([np.ones(n), minority])
XtX_inv = np.linalg.inv(X.T @ X)
bhat = XtX_inv @ (X.T @ denied)
resid = denied - X @ bhat
s2 = (resid @ resid) / (n - 2)
se_b1 = np.sqrt(s2 * XtX_inv[1, 1])
tstat = bhat[1] / se_b1
if alternative == "greater":
reject = tstat > z
elif alternative == "less":
reject = tstat < -z
else:
reject = abs(tstat) > z
rejects += int(reject)
return rejects / n_sims
# Size check first: beta = 0 must reject ~alpha of the time.
size_hat = simulate_power_regression(0.0, P_BAR, N_REALISTIC, MINORITY_SHARE,
n_sims=3000, alternative="greater", r=rng)
print(f"MC size (beta=0): {size_hat:.4f} (target {ALPHA})")
assert abs(size_hat - ALPHA) < 0.02, "Spec is not calibrated under the null!"
# Power under the alternative beta = 2 pp.
power_mc = simulate_power_regression(delta, P_BAR, N_REALISTIC, MINORITY_SHARE,
n_sims=3000, alternative="greater", r=rng)
# Analytic comparison: for a SIMPLE regression on a binary x with no other controls,
# the residualized regressor SD equals the raw SD, and resid_sd of y ~ sqrt(p(1-p)).
power_an = power_regression_coef(delta, RESID_SD, N_REALISTIC, sd_x_raw,
alpha=ALPHA, alternative="greater")
print(f"MC power (beta={MDE_PP}pp): {power_mc:.4f}")
print(f"Analytic power (same setup) : {power_an:.4f}")
print(f"Difference: {abs(power_mc - power_an):.4f}")
assert abs(power_mc - power_an) < 0.05, "Analytic and simulated power disagree -- investigate!"
print("\nAnalytic ~ simulated: the normal-approx power calc is trustworthy for this design.")5. The power curve vs. N — read off where 80% is crossed
Power is a curve, not a number (the nb1.5 lesson). Sweep the confirmation-set size N at the fixed effect Maya cares about, overlay the analytic curve on the simulated points, and mark the 80% line. Where the curve crosses 80% is exactly the N that section 2 solved for in closed form — a visual confirmation that the inversion and the curve agree.
N_grid = np.array([500, 1000, 2000, 4000, 6000, 8000, 12000, 16000, 24000])
sim_curve = np.array([
simulate_power_regression(delta, P_BAR, int(n), MINORITY_SHARE,
n_sims=1500, alternative="greater", r=rng)
for n in N_grid
])
ana_curve = np.array([
power_regression_coef(delta, RESID_SD, int(n), sd_x_raw, alpha=ALPHA, alternative="greater")
for n in N_grid
])
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(N_grid, sim_curve, "o-", color="darkorange", label="simulated (LPM)")
ax.plot(N_grid, ana_curve, "--", color="navy", label="analytic (normal approx)")
ax.axhline(TARGET_POWER, color="green", linestyle=":", label=f"{int(TARGET_POWER*100)}% power")
ax.axvline(N_needed, color="grey", linestyle=":", label=f"N for 80% = {N_needed}")
ax.set_xlabel("confirmation-set size N")
ax.set_ylabel("power = P(reject H0: beta_1 <= 0)")
ax.set_title(f"Power vs. N for a {MDE_PP} pp denial gap (one-sided, alpha={ALPHA})")
ax.set_ylim(0, 1.02)
ax.legend()
fig.tight_layout()
print(f"The simulated curve crosses 80% near N={N_needed}, matching the closed-form inversion.")6. Minimum detectable effect (MDE) at the N you actually have
Often you cannot choose N — Maya’s confirmation set is whatever the HMDA extract gives her. The right question flips again: given my realistic N, what is the smallest effect I can detect at 80% power? That number is the minimum detectable effect (MDE). Inverting the one-sided power formula for \Delta instead of N:
\text{MDE}=(z_\alpha+z_{1-\beta})\;\sigma\sqrt{\tfrac{1}{n_1}+\tfrac{1}{n_0}} \qquad\text{(two-group)},\qquad \text{MDE}=(z_\alpha+z_{1-\beta})\;\frac{\sigma_\varepsilon}{\sqrt N\,\operatorname{sd}(\tilde x_1)} \qquad\text{(regression coef)}.
If the MDE is larger than the effect Maya says would matter (her 2 pp threshold), her study is underpowered for the question — and the PAP-honest response is one of the options in Reflection Prompt 3 of the chapter, not “run it anyway and hope.”
def mde_two_group(sigma, n1, n0, alpha=ALPHA, power=TARGET_POWER, alternative="greater"):
z_a = stats.norm.ppf(1 - alpha / 2) if alternative == "two-sided" else stats.norm.ppf(1 - alpha)
z_b = stats.norm.ppf(power)
return (z_a + z_b) * sigma * np.sqrt(1.0 / n1 + 1.0 / n0)
def mde_regression(resid_sd, n, sd_x_resid, alpha=ALPHA, power=TARGET_POWER, alternative="greater"):
z_a = stats.norm.ppf(1 - alpha / 2) if alternative == "two-sided" else stats.norm.ppf(1 - alpha)
z_b = stats.norm.ppf(power)
return (z_a + z_b) * resid_sd / (np.sqrt(n) * sd_x_resid)
mde_tg = mde_two_group(RESID_SD, n1, n0, alpha=ALPHA, power=TARGET_POWER, alternative="greater")
mde_reg = mde_regression(resid_sd_y, N_REALISTIC, sd_x_resid, alpha=ALPHA,
power=TARGET_POWER, alternative="greater")
print(f"At N={N_REALISTIC}:")
print(f" MDE (two-group) = {mde_tg*100:.3f} pp")
print(f" MDE (regression coef) = {mde_reg*100:.3f} pp (higher: controls inflate the regressor variance)")
print(f" Effect Maya cares about = {MDE_PP:.3f} pp")
# Cross-check: the MDE fed back into the power calc must return ~80% power.
pw_at_mde = power_two_group(mde_tg, RESID_SD, n1, n0, alpha=ALPHA, alternative="greater")
print(f"\nCheck: power at the two-group MDE = {pw_at_mde:.4f} (target {TARGET_POWER})")
assert abs(pw_at_mde - TARGET_POWER) < 0.01, "MDE inversion is inconsistent with power calc!"
verdict = "ADEQUATELY POWERED" if mde_tg * 100 <= MDE_PP else "UNDERPOWERED"
print(f"\nVerdict (two-group): study is {verdict} to detect a {MDE_PP} pp gap at "
f"{int(TARGET_POWER*100)}% power.")7. Benjamini–Hochberg FDR — input p-values, output discoveries
The multiple-testing plan (Component 3) quarantines Maya’s subgroup analyses into a family and controls the false discovery rate — the expected fraction of declared discoveries that are false — with the Benjamini–Hochberg procedure from §7.3.2. The recipe, exactly as the chapter states it:
Sort the m p-values ascending p_{(1)}\le\dots\le p_{(m)}. Find the largest rank k with p_{(k)}\le\frac{k}{m}\,q. Reject every hypothesis with rank \le k.
We implement it as a function that returns a boolean reject-mask aligned to the original input order (so you can map discoveries straight back to your hypotheses), plus the BH-adjusted q-values. Then we run the worked income-band family — \{0.008,0.012,0.039,0.041,0.330\} — and assert it reproduces the chapter’s answer: BH keeps the two smallest, where Bonferroni keeps only one.
def benjamini_hochberg(pvals, q=0.05):
"""Benjamini-Hochberg FDR procedure.
pvals : 1-D array of p-values (any order).
q : target false discovery rate.
Returns (reject_mask, qvalues), both aligned to the INPUT order.
reject_mask[i] = True if hypothesis i is a declared discovery.
qvalues[i] = BH-adjusted p-value (monotone, capped at 1).
"""
p = np.asarray(pvals, dtype=float)
m = p.size
order = np.argsort(p) # indices that sort p ascending
p_sorted = p[order]
ranks = np.arange(1, m + 1)
thresholds = ranks / m * q # the BH staircase (k/m)*q
below = p_sorted <= thresholds
reject_sorted = np.zeros(m, dtype=bool)
if below.any():
k = np.max(np.where(below)[0]) # largest rank (0-indexed) satisfying the test
reject_sorted[: k + 1] = True # reject ranks 1..k
# BH-adjusted q-values: enforce monotonicity from the largest rank down.
q_sorted = p_sorted * m / ranks
q_sorted = np.minimum.accumulate(q_sorted[::-1])[::-1]
q_sorted = np.clip(q_sorted, 0, 1)
# Map back to original input order.
reject = np.empty(m, dtype=bool); reject[order] = reject_sorted
qvalues = np.empty(m); qvalues[order] = q_sorted
return reject, qvalues
# The chapter's worked income-band family (deliberately scrambled to test order-independence).
labels = ["band_A", "band_B", "band_C", "band_D", "band_E"]
pvals = np.array([0.039, 0.008, 0.330, 0.012, 0.041]) # NOT sorted
reject, qvals = benjamini_hochberg(pvals, q=ALPHA)
print("BH-FDR at q = 0.05:")
for lab, p, r, qv in zip(labels, pvals, reject, qvals):
print(f" {lab}: p={p:.3f} q={qv:.3f} {'DISCOVERY' if r else '--'}")
print(f"\nTotal discoveries: {reject.sum()}")
# The two smallest p-values (0.008 and 0.012) should survive; the rest should not.
expected = (pvals <= 0.012)
assert np.array_equal(reject, expected), "BH did not select the chapter's two bands!"
# Contrast with Bonferroni at the same level: each held to q/m = 0.01 -> only 0.008 survives.
bonf = pvals <= ALPHA / pvals.size
print(f"Bonferroni (q/m = {ALPHA/pvals.size:.3f}) discoveries: {bonf.sum()} "
f"({[labels[i] for i in np.where(bonf)[0]]})")
print("BH keeps the 0.012 band that Bonferroni throws away -- the payoff of FDR over FWER.")A quick calibration sanity check for the BH routine
One more reassurance, in the nb1.5 spirit: if a whole family is null (all p-values uniform on [0,1]), BH at q=0.05 should make false discoveries rarely — its expected false discovery proportion is bounded by q. We simulate many all-null families and confirm the realized FDR sits at or below q.
m_family = 10
n_trials = 4000
false_discovery_props = np.empty(n_trials)
any_reject = 0
for t in range(n_trials):
null_p = rng.random(m_family) # all nulls true -> p ~ Uniform(0,1)
rej, _ = benjamini_hochberg(null_p, q=ALPHA)
n_rej = rej.sum()
false_discovery_props[t] = (n_rej / n_rej) if n_rej > 0 else 0.0 # all rejections are false here
any_reject += int(n_rej > 0)
realized_fdr = false_discovery_props.mean() # E[false discoveries / discoveries]
print(f"All-null families, m={m_family}, {n_trials} trials:")
print(f" Realized FDR (mean false-discovery proportion): {realized_fdr:.4f} (bound = q = {ALPHA})")
print(f" Fraction of families with >=1 false discovery: {any_reject / n_trials:.4f}")
assert realized_fdr <= ALPHA + 0.01, "BH is not controlling FDR at q!"
print("BH holds the false discovery rate at or below q -- the routine is calibrated.")8. The PAP-template emitter
The payoff. This function takes the design choices you have made above — your hypotheses, the seven slots of your primary specification, your multiple-testing families, your falsification conditions — plus the power numbers this notebook just computed — and prints a filled-in short-PAP skeleton in the exact form of §7.3.4. You paste it into your PAP.md, edit the prose, and tag the commit pap-filed. The emitter does not invent your design; it enforces its shape, so no slot (especially the identifying assumption and the falsification line) can be silently left blank — the vagueness peril of §7.3.5 made mechanically hard to commit.
def emit_pap_template(*, title, author, data_note,
hypotheses, spec, mt_families, falsification, robustness,
holdout, power_summary):
"""Print a filled short-form PAP skeleton (§7.3.4 form). All fields required."""
lines = []
A = lines.append
A("=" * 78)
A(f"PRE-ANALYSIS PLAN — {title}")
A(f"Author: {author}")
A(f"Filed (tagged commit `pap-filed`): before any confirmatory regression")
A(f"Data: {data_note}")
A("=" * 78)
A("\n§1. Hypotheses (directional).")
for h in hypotheses:
A(f" - {h}")
A("\n§2. Primary specification (CONVENTIONS §4 form).")
A(f" Equation: {spec['equation']}")
for slot in ["outcome", "treatment", "controls", "fixed_effects",
"clustering", "sample", "identifying_assumption"]:
A(f" - {slot.replace('_', ' ').title()}: {spec[slot]}")
A("\n§3. Multiple-testing plan.")
for fam in mt_families:
A(f" - {fam}")
A("\n§4. Falsification.")
for f in falsification:
A(f" - {f}")
A("\n§5. Planned robustness (reported whichever way they come out).")
for rob in robustness:
A(f" - {rob}")
A("\n§6. Hold-out & registration.")
A(f" - {holdout}")
A("\n§7. Power (computed in nb7.3 on the EXPLORATION set; pasted here as evidence).")
for k, v in power_summary.items():
A(f" - {k}: {v}")
A("=" * 78)
print("\n".join(lines))
# Assemble Maya's PAP from the constants and results computed above.
emit_pap_template(
title="Conditional Racial Disparity in Mortgage Denial, and the FinTech Margin",
author="Maya",
data_note="HMDA 2019-2021 extract, snapshot pinned in data-cards/hmda.md",
hypotheses=[
"H1 (primary). Conditional on observables, lender, and tract, minority applicants face "
"a HIGHER denial rate (beta_1 > 0). One-sided.",
"H2 (secondary). The conditional gap is SMALLER for FinTech lenders (interaction < 0). One-sided.",
"H3 (secondary, exploratory family). The gap differs across {Black, Hispanic, Asian} x "
"{5 income bands}. No directional prior; two-sided.",
],
spec={
"equation": "denied_i = b0 + b1*minority_i + x_i'gamma + alpha_lender + delta_tract + eps_i",
"outcome": "denied = 1 if HMDA action_taken in {3,7}, else 0",
"treatment": "minority (1 if applicant is Black or Hispanic); coefficient b1 answers H1",
"controls": "log loan amount, log applicant income, loan-to-value, loan term. "
"DELIBERATELY EXCLUDED: loan-product type, high-cost flag (on the steering pathway)",
"fixed_effects": "lender (alpha) and census tract (delta)",
"clustering": "by lender (errors plausibly correlated within a lender's underwriting style)",
"sample": "first-lien, owner-occupied, 1-4-family, conventional conforming PURCHASE apps, "
f"2019-2021; drop withdrawn/incomplete/preapproval-only; expected N ~ {N_REALISTIC}",
"identifying_assumption": "Conditional on controls + both FE, minority is uncorrelated with the "
"unobserved creditworthiness HMDA omits; beta_1 is an UPPER BOUND on "
"differential treatment, not a clean causal effect",
},
mt_families=[
"Family A (primary): H1 alone, alpha=0.05, one-sided. No correction (single primary test).",
"Family B: H2 FinTech interaction, alpha=0.05, one-sided.",
"Family C (exploratory): {Black,Hispanic,Asian} x {5 income bands}. Benjamini-Hochberg FDR "
"at q=0.05 WITHIN the family; labeled exploratory regardless of outcome.",
],
falsification=[
f"H1 falsified if the 95% CI for beta_1 contains zero AND excludes economically meaningful "
f"values (gap >= {MDE_PP} pp) -- a tight bound near zero. A WIDE CI including zero is a power "
f"failure, reported as inconclusive.",
"H2 falsified if the interaction is indistinguishable from zero or positive (FinTech worse).",
],
robustness=[
"Re-cluster by county. Drop 2020 (COVID). Add continuous rate-spread outcome. "
"Restrict to 50 largest lenders. Re-estimate on full sample as a precision check (labeled).",
],
holdout="Data split 30% exploration / 70% confirmation on hash(loan_id) mod 10 (pinned, "
"CONVENTIONS §5). Spec designed ONLY on the 30%; primary regression runs ONCE on the "
"untouched 70% after this PAP is tagged. Deviations logged in DEVIATIONS.md.",
power_summary={
"Effect of interest (MDE target)": f"{MDE_PP} pp denial gap",
"Residual SD / baseline rate": f"{RESID_SD:.3f} (p_bar={P_BAR})",
"N for 80% power (one-sided, two-group)": f"{N_needed}",
"Realistic confirmation N": f"{N_REALISTIC}",
"MDE at realistic N (two-group, 80%)": f"{mde_tg*100:.2f} pp",
"Powered for the 2 pp effect?": ("YES" if mde_tg * 100 <= MDE_PP else "NO -- underpowered"),
},
)9. Summary — what you carry into Lab 7
You have built, in one notebook, the quantitative spine of a pre-analysis plan:
- Analytic power for a two-group mean comparison and for a regression coefficient, the latter using the residualized regressor SD so variance inflation does not flatter you;
- the inversion to the N for 80% power, self-checked by feeding the answer back in;
- a Monte Carlo simulator of your actual spec, confirming analytic \approx simulated so you can trust the shortcut;
- the MDE at your realistic N — the number that tells you, before peeking, whether your study can even answer its question;
- a calibrated Benjamini–Hochberg routine that selects the right discoveries from a family;
- a PAP emitter that forces every slot of §7.3.4 to be filled.
The deepest connection back to Week 1: power and size are two readings of the same sampling distribution you simulated in nb1.5 — size is the rejection rate when the effect is zero, power is the rejection rate when it is not — and the multiple-testing arithmetic of §1.5.10 (1-0.95^{20}\approx0.64) is exactly what the BH family control is built to tame. Run the power calc on your exploration set, paste the emitter’s output into your PAP.md, and tag pap-filed. That tag is your promise that you did not fish.
Your Turn
Re-run the power calc for your outcome. Replace
P_BAR,RESID_SD,MINORITY_SHARE, andN_REALISTICwith the numbers from your own exploration set (for a continuous outcome, setRESID_SDto the residual SD of your regression, not \sqrt{p(1-p)}). Report your MDE and whether your design is powered for the smallest effect you care about.Two-sided cost. Recompute
N_for_powerandMDEwithalternative="two-sided"and explain, in one sentence tying back to nb1.5, why a two-sided test demands a larger N for the same power (hint: z_{\alpha/2}=1.96 vs z_\alpha=1.645).Design your family. Build the vector of p-values you expect from your secondary family (you can simulate them: a few real effects buried among nulls), run
benjamini_hochberg, and compare its discoveries to plain Bonferroni. How many real effects does FDR rescue that FWER would have discarded?Emit and tag. Fill
emit_pap_templatewith your design, paste the output intoPAP.md, write yourDEVIATIONS.mdstub, and make thepap-filedtagged commit. That is the deliverable for PS 7.3 and Lab 7.