Experimental Design, Randomization, and Blocking

import numpy as np, pandas as pd, wooldridge as wd, statsmodels.formula.api as smf
d = wd.data("jtrain2")
X = ["age", "educ", "black", "hisp", "married", "nodegree", "re74", "re75"]
g = d.groupby("train")[X]
tab = pd.DataFrame({"control": g.mean().loc[0], "treated": g.mean().loc[1]})
tab["std diff"] = (tab["treated"] - tab["control"]) / np.sqrt((g.var().loc[0] + g.var().loc[1]) / 2)
print(tab.round(3).to_string())
r = smf.ols("train ~ " + " + ".join(X), data=d).fit(cov_type="HC1")
w = r.wald_test(", ".join(f"{x} = 0" for x in X), scalar=True)
print(f"\nomnibus balance test: chi2({len(X)}) = {w.statistic:.2f}, p = {w.pvalue:.3f}")
y, t = d["re78"].to_numpy(), d["train"].to_numpy().astype(bool)
obs = y[t].mean() - y[~t].mean(); rng = np.random.default_rng(11); n1 = t.sum(); hits = 0
for _ in range(5000):
    s = np.zeros(len(y), bool); s[rng.choice(len(y), n1, replace=False)] = True
    hits += abs(y[s].mean() - y[~s].mean()) >= abs(obs)
print(f"randomization-inference p-value for re78: {hits/5000:.4f}  (observed difference {obs:.3f})")

Output:

control  treated  std diff
age        25.054   25.816     0.107
educ       10.088   10.346     0.141
black       0.827    0.843     0.044
hisp        0.108    0.059    -0.175
married     0.154    0.189     0.094
nodegree    0.835    0.708    -0.304
re74        2.107    2.096    -0.002
re75        1.267    1.532     0.084

omnibus balance test: chi2(8) = 17.99, p = 0.021
randomization-inference p-value for re78: 0.0042  (observed difference 1.794)

Most standardized differences are small, but nodegree differs by 0.30 SD and hisp by 0.18, and the omnibus test rejects balance at 5% (p ≈ 0.02). In a genuinely randomized experiment, imbalance on some covariate happens by chance; it is a reason to adjust for these covariates for precision, not evidence that randomization failed. The randomization inference p-value (0.004) is valid by design, with no distributional assumptions.

import numpy as np, wooldridge as wd
c = wd.data("jtrain2").query("train == 0").reset_index(drop=True)
y0 = c["re78"].to_numpy(); y1 = y0 + 1.0; n = len(y0)
block = ((c["re75"] > 0).astype(int) * 2 + c["nodegree"]).to_numpy()
rng = np.random.default_rng(5); cr, bl = [], []
for _ in range(5000):
    t = np.zeros(n, bool); t[rng.choice(n, n // 2, replace=False)] = True
    cr.append(y1[t].mean() - y0[~t].mean())
    est = 0.0
    for b in np.unique(block):
        idx = np.where(block == b)[0]; tb = np.zeros(len(idx), bool); tb[rng.choice(len(idx), len(idx) // 2, replace=False)] = True
        est += len(idx) / n * (y1[idx][tb].mean() - y0[idx][~tb].mean())
    bl.append(est)
for lab, v in (("complete randomization", cr), ("blocked randomization", bl)):
    print(f"{lab:23s}: mean estimate {np.mean(v):.3f} (true 1.000), SD {np.std(v):.3f}")
print(f"variance reduction from blocking: {1 - np.var(bl)/np.var(cr):.1%}")
print("block sizes:", np.bincount(block), " block means of Y(0):", [round(y0[block == b].mean(), 2) for b in range(4)])

Output:

complete randomization : mean estimate 0.993 (true 1.000), SD 0.675
blocked randomization  : mean estimate 0.984 (true 1.000), SD 0.683
variance reduction from blocking: -2.2%
block sizes: [ 25 153  18  64]  block means of Y(0): [np.float64(4.62), np.float64(4.15), np.float64(5.19), np.float64(5.32)]

Both designs are unbiased. Here blocking does not reduce variance — in this simulation it is even slightly worse — because the blocking variables barely predict 1978 earnings (block means range from 4.2 to 5.3, while the within-block SD is about 5.5). Blocking pays off only when blocks are homogeneous in the outcome; blocking on weak predictors costs a little precision through smaller within-block samples. In a real design, choose blocks by their predictive power for the outcome, for example from pilot data.