#!/usr/bin/env python3 """ Building a Dataset Where Every Group Trends Up and the Whole Trends Down ======================================================================= Science Journaling Club, Volume 2 Issue 2, Winter 2026, "Statistics That Fool You". QUESTION -------- Simpson's paradox is usually shown once, with one famous table, and then left alone as a curiosity. We ask three quantitative questions instead. (1) Can a dataset showing the paradox be CONSTRUCTED to order? That is, given a target within-group slope b_w and a target pooled slope b_p of the opposite sign, can we write down a dataset that hits both exactly? (2) How often does the paradox appear BY ACCIDENT? Simulate many grouped datasets in which a group-level lurking variable shifts both x and y, with the lurking variable's effect on y drawn at random each time, and measure the frequency of a sign reversal as a function of confounder strength and group-size imbalance. (3) What happens in the CONTINUOUS case, with no groups at all, only a lurking variable Z? Measure how often the marginal slope of Y on X has the opposite sign to the slope adjusted for Z. MODEL ----- Everything here is a computation. There is no experiment, no survey, no measurement of anything outside this file. The club has no laboratory; the computation IS the experiment, and every number printed below came out of this program. Where real published data appear (Berkeley admissions 1973, the London kidney-stone series) they are transcribed from the cited papers and are marked as such; we recompute the arithmetic, we did not collect the data. The core object is a grouped linear model. Group g (g = 1..G) holds n_g points: x_gi = mu_g + e_gi, e_gi ~ N(0, sx^2) y_gi = nu_g + b_w (x_gi - mu_g) + d_gi, d_gi ~ N(0, sigma^2) so every group has the SAME internal slope b_w, and the groups are displaced from one another in x by mu_g and in y by nu_g. The identity the whole study turns on is exact and elementary. Write T_w = sum_g sum_i (x_gi - xbar_g)^2 (within-group spread of x) T_b = sum_g n_g (xbar_g - xbar)^2 (between-group spread of x) w = T_w / (T_w + T_b) then the ordinary least squares slope fitted to the pooled data is b_pooled = w * b_w + (1 - w) * b_b where b_b is the least-squares slope of the G group means (xbar_g, ybar_g), weighted by n_g. The pooled slope is a WEIGHTED AVERAGE of the within-group slope and the between-group slope, and the weight is decided entirely by how much of the spread in x is inside groups versus between them. That is the whole trick. Simpson's paradox is what a weighted average does when the two things being averaged have opposite signs and the weight is large enough. The constructor inverts this: choose b_w, choose the target b_p, lay out the group means (which fixes T_b) and the within-group spread (which fixes T_w), then solve for the required between-group slope b_b = (b_p - w * b_w) / (1 - w) and place the group means on a line of that slope. No search, no fitting. ASSUMPTIONS ----------- * Linear conditional means, homoscedastic Gaussian noise, independent points. * Every group shares one within-group slope b_w. Real data almost never do. * The lurking variable is group membership and nothing else; there is no measurement error in x, no selection on y, no time ordering. * In the frequency experiment the lurking variable's loading on y, theta, is drawn N(0, tau^2) independently of its loading on x. "By accident" means exactly that: the confounder is as likely to push with the within-group effect as against it. This is a modelling CHOICE and it sets the ceiling of the reversal curve at 1/2. Section 10 shows what other choices would do. * Sign reversal is judged on point estimates unless the output explicitly says "significance-gated", in which case both the within-group slopes and the pooled slope must clear a two-sided t test at alpha = 0.05. LIMITATIONS, i.e. what this model leaves out -------------------------------------------- * Nothing here tells you WHICH answer is correct. The paradox is a statement about arithmetic, not about causation. Deciding whether to report the within-group number or the pooled number requires a causal model that no amount of simulation can supply (Pearl 2014; Hernan, Clayton & Keiding 2011). * Our groups are known and observed. The frightening real case is a lurking variable nobody recorded, and this study cannot estimate how often THAT happens, because that frequency is a property of scientific practice rather than of probability. * The Gaussian, equal-slope, equal-variance world is kind. Heavy tails, unequal within-group slopes and non-linearity all make reversals easier, not harder, so our frequencies should be read as a floor. * The 2x2x2 table experiment uses a uniform Dirichlet prior on cell probabilities, which is a mathematical convention with no claim to describe any real corpus of published tables. SEED ---- Master seed 20260214, fixed at the top of the file and never touched again. Every sub-experiment draws an independent stream from it through numpy's SeedSequence.spawn, so the whole output is deterministic. Run: python simpsons-paradox.py > simpsons-paradox-output.txt Needs: Python 3.12, numpy. Nothing else. """ import math import sys import time from fractions import Fraction import numpy as np MASTER_SEED = 20260214 try: sys.stdout.reconfigure(encoding="utf-8") except Exception: pass T0 = time.time() def rule(title=""): if title: print("\n" + "=" * 78) print(title) print("=" * 78) else: print("-" * 78) def fmt(v): if v == 0: return "0" a = abs(v) if a < 1e-6 or a >= 1e7: return "%.6e" % v return "%.10g" % v def cmp_line(label, club, ref, unit="", tol=None): """Print the club's value beside the accepted or analytic value, with the difference. Every validation in this file goes through here.""" d = club - ref s = " %-44s club %-16s ref %-16s diff %s" % ( label, fmt(club) + unit, fmt(ref) + unit, fmt(d) + unit) if tol is not None: s += " [%s]" % ("OK" if abs(d) <= tol else "MISMATCH") print(s) # =========================================================================== # 0. The club's own regression routine. Everything downstream uses it. # =========================================================================== def ols(x, y): """Simple least squares fit of y on x. Returns (slope, intercept, se_slope, t, df, r). Written from the normal equations, not from a library fitter.""" x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) n = x.size xb = x.mean() yb = y.mean() dx = x - xb dy = y - yb sxx = float(np.dot(dx, dx)) sxy = float(np.dot(dx, dy)) syy = float(np.dot(dy, dy)) slope = sxy / sxx intercept = yb - slope * xb df = n - 2 sse = syy - slope * sxy if df > 0 and sxx > 0: s2 = max(sse, 0.0) / df se = math.sqrt(s2 / sxx) t = slope / se if se > 0 else float("inf") else: se, t = float("nan"), float("nan") r = sxy / math.sqrt(sxx * syy) if sxx > 0 and syy > 0 else float("nan") return slope, intercept, se, t, df, r # --- Student t quantile, needed for the significance-gated reversal rate ---- def _betacf(a, b, x): """Continued fraction for the incomplete beta function (modified Lentz).""" TINY, EPS, ITMAX = 1e-300, 3e-16, 500 qab, qap, qam = a + b, a + 1.0, a - 1.0 c = 1.0 d = 1.0 - qab * x / qap if abs(d) < TINY: d = TINY d = 1.0 / d h = d for m in range(1, ITMAX + 1): m2 = 2 * m aa = m * (b - m) * x / ((qam + m2) * (a + m2)) d = 1.0 + aa * d if abs(d) < TINY: d = TINY c = 1.0 + aa / c if abs(c) < TINY: c = TINY d = 1.0 / d h *= d * c aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)) d = 1.0 + aa * d if abs(d) < TINY: d = TINY c = 1.0 + aa / c if abs(c) < TINY: c = TINY d = 1.0 / d de = d * c h *= de if abs(de - 1.0) < EPS: break return h def betai(a, b, x): """Regularised incomplete beta I_x(a, b).""" if x <= 0.0: return 0.0 if x >= 1.0: return 1.0 lbeta = (math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) + a * math.log(x) + b * math.log1p(-x)) if x < (a + 1.0) / (a + b + 2.0): return math.exp(lbeta) * _betacf(a, b, x) / a return 1.0 - math.exp(lbeta) * _betacf(b, a, 1.0 - x) / b def t_sf(t, df): """Upper tail P(T > t) for Student t with df degrees of freedom.""" if df <= 0: return float("nan") x = df / (df + t * t) p = 0.5 * betai(0.5 * df, 0.5, x) return p if t > 0 else 1.0 - p _TCACHE = {} def t_crit(df, alpha=0.05): """Two-sided critical value by bisection on the survival function.""" key = (int(df), alpha) if key in _TCACHE: return _TCACHE[key] lo, hi = 0.0, 400.0 target = alpha / 2.0 for _ in range(200): mid = 0.5 * (lo + hi) if t_sf(mid, df) > target: lo = mid else: hi = mid v = 0.5 * (lo + hi) _TCACHE[key] = v return v def ncdf(z): return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0))) # =========================================================================== rule("SIMPSON'S PARADOX: CONSTRUCTION, FREQUENCY, AND THE CONTINUOUS CASE") print("Science Journaling Club, Volume 2 Issue 2, Winter 2026") print("Master seed: %d | numpy %s | Python %s" % (MASTER_SEED, np.__version__, sys.version.split()[0])) print("No data were collected. The computation is the experiment.") ss = np.random.SeedSequence(MASTER_SEED) (S_HAND, S_IDENT, S_CONSTRUCT, S_CONTROL, S_GRID, S_CAT, S_CONT, S_TABLE, S_SENS) = ss.spawn(9) # =========================================================================== rule("1. VALIDATION: the regression routine against a closed-form solution") # =========================================================================== print(""" A five-point dataset small enough to check by hand, and small enough that the closed form is exact in rational arithmetic. x = 1,2,3,4,5; y = 2,4,5,4,5. Normal equations: b = Sxy/Sxx, a = ybar - b*xbar. """) hx = [1, 2, 3, 4, 5] hy = [2, 4, 5, 4, 5] fx = [Fraction(v) for v in hx] fy = [Fraction(v) for v in hy] n_h = len(fx) xb_f = sum(fx) / n_h yb_f = sum(fy) / n_h Sxx_f = sum((v - xb_f) ** 2 for v in fx) Sxy_f = sum((a - xb_f) * (b - yb_f) for a, b in zip(fx, fy)) slope_f = Sxy_f / Sxx_f int_f = yb_f - slope_f * xb_f print(" exact rational xbar = %s, ybar = %s" % (xb_f, yb_f)) print(" exact rational Sxx = %s, Sxy = %s" % (Sxx_f, Sxy_f)) print(" exact rational slope = %s = %.17g" % (slope_f, float(slope_f))) print(" exact rational intercept = %s = %.17g" % (int_f, float(int_f))) print() cb, ca, cse, ct, cdf, cr = ols(hx, hy) cmp_line("slope, club ols() vs exact rational", cb, float(slope_f), tol=1e-15) cmp_line("intercept, club ols() vs exact rational", ca, float(int_f), tol=1e-15) np_b, np_a = np.polyfit(np.array(hx, float), np.array(hy, float), 1) cmp_line("slope, club ols() vs numpy.polyfit", cb, float(np_b), tol=1e-12) cmp_line("intercept, club ols() vs numpy.polyfit", ca, float(np_a), tol=1e-12) A = np.vstack([np.array(hx, float), np.ones(5)]).T sol, *_ = np.linalg.lstsq(A, np.array(hy, float), rcond=None) cmp_line("slope, club ols() vs numpy.linalg.lstsq", cb, float(sol[0]), tol=1e-12) syy_f = sum((b - yb_f) ** 2 for b in fy) sse_f = syy_f - slope_f * Sxy_f se_exact = math.sqrt(float(sse_f) / 3.0 / float(Sxx_f)) print() print(" residual sum of squares, exact: Syy - b*Sxy = %s = %.10g" % (sse_f, float(sse_f))) cmp_line("se(slope), club vs closed form", cse, se_exact, tol=1e-14) print(" club routine also reports t = %.6f on df = %d, r = %.10f" % (ct, cdf, cr)) print(""" The t-distribution quantiles, needed later for the significance-gated reversal rate, set against published two-sided 5 percent table values. The club value comes from bisecting our own incomplete-beta implementation.""") for df, tab in [(1, 12.706), (2, 4.303), (5, 2.571), (10, 2.228), (30, 2.042), (60, 2.000), (100, 1.984), (1000, 1.962)]: cmp_line("t_{0.975, df=%d}" % df, t_crit(df), tab, tol=6e-4) print(" Differences are at the rounding of the published tables (3 d.p.).") # =========================================================================== rule("2. VALIDATION: the weighted-average identity, to machine precision") # =========================================================================== print(""" Claim: b_pooled = w * b_within + (1 - w) * b_between, w = T_w/(T_w + T_b), with b_between the n_g-weighted least-squares slope of the group means and b_within the T_w-weighted average of the per-group slopes. Checked on a random grouped dataset with unequal group sizes, unequal within-group spreads and deliberately unequal within-group slopes, none of which the identity requires. """) rng = np.random.default_rng(S_IDENT) sizes_id = [37, 12, 61, 25, 9] xs_id, ys_id, gid = [], [], [] for g, n in enumerate(sizes_id): mu = rng.normal(0, 3) nu = rng.normal(0, 3) bg = rng.normal(0, 2) xg = mu + rng.normal(0, 1.0 + 0.4 * g, n) yg = nu + bg * (xg - xg.mean()) + rng.normal(0, 0.7, n) xs_id.append(xg) ys_id.append(yg) gid.append(np.full(n, g)) X = np.concatenate(xs_id) Y = np.concatenate(ys_id) G_id = np.concatenate(gid) Tw = 0.0 Sw = 0.0 for g in range(len(sizes_id)): m = G_id == g dx = X[m] - X[m].mean() dy = Y[m] - Y[m].mean() Tw += float(np.dot(dx, dx)) Sw += float(np.dot(dx, dy)) gm_x = np.array([X[G_id == g].mean() for g in range(len(sizes_id))]) gm_y = np.array([Y[G_id == g].mean() for g in range(len(sizes_id))]) nn = np.array(sizes_id, float) xbar = float(np.dot(nn, gm_x) / nn.sum()) ybar = float(np.dot(nn, gm_y) / nn.sum()) Tb = float(np.dot(nn, (gm_x - xbar) ** 2)) Sb = float(np.dot(nn, (gm_x - xbar) * (gm_y - ybar))) b_w_eff = Sw / Tw b_b_eff = Sb / Tb w_id = Tw / (Tw + Tb) pred = w_id * b_w_eff + (1 - w_id) * b_b_eff b_pool = ols(X, Y)[0] print(" group sizes %s (N = %d)" % (sizes_id, int(nn.sum()))) print(" T_w = %.10f T_b = %.10f w = %.12f" % (Tw, Tb, w_id)) print(" effective within slope b_w = %.12f" % b_w_eff) print(" between-group slope b_b = %.12f" % b_b_eff) cmp_line("pooled slope, fitted vs identity", b_pool, pred, tol=1e-11) print(" relative difference: %.3e" % (abs(b_pool - pred) / abs(b_pool))) print(" The identity is algebra, so the only thing this can catch is a coding") print(" error. It caught two while we were writing the file.") # =========================================================================== rule("3. THE CONSTRUCTOR") # =========================================================================== print(""" construct(b_within, b_pooled_target, group means, group sizes) places the group means on a line of slope b_b = (b_p - w*b_w)/(1 - w) and draws each group with the requested internal slope. No search is involved; the answer is solved for. The table below is the noiseless construction, where recovery should be at machine precision. Honest i.i.d. noise is added afterwards. """) def construct(b_w, b_p_target, mus, sizes, sx=1.0, noise=0.0, rng=None, y0=0.0): """Build a grouped dataset with within-group slope b_w and pooled OLS slope b_p_target. Returns (x, y, group index, diagnostics).""" mus = np.asarray(mus, float) sizes = np.asarray(sizes, int) G = mus.size N = int(sizes.sum()) xs, gs = [], [] for g in range(G): n = int(sizes[g]) if rng is None: off = np.linspace(-1.0, 1.0, n) * sx * math.sqrt(3.0) else: off = rng.normal(0.0, sx, n) off = off - off.mean() # group mean is exactly mus[g] xs.append(mus[g] + off) gs.append(np.full(n, g)) x = np.concatenate(xs) gidx = np.concatenate(gs) Tw_ = 0.0 for g in range(G): d = x[gidx == g] - mus[g] Tw_ += float(np.dot(d, d)) xbar_ = float(x.mean()) Tb_ = float(np.dot(sizes.astype(float), (mus - xbar_) ** 2)) if Tb_ <= 0: raise ValueError("group means coincide: T_b = 0, no lever to pull") w_ = Tw_ / (Tw_ + Tb_) b_b = (b_p_target - w_ * b_w) / (1.0 - w_) y = np.empty(N) for g in range(G): m = gidx == g nu = y0 + b_b * (mus[g] - xbar_) y[m] = nu + b_w * (x[m] - mus[g]) if noise > 0 and rng is not None: for g in range(G): m = gidx == g e = rng.normal(0.0, noise, int(m.sum())) y[m] = y[m] + (e - e.mean()) # group means stay exact diag = dict(T_w=Tw_, T_b=Tb_, w=w_, b_between=b_b, xbar=xbar_, N=N) return x, y, gidx, diag rng_c = np.random.default_rng(S_CONSTRUCT) TARGETS = [ ("flagship 5 groups", 1.0, -1.0, [-4, -2, 0, 2, 4], [15] * 5), ("mild reversal", 0.5, -0.10, [-3, -1, 1, 3], [40, 40, 40, 40]), ("violent reversal", 0.5, -2.00, [-2, 0, 2], [50, 50, 50]), ("downhill groups, uphill whole", -1.0, 1.0, [-3, 0, 3], [30, 30, 30]), ("nearly flat both ways", 0.1, -0.1, [-6, -2, 2, 6], [25] * 4), ("wildly unequal sizes", 2.0, -0.25, [-5, -1, 3, 7], [8, 60, 12, 120]), ("two groups only", 1.0, -1.0, [-3, 3], [80, 80]), ("ten groups", 0.8, -0.4, list(range(-9, 11, 2)), [12] * 10), ] print(" %-30s %3s %8s %9s %9s %10s %10s" % ("dataset", "G", "b_w tgt", "b_p tgt", "w", "b_between", "max|dev|")) rule() constructed = {} max_dev_all = 0.0 for label, bw, bp, mus, sizes in TARGETS: xx, yy, gg, dd = construct(bw, bp, mus, sizes, sx=1.0) within = [ols(xx[gg == g], yy[gg == g])[0] for g in range(len(sizes))] pooled = ols(xx, yy)[0] dev = max([abs(v - bw) for v in within] + [abs(pooled - bp)]) max_dev_all = max(max_dev_all, dev) print(" %-30s %3d %8.3f %9.3f %9.5f %10.4f %10.2e" % (label, len(sizes), bw, bp, dd["w"], dd["b_between"], dev)) constructed[label] = (xx, yy, gg, dd, within, pooled) rule() cmp_line("worst deviation over all 8 constructions", max_dev_all, 0.0, tol=1e-11) print(" The construction is algebra, so this is a check on the code. The idea") print(" is checked by the fits themselves, printed next.") print(""" The flagship dataset in full. Five groups, within-group slope +1.000 by design, pooled slope -1.000 by design. This is Figure 1.""") x, y, gi, d, within, pooled = constructed["flagship 5 groups"] print(" T_w = %.6f T_b = %.6f w = %.8f" % (d["T_w"], d["T_b"], d["w"])) print(" required between-group slope b_b = %.8f" % d["b_between"]) print(" group | n | xbar | ybar | fitted within slope | target") for g in range(5): m = gi == g b, a, se, t, dfg, r = ols(x[m], y[m]) print(" %5d | %3d | %6.3f | %6.3f | %19.12f | %+0.1f" % (g + 1, int(m.sum()), x[m].mean(), y[m].mean(), b, 1.0)) bp_f, ap_f, sep_f, tp_f, dfp_f, rp_f = ols(x, y) print(" pooled| %3d | %6.3f | %6.3f | %19.12f | %+0.1f" % (x.size, x.mean(), y.mean(), bp_f, -1.0)) cmp_line("flagship within slope (all five identical)", within[0], 1.0, tol=1e-12) cmp_line("flagship pooled slope", pooled, -1.0, tol=1e-12) print(" pooled correlation r = %+.6f; within-group correlation r = %+.6f" % (rp_f, ols(x[gi == 0], y[gi == 0])[5])) print(""" The same construction with honest i.i.d. Gaussian noise, sigma = 0.8, added inside each group. The targets are no longer hit exactly; they are hit within the standard error of the fit, which is all that can be asked of a noisy dataset. z = (fitted - target)/se.""") xn, yn, gin, dn = construct(1.0, -1.0, [-4, -2, 0, 2, 4], [15] * 5, sx=1.0, noise=0.8, rng=rng_c) zs = [] wsl = [] print(" group | fitted | se | z") for g in range(5): m = gin == g b, a, se, t, dfg, r = ols(xn[m], yn[m]) wsl.append(b) z = (b - 1.0) / se zs.append(z) print(" %5d | %9.5f | %8.5f | %+6.2f" % (g + 1, b, se, z)) bpn, apn, sepn, tpn, dfpn, rpn = ols(xn, yn) zp = (bpn - (-1.0)) / sepn zs.append(zp) print(" pooled| %9.5f | %8.5f | %+6.2f" % (bpn, sepn, zp)) print(" largest |z| across the six fits: %.2f" % max(abs(v) for v in zs)) print(" (about 1 is expected; above 3 would mean the constructor is broken)") print(" all five within-group slopes positive: %s" % all(v > 0 for v in wsl)) print(" pooled slope negative: %s" % (bpn < 0)) print(""" How extreme can the target be? b_b = (b_p - w*b_w)/(1 - w) always has a solution, so ANY pair of signs is constructible. The price is the steepness of the between-group line. Holding b_w = +1 and the layout fixed:""") print(" b_p target | required b_between") for bp_t in [-0.1, -0.5, -1.0, -5.0, -20.0]: _, _, _, dz = construct(1.0, bp_t, [-4, -2, 0, 2, 4], [15] * 5, sx=1.0) print(" %10.2f | %19.4f" % (bp_t, dz["b_between"])) print(" w = %.5f for this layout, so the lever arm 1/(1-w) = %.3f" % (dz["w"], 1.0 / (1.0 - dz["w"]))) # =========================================================================== rule("4. THE CONTROL: remove the confounder, balance the groups") # =========================================================================== print(""" The argument needs a control, and here it is a 2x2. Turn the confounder off (c = 0, so all group means coincide in x) and turn imbalance off (equal n_g), then switch each on separately. G = 4 groups, N = 240 points, true within slope b_w = +1, residual sigma = 1. A 'reversal' means every within-group slope has one sign and the pooled slope has the other. """) def group_sizes(N, G, ratio, floor=5): """G sizes summing to N, log-spaced, largest/smallest approx = ratio.""" wts = np.array([float(ratio) ** (k / (G - 1.0)) for k in range(G)]) raw = wts / wts.sum() * N s = np.maximum(np.floor(raw).astype(int), floor) guard = 0 while int(s.sum()) != N and guard < 10000: guard += 1 if int(s.sum()) < N: s[int(np.argmax(raw - s))] += 1 else: cand = np.where(s > floor)[0] j = cand[int(np.argmin((raw - s)[cand]))] s[j] -= 1 return s def run_cell(c, sizes, mode, trials, seed, b_w=1.0, sigma=1.0, tau=1.0, sx=1.0, alpha=0.05, chunk=4000, track=False): """Simulate `trials` grouped datasets and count sign reversals. c : confounder strength, the between-group SD of the group x means, in units of the within-group x SD. mode : 'random' -> group sizes assigned independently of the confounder 'coupled' -> largest group sits at the largest confounder value tau : SD of theta, the confounder's loading on y per unit of its loading on x. theta is drawn fresh for every dataset, so the confounder is as likely to help as to hurt. """ rg = np.random.default_rng(seed) G = len(sizes) N = int(np.sum(sizes)) sizes_sorted = np.sort(np.asarray(sizes)).astype(int) tw_all = np.array([t_crit(int(n) - 2, alpha) for n in sizes_sorted]) tp_all = t_crit(N - 2, alpha) n_rev = 0 n_rev_sig = 0 n_allpos = 0 sum_w = 0.0 done = 0 running = [] while done < trials: m = min(chunk, trials - done) z = rg.normal(0.0, 1.0, (m, G)) theta = rg.normal(0.0, tau, m) order = np.argsort(z, axis=1) if mode == "coupled" else None Sx = np.zeros(m); Sy = np.zeros(m) Sxx = np.zeros(m); Sxy = np.zeros(m); Syy = np.zeros(m) slopes = np.empty((m, G)) tvals = np.empty((m, G)) Tw_acc = np.zeros(m) for k in range(G): n_k = int(sizes_sorted[k]) if mode == "coupled": # slot k holds the k-th smallest group, attached to the k-th # smallest confounder value zk = np.take_along_axis(z, order[:, [k]], axis=1)[:, 0] else: zk = z[:, k] mu = c * zk nu = c * theta * zk xs = mu[:, None] + rg.normal(0.0, sx, (m, n_k)) ys = (nu[:, None] + b_w * (xs - mu[:, None]) + rg.normal(0.0, sigma, (m, n_k))) xm = xs.mean(axis=1, keepdims=True) ym = ys.mean(axis=1, keepdims=True) dx = xs - xm dy = ys - ym sxx = (dx * dx).sum(axis=1) sxy = (dx * dy).sum(axis=1) syy = (dy * dy).sum(axis=1) b = sxy / sxx slopes[:, k] = b dfk = n_k - 2 sse = np.maximum(syy - b * sxy, 1e-300) se = np.sqrt((sse / dfk) / sxx) tvals[:, k] = b / se Tw_acc += sxx Sx += xs.sum(axis=1); Sy += ys.sum(axis=1) Sxx += (xs * xs).sum(axis=1); Sxy += (xs * ys).sum(axis=1) Syy += (ys * ys).sum(axis=1) SXX = Sxx - Sx * Sx / N SXY = Sxy - Sx * Sy / N SYY = Syy - Sy * Sy / N bp = SXY / SXX ssep = np.maximum(SYY - bp * SXY, 1e-300) sep = np.sqrt((ssep / (N - 2)) / SXX) tp = bp / sep allpos = np.all(slopes > 0, axis=1) allneg = np.all(slopes < 0, axis=1) rev = (allpos & (bp < 0)) | (allneg & (bp > 0)) sig_w_pos = np.all(tvals > tw_all[None, :], axis=1) sig_w_neg = np.all(tvals < -tw_all[None, :], axis=1) rev_sig = (sig_w_pos & (tp < -tp_all)) | (sig_w_neg & (tp > tp_all)) n_rev += int(rev.sum()) n_rev_sig += int(rev_sig.sum()) n_allpos += int(allpos.sum()) sum_w += float((Tw_acc / SXX).sum()) if track: running.append(rev.astype(np.int8)) done += m p = n_rev / trials se_p = math.sqrt(max(p * (1 - p), 0.0) / trials) ps = n_rev_sig / trials se_ps = math.sqrt(max(ps * (1 - ps), 0.0) / trials) out = dict(p=p, se=se_p, n_rev=n_rev, p_sig=ps, se_sig=se_ps, n_sig=n_rev_sig, trials=trials, mean_w=sum_w / trials, frac_allpos=n_allpos / trials) if track: out["stream"] = np.concatenate(running) return out N_TOT, G_TOT = 240, 4 TRIALS_CTL = 40000 ctl_seeds = S_CONTROL.spawn(4) bal = group_sizes(N_TOT, G_TOT, 1) imb = group_sizes(N_TOT, G_TOT, 30) print(" balanced sizes : %s (sum %d)" % (list(bal), int(bal.sum()))) print(" imbalanced sizes: %s (sum %d)" % (list(imb), int(imb.sum()))) print() print(" %-44s %10s %12s %10s" % ("cell", "reversals", "rate", "SE")) rule() ctl_rows = [] for i, (lab, c, sz) in enumerate([ ("no confounder (c=0), balanced [CONTROL]", 0.0, bal), ("no confounder (c=0), imbalance 30:1", 0.0, imb), ("confounder c=1, balanced", 1.0, bal), ("confounder c=1, imbalance 30:1", 1.0, imb)]): r = run_cell(c, sz, "random", TRIALS_CTL, ctl_seeds[i]) ctl_rows.append((lab, r)) print(" %-44s %10d %12.6f %10.6f" % (lab, r["n_rev"], r["p"], r["se"])) rule() p0 = ctl_rows[0][1]["p"] n0 = ctl_rows[0][1]["n_rev"] print(" The control cell produced %d reversals in %d datasets." % (n0, TRIALS_CTL)) cmp_line("control reversal rate vs analytic 0", p0, 0.0, tol=1e-12) if n0 == 0: print(" Zero events in %d trials. By the rule of three the one-sided 95%%" % TRIALS_CTL) print(" upper bound on the true rate is 3/%d = %.2e." % (TRIALS_CTL, 3.0 / TRIALS_CTL)) print(""" Why the control must give exactly zero and not merely a small number: with c = 0 every group mean in x is the same, so T_b = 0, so w = 1 exactly, so b_pooled = b_within as an algebraic identity rather than an approximation. Pooling cannot reverse anything unless it is given a between-group displacement to work with. Imbalance on its own (row 2) also gives zero, for the same reason. The confounder is doing all the work.""") # =========================================================================== rule("5. HOW OFTEN BY ACCIDENT: the confounder / imbalance grid") # =========================================================================== print(""" G = 4 groups, N = 240 points, true within-group slope b_w = +1, residual sigma = 1, within-group x spread sx = 1. The group-level lurking variable z_g is standard normal. It shifts the group's x mean by c*z_g and the group's y level by c*theta*z_g, with theta ~ N(0,1) drawn fresh for each dataset, so the confounder pushes against the within-group effect exactly half the time. 'random' : group sizes assigned independently of the confounder. 'coupled' : the largest group sits at the largest value of the confounder, which is the structure of the Berkeley admissions data, where the biggest departments were also the most selective. rate : fraction of datasets where all within-group slopes share one sign and the pooled slope takes the other. gated : the same event, but every within-group slope and the pooled slope must also clear a two-sided t test at alpha = 0.05. mean w : mean of T_w/(T_w+T_b), the weight the pooled fit gives within-group information. """) C_GRID = [0.0, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0] R_GRID = [1, 3, 10, 30] TRIALS_GRID = 20000 grid_seeds = S_GRID.spawn(len(C_GRID) * len(R_GRID) * 2) grid = {} k = 0 print(" %-8s %6s %8s %10s %10s %10s %10s %8s" % ("mode", "c", "ratio", "rate", "SE", "gated", "SE", "mean w")) rule() for mode in ("random", "coupled"): for r in R_GRID: sz = group_sizes(N_TOT, G_TOT, r) for c in C_GRID: res = run_cell(c, sz, mode, TRIALS_GRID, grid_seeds[k]) k += 1 grid[(mode, r, c)] = res print(" %-8s %6.2f %8s %10.5f %10.5f %10.5f %10.5f %8.4f" % (mode, c, "%d:1" % r, res["p"], res["se"], res["p_sig"], res["se_sig"], res["mean_w"])) print() rule() print(""" Read the table twice. The first reading is the obvious one: reversal frequency climbs with confounder strength, from exactly zero at c = 0 to something close to one half once the between-group spread of x is several times the within-group spread. The ceiling of 1/2 is not a discovery; it is the modelling choice that theta is symmetric about zero. WHERE the curve reaches that ceiling is the finding.""") print() for c in C_GRID: if c > 0: res = grid[("random", 1, c)] print(" c = %-5.2f rate = %.4f +/- %.4f gated = %.4f +/- %.4f" % (c, res["p"], res["se"], res["p_sig"], res["se_sig"])) print(""" The second reading surprised us. Compare the 'random' rows across imbalance ratios at fixed confounder strength:""") for c in [0.5, 1.0, 2.0]: row = " c = %.2f :" % c for r in R_GRID: row += " %2d:1 -> %.4f" % (r, grid[("random", r, c)]["p"]) print(row) print(""" Imbalance in GROUP SIZE, on its own, makes the regression paradox LESS likely, not more. The reason is in the arithmetic: for two groups of sizes n1 and n2 separated in x by D, the between-group sum of squares is T_b = (n1*n2/N)*D^2, which is largest when the groups are the same size. A dataset dominated by one huge group has very little between-group leverage, so the pooled line is pinned to that group's own slope. Watch the mean-w column climb as the ratio grows: at c = 1 it goes 0.62, 0.63, 0.68, 0.73. We expected coupling size to the confounder to undo that, on the Berkeley analogy where the biggest departments were also the most selective. It does not. It makes reversals rarer still:""") for c in [0.5, 1.0, 2.0]: row = " c = %.2f :" % c for r in R_GRID: row += " %2d:1 -> %.4f" % (r, grid[("coupled", r, c)]["p"]) print(row) print(" difference (coupled minus random), with the SE of the difference:") for c in [0.5, 1.0, 2.0]: for r in [10, 30]: a = grid[("coupled", r, c)] b = grid[("random", r, c)] dd2 = a["p"] - b["p"] sd = math.sqrt(a["se"] ** 2 + b["se"] ** 2) print(" c=%.2f, %2d:1 : %+0.5f +/- %.5f (%+.2f SE)" % (c, r, dd2, sd, dd2 / sd if sd > 0 else float("nan"))) print(""" Every difference is negative and several are more than ten standard errors from zero, so this is not noise. The mechanism is the grand mean. Put the 165-point group at the extreme of the confounder and xbar follows it there, which shrinks that group's own deviation (xbar_g - xbar) to almost nothing and leaves the between-group sum of squares smaller than before. So in the regression form of the paradox, unequal group sizes are a protection rather than a hazard. That is the opposite of what we were told, and it is the single result in this study we would most like someone to check. Section 5b shows where the folklore is actually right.""") hl = grid[("random", 1, 1.0)] print(""" HEADLINE. With a confounder as strong as the within-group spread (c = 1), four balanced groups and 240 points, %.2f%% of datasets showed a full sign reversal, and %.2f%% survived a significance gate on every slope involved. The second number is the one to be frightened of. It is the rate at which this model produces a dataset in which a careful analyst, testing each group and then the pooled data, would find significant effects pointing in opposite directions.""" % (100 * hl["p"], 100 * hl["p_sig"])) print(" plain: %d reversals / %d trials, SE %.5f" % (hl["n_rev"], hl["trials"], hl["se"])) print(" gated: %d reversals / %d trials, SE %.5f" % (hl["n_sig"], hl["trials"], hl["se_sig"])) # =========================================================================== rule("5b. THE CATEGORICAL FORM: where imbalance really is the culprit") # =========================================================================== print(""" The regression paradox and the rate-table paradox are not the same arithmetic, and they do not have the same risk factors. In a rate table there is no leverage weighting. The pooled rate for an arm is just its stratum rates averaged with the arm's OWN allocation weights, so a reversal needs the two arms to be allocated DIFFERENTLY across the strata. That, and not group size, is the imbalance the folklore means. Two strata, two arms, n = 350 per arm, matching the kidney-stone series. q_k ~ U(0.05, 0.85) stratum-k success probability for arm B delta ~ U(0, 0.10) arm A's advantage, applied in BOTH strata f_A = 0.5 + s*h/2, f_B = 0.5 - s*h/2, s = +/-1 at random h is the allocation gap: h = 0 means both arms send half their patients to each stratum, h = 0.9 means one arm sends 95 percent of its patients to a stratum the other arm barely touches. The population reversal condition is exactly s*h*(q1 - q2) + delta < 0, which integrates to a closed form: P(h) = 5h * (L/3) * [1 - ((L - U)/L)^3], L = 0.8, U = min(0.1/h, L) The upper end of q is 0.85 rather than 0.95 so that q_k + delta never needs clipping at 1. Our first version used 0.95, the clipping bit on about 3 percent of draws, and the Monte Carlo sat 3.4 standard errors off the closed form. The closed form was right and the simulation was wrong. With the clipping gone the worst cell sits inside 2.1 standard errors. This experiment, like the continuous one, has an analytic answer to be checked against. """) L_UNIF = 0.8 DELTA_MAX = 0.10 def analytic_cat(h, L=L_UNIF, dmax=DELTA_MAX): if h <= 0: return 0.0 U = min(dmax / h, L) return 0.5 * (h / dmax) * (L / 3.0) * (1.0 - ((L - U) / L) ** 3) H_GRID = [0.0, 0.1, 0.2, 0.3, 0.5, 0.7, 0.9] TRIALS_CAT = 200000 N_ARM = 350 cat_seeds = S_CAT.spawn(len(H_GRID)) cat = {} print(" %6s %12s %10s %12s %11s %8s %11s %10s" % ("h", "MC pop", "SE", "closed form", "diff", "z", "observable", "SE")) rule() for j, h in enumerate(H_GRID): rgc = np.random.default_rng(cat_seeds[j]) q1 = rgc.uniform(0.05, 0.85, TRIALS_CAT) q2 = rgc.uniform(0.05, 0.85, TRIALS_CAT) delta = rgc.uniform(0.0, DELTA_MAX, TRIALS_CAT) sgn = rgc.choice(np.array([-1.0, 1.0]), TRIALS_CAT) fA = 0.5 + sgn * h / 2.0 fB = 0.5 - sgn * h / 2.0 pA1 = np.clip(q1 + delta, 0.0, 1.0) pA2 = np.clip(q2 + delta, 0.0, 1.0) poolA = fA * pA1 + (1 - fA) * pA2 poolB = fB * q1 + (1 - fB) * q2 rev = poolA < poolB p_pop = float(rev.mean()) se_pop = math.sqrt(max(p_pop * (1 - p_pop), 0.0) / TRIALS_CAT) aval = analytic_cat(h) zc = (p_pop - aval) / se_pop if se_pop > 0 else float("nan") # what an analyst would see: integer counts, judged on observed rates nA1 = np.rint(fA * N_ARM).astype(int) nB1 = np.rint(fB * N_ARM).astype(int) nA2 = N_ARM - nA1 nB2 = N_ARM - nB1 sA1 = rgc.binomial(nA1, pA1); sA2 = rgc.binomial(nA2, pA2) sB1 = rgc.binomial(nB1, q1); sB2 = rgc.binomial(nB2, q2) with np.errstate(invalid="ignore", divide="ignore"): oA1 = np.where(nA1 > 0, sA1 / np.maximum(nA1, 1), np.nan) oA2 = np.where(nA2 > 0, sA2 / np.maximum(nA2, 1), np.nan) oB1 = np.where(nB1 > 0, sB1 / np.maximum(nB1, 1), np.nan) oB2 = np.where(nB2 > 0, sB2 / np.maximum(nB2, 1), np.nan) okA = (oA1 > oB1) & (oA2 > oB2) okB = (oA1 < oB1) & (oA2 < oB2) poolobsA = (sA1 + sA2) / float(N_ARM) poolobsB = (sB1 + sB2) / float(N_ARM) obs = ((okA & (poolobsA < poolobsB)) | (okB & (poolobsA > poolobsB))) obs = np.where(np.isnan(oA1) | np.isnan(oA2) | np.isnan(oB1) | np.isnan(oB2), False, obs) p_obs = float(obs.mean()) se_obs = math.sqrt(max(p_obs * (1 - p_obs), 0.0) / TRIALS_CAT) cat[h] = dict(p_pop=p_pop, se_pop=se_pop, analytic=aval, z=zc, p_obs=p_obs, se_obs=se_obs) print(" %6.2f %12.6f %10.6f %12.6f %+11.6f %8.2f %11.6f %10.6f" % (h, p_pop, se_pop, aval, p_pop - aval, zc, p_obs, se_obs)) rule() zc_list = [abs(cat[h]["z"]) for h in H_GRID if cat[h]["se_pop"] > 0] cmp_line("largest |z|, categorical MC vs closed form", max(zc_list), 0.0, tol=3.0) print(" h = 0 gives exactly %d reversals in %d trials, which is the" % (int(round(cat[0.0]["p_pop"] * TRIALS_CAT)), TRIALS_CAT)) print(" categorical control: equal allocation, no paradox, ever.") print(""" This is where the folklore is right. Nothing about the SIZE of the strata matters here. What matters is that the two arms are allocated differently across them, which is exactly what happened to the kidney-stone series in section 9: the surgeons sent 75 percent of their open-surgery cases to the hard stratum and only 23 percent of their keyhole cases, a gap of h = 0.52. Reading that gap off our curve gives a reversal probability of about %.2f for a treatment whose real advantage is somewhere in 0 to 10 points.""" % analytic_cat(0.523)) # =========================================================================== rule("6. THE CONTINUOUS CASE: one lurking variable, no groups at all") # =========================================================================== print(""" Drop the groups entirely. Structural model: Z ~ N(0,1) X = a*Z + u, u ~ N(0,1) Y = beta*X + b*Z + v, v ~ N(0,1) beta = +0.5 is the effect of X on Y holding Z fixed, positive by construction. The MARGINAL slope of Y on X ignores Z: b_marg = Cov(X,Y)/Var(X) = beta + a*b/(a^2 + 1) so a reversal happens exactly when a*b/(a^2+1) < -beta. The loadings a and b are drawn N(0, s^2) with s the 'lurking strength', independently, so once again the lurking variable is as likely to push with the effect as against it. That reversal probability is a one-dimensional integral, evaluable to high precision, which gives an analytic value to check the Monte Carlo against: P = 2 * Int_0^inf phi_s(a) [1 - Phi(beta(a^2+1)/(a*s))] da. """) BETA_TRUE = 0.5 def analytic_reversal(s, beta=BETA_TRUE, A_mult=14.0, nseg=200000): """Composite Simpson on a in (0, A].""" A = A_mult * s h = A / nseg def f(a): if a <= 0: return 0.0 phi = math.exp(-0.5 * (a / s) ** 2) / (s * math.sqrt(2 * math.pi)) return phi * (1.0 - ncdf(beta * (a * a + 1.0) / (a * s))) tot = f(1e-12) + f(A) for i in range(1, nseg): tot += (4.0 if i % 2 else 2.0) * f(i * h) return 2.0 * tot * h / 3.0 S_LIST = [0.25, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0, 8.0] print(" Analytic reversal probability, with a coarser Simpson grid as a") print(" convergence check on the quadrature itself:") print(" %8s %14s %14s %12s" % ("s", "P (200k seg)", "P (50k seg)", "diff")) ana = {} for s in S_LIST: p_fine = analytic_reversal(s) p_coarse = analytic_reversal(s, nseg=50000) ana[s] = p_fine print(" %8.2f %14.9f %14.9f %12.2e" % (s, p_fine, p_coarse, p_fine - p_coarse)) print(""" Now the Monte Carlo, run two ways. 'population' draws a and b and evaluates the exact population marginal slope, so it must converge to the integral above: a direct check on the simulation machinery. 'sample' draws n = 200 observations per dataset and judges the reversal from the FITTED marginal and Z-adjusted slopes, which is what an analyst would actually see. The adjusted slope comes from residualising both X and Y on Z, the Frisch-Waugh route through the club's own least squares. """) N_OBS = 200 TRIALS_POP = 400000 TRIALS_SAMP = 30000 cont_seeds = S_CONT.spawn(len(S_LIST) * 2 + 1) cont = {} print(" %6s %12s %10s %12s %11s %8s %10s %10s" % ("s", "MC pop", "SE", "analytic", "diff", "z", "MC samp", "SE")) rule() ci = 0 for s in S_LIST: rg1 = np.random.default_rng(cont_seeds[ci]); ci += 1 a1 = rg1.normal(0.0, s, TRIALS_POP) b1 = rg1.normal(0.0, s, TRIALS_POP) bm = BETA_TRUE + a1 * b1 / (a1 * a1 + 1.0) p_mc = float((bm < 0).mean()) se_mc = math.sqrt(p_mc * (1 - p_mc) / TRIALS_POP) if se_mc > 0: zsc = (p_mc - ana[s]) / se_mc else: # zero events: compare against the Poisson expectation instead, so the # cell is still checked rather than quietly skipped zsc = float("nan") rg2 = np.random.default_rng(cont_seeds[ci]); ci += 1 a2 = rg2.normal(0.0, s, TRIALS_SAMP) b2 = rg2.normal(0.0, s, TRIALS_SAMP) Zs = rg2.normal(0.0, 1.0, (TRIALS_SAMP, N_OBS)) Xs = a2[:, None] * Zs + rg2.normal(0.0, 1.0, (TRIALS_SAMP, N_OBS)) Ys = (BETA_TRUE * Xs + b2[:, None] * Zs + rg2.normal(0.0, 1.0, (TRIALS_SAMP, N_OBS))) Zc = Zs - Zs.mean(1, keepdims=True) Xc = Xs - Xs.mean(1, keepdims=True) Yc = Ys - Ys.mean(1, keepdims=True) zz = (Zc * Zc).sum(1) rx = Xc - ((Xc * Zc).sum(1) / zz)[:, None] * Zc ry = Yc - ((Yc * Zc).sum(1) / zz)[:, None] * Zc b_adj = (rx * ry).sum(1) / (rx * rx).sum(1) b_marg = (Xc * Yc).sum(1) / (Xc * Xc).sum(1) p_s = float((np.sign(b_marg) != np.sign(b_adj)).mean()) se_s = math.sqrt(p_s * (1 - p_s) / TRIALS_SAMP) cont[s] = dict(p_pop=p_mc, se_pop=se_mc, p_samp=p_s, se_samp=se_s, analytic=ana[s], z=zsc) print(" %6.2f %12.6f %10.6f %12.6f %+11.6f %8.2f %10.6f %10.6f" % (s, p_mc, se_mc, ana[s], p_mc - ana[s], zsc, p_s, se_s)) rule() zlist = [abs(cont[s]["z"]) for s in S_LIST if cont[s]["se_pop"] > 0] zmax = max(zlist) cmp_line("largest |z|, Monte Carlo vs analytic integral", zmax, 0.0, tol=3.0) if zmax < 3.0: print(" Every cell with a non-zero count agrees with the closed form") print(" inside 3 standard errors.") else: print(" WARNING: a cell disagrees by more than 3 SE. Investigate.") s0 = S_LIST[0] exp0 = ana[s0] * TRIALS_POP print(" The s = %.2f cell produced zero events, so it has no usable normal" % s0) print(" standard error. The integral predicts %.4f events in %d draws, and" % (exp0, TRIALS_POP)) print(" the Poisson probability of seeing none of them is %.4f, so zero is" % math.exp(-exp0)) print(" the expected outcome rather than a failure.") print(""" The two columns are not the same thing and the difference has a shape. Where the population rate is SMALL, the sample column sits above it: with n = 200 the fitted marginal slope scatters around its population value, so datasets whose true marginal slope is barely positive come out negative often enough to matter, and those extra reversals are pure sampling noise laid on top of the confounding. Where the population rate is already large the two agree inside a standard error, because a reversal that is going to happen at all happens whether or not the sample is noisy. Sampling error manufactures reversals only where confounding on its own would not have produced any.""") for s in [0.5, 1.0, 2.0, 5.0]: dq = cont[s] dd3 = dq["p_samp"] - dq["p_pop"] sd = math.sqrt(dq["se_samp"] ** 2 + dq["se_pop"] ** 2) print(" s = %.2f : sample %.5f, population %.5f, excess %+0.5f (%.1f SE)" % (s, dq["p_samp"], dq["p_pop"], dd3, dd3 / sd)) # =========================================================================== rule("7. A PUBLISHED NUMBER TO CHECK AGAINST: random 2x2x2 tables") # =========================================================================== print(""" Pavlides and Perlman (2009) give, via a proof due to Hadjicostas, the probability that a random 2x2x2 table exhibits Simpson's paradox when the eight cell probabilities are drawn uniformly from the simplex: exactly 1/60. That is a genuine published number with an exact value, so it is the best check available on our machinery for the categorical form of the paradox. We draw Dirichlet(1,...,1) tables, compute the two stratum-specific success rates and the pooled success rates, and count. """) TRIALS_TAB = 3000000 rgt = np.random.default_rng(S_TABLE) CHUNK = 500000 cnt_dir = 0 cnt_any = 0 done = 0 track_n, track_p, track_se = [], [], [] run_dir = 0 checkpoints = sorted(set([int(round(10 ** (e / 8.0))) for e in range(24, 8 * 7 + 1)] + [TRIALS_TAB])) checkpoints = [cq for cq in checkpoints if 1000 <= cq <= TRIALS_TAB] cp_i = 0 while done < TRIALS_TAB: m = min(CHUNK, TRIALS_TAB - done) ex = rgt.exponential(1.0, (m, 8)) p = ex / ex.sum(axis=1, keepdims=True) # columns: A-success-1, A-fail-1, A-success-2, A-fail-2, then the same for B aS1, aF1, aS2, aF2, bS1, bF1, bS2, bF2 = (p[:, i] for i in range(8)) rA1 = aS1 / (aS1 + aF1); rB1 = bS1 / (bS1 + bF1) rA2 = aS2 / (aS2 + aF2); rB2 = bS2 / (bS2 + bF2) RA = (aS1 + aS2) / (aS1 + aF1 + aS2 + aF2) RB = (bS1 + bS2) / (bS1 + bF1 + bS2 + bF2) d1 = (rA1 > rB1) & (rA2 > rB2) & (RA < RB) d2 = (rA1 < rB1) & (rA2 < rB2) & (RA > RB) cnt_dir += int(d1.sum()) cnt_any += int((d1 | d2).sum()) cs = np.cumsum(d1) while cp_i < len(checkpoints) and checkpoints[cp_i] <= done + m: idx = checkpoints[cp_i] - done - 1 tot = run_dir + int(cs[idx]) nn_ = checkpoints[cp_i] ph = tot / nn_ track_n.append(nn_) track_p.append(ph) track_se.append(math.sqrt(max(ph * (1 - ph), 0.0) / nn_)) cp_i += 1 run_dir += int(cs[-1]) done += m p_dir = cnt_dir / TRIALS_TAB se_dir = math.sqrt(p_dir * (1 - p_dir) / TRIALS_TAB) p_any = cnt_any / TRIALS_TAB se_any = math.sqrt(p_any * (1 - p_any) / TRIALS_TAB) PP = 1.0 / 60.0 print(" trials: %d" % TRIALS_TAB) print(" one-directional paradox (A better in both strata, worse pooled):") print(" count %d, rate %.7f, SE %.7f" % (cnt_dir, p_dir, se_dir)) print(" either direction:") print(" count %d, rate %.7f, SE %.7f" % (cnt_any, p_any, se_any)) print() cmp_line("directional rate vs published 1/60", p_dir, PP, tol=4 * se_dir) print(" %+.2f standard errors from 1/60" % ((p_dir - PP) / se_dir)) cmp_line("either-direction rate vs published 1/60", p_any, PP, tol=4 * se_any) print(" %+.2f standard errors from 1/60" % ((p_any - PP) / se_any)) cmp_line("either-direction rate vs 2/60", p_any, 2 * PP, tol=4 * se_any) print(" %+.2f standard errors from 2/60" % ((p_any - 2 * PP) / se_any)) cmp_line("directional rate vs 1/120", p_dir, PP / 2.0, tol=4 * se_dir) print(" %+.2f standard errors from 1/120" % ((p_dir - PP / 2) / se_dir)) print(" Two of those four lines are MEANT to fail. They are printed so that") print(" the reader can see which convention is ruled out rather than being") print(" shown only the one that agrees.") print(""" Whichever of the two conventions the published 1/60 refers to, one of our two counts lands on it and the other lands on exactly twice it, as the symmetry between the two treatment labels requires. We print both and say which is which rather than quietly choosing the one that agrees.""") # =========================================================================== rule("8. CONVERGENCE") # =========================================================================== print(""" The running estimate of the directional 2x2x2 rate as trials accumulate, with its own standard error, against its exact value. This is Figure 4. It is the only quantity in the study with an exact known answer AND a Monte Carlo estimate, so it is the one place convergence can be watched against truth rather than against itself. """) print(" The quantity tracked is the ONE-DIRECTIONAL rate, whose exact value") print(" is 1/120 = %.7f, half the published 1/60 because 1/60 counts both" % (PP / 2)) print(" directions. The z column is against 1/120.") print(" %10s %12s %12s %12s %8s" % ("trials", "estimate", "SE", "1/120", "z")) show = set([1000, 3162, 10000, 31623, 100000, 316228, 1000000, 3000000]) for nn_, ph, sh in zip(track_n, track_p, track_se): if nn_ in show: zq = (ph - PP / 2) / sh if sh > 0 else float("nan") print(" %10d %12.7f %12.7f %12.7f %+8.2f" % (nn_, ph, sh, PP / 2, zq)) print("\n FIGDATA convergence (trials, estimate, se):") for nn_, ph, sh in zip(track_n, track_p, track_se): print(" CONV %d %.8f %.8f" % (nn_, ph, sh)) # =========================================================================== rule("9. TWO REAL PUBLISHED TABLES, RECOMPUTED") # =========================================================================== print(""" The club collected no data. The two tables below are transcribed from the cited papers; what is ours is the arithmetic performed on them, redone here so that the constructed and simulated results above can be set beside something that actually happened to real people. """) print(" (a) Berkeley graduate admissions, autumn 1973, six largest") print(" departments. Transcribed from Bickel, Hammel & O'Connell (1975),") print(" Science 187, 398-404.") berk = [("A", 825, 62, 108, 82), ("B", 560, 63, 25, 68), ("C", 325, 37, 593, 34), ("D", 417, 33, 375, 35), ("E", 191, 28, 393, 24), ("F", 373, 6, 341, 7)] print(" dept men apps men % women apps women % women-men (pp)") mm = wm = 0.0 ma = wa = 0 for dp, na, pa, nb, pb in berk: mm += na * pa / 100.0 wm += nb * pb / 100.0 ma += na wa += nb print(" %-4s %9d %6d %11d %8d %+15d" % (dp, na, pa, nb, pb, pb - pa)) rm = mm / ma rw = wm / wa print(" pooled: men %d/%d = %.4f (%.1f%%), women %d/%d = %.4f (%.1f%%)" % (round(mm), ma, rm, 100 * rm, round(wm), wa, rw, 100 * rw)) print(" pooled difference (women minus men): %+.1f percentage points" % (100 * (rw - rm))) nd = sum(1 for dp, na, pa, nb, pb in berk if pb > pa) print(" departments admitting women at the higher rate: %d of 6" % nd) tot = [(na + nb) for _, na, _, nb, _ in berk] Tt = sum(tot) std_m = sum(t * e[2] / 100.0 for t, e in zip(tot, berk)) / Tt std_w = sum(t * e[4] / 100.0 for t, e in zip(tot, berk)) / Tt print(" directly standardised to the combined applicant pool:") print(" men %.4f (%.1f%%), women %.4f (%.1f%%), difference %+.1f pp" % (std_m, 100 * std_m, std_w, 100 * std_w, 100 * (std_w - std_m))) print(" SIGN FLIP: pooled %+.1f pp, standardised %+.1f pp." % (100 * (rw - rm), 100 * (std_w - std_m))) print(" The campus-wide figures quoted in the same paper are 8442 male") print(" applicants at 44%, 4321 female applicants at 35%.") print(""" (b) Renal calculi, open surgery against percutaneous nephrolithotomy. Transcribed from Charig, Webb, Payne & Wickham (1986), BMJ 292, 879-882, in the arrangement popularised by Julious & Mullee (1994).""") stones = [("small (<2 cm)", 81, 87, 234, 270), ("large (>=2 cm)", 192, 263, 55, 80)] oa = ob = ta = tb = 0 print(" stratum open surgery PCNL difference") for lab, sa, na, sb, nb in stones: oa += sa; ta += na; ob += sb; tb += nb print(" %-15s %5d/%-4d %5.1f%% %5d/%-4d %5.1f%% %+7.1f pp" % (lab, sa, na, 100 * sa / na, sb, nb, 100 * sb / nb, 100 * (sa / na - sb / nb))) print(" %-15s %5d/%-4d %5.1f%% %5d/%-4d %5.1f%% %+7.1f pp" % ("pooled", oa, ta, 100 * oa / ta, ob, tb, 100 * ob / tb, 100 * (oa / ta - ob / tb))) print(" SIGN FLIP: open surgery wins both strata by +%.1f and +%.1f pp" % (100 * (81 / 87 - 234 / 270), 100 * (192 / 263 - 55 / 80))) print(" and loses the pooled comparison by %.1f pp." % (100 * (ob / tb - oa / ta))) print(" The confounder is stone size: %.1f%% of the open-surgery cases" % (100 * 263 / 350)) print(" were large stones against %.1f%% of the PCNL cases, a gap of" % (100 * 80 / 350)) print(" %.1f percentage points in exposure to the harder problem." % (100 * (263 / 350 - 80 / 350))) # =========================================================================== rule("10. WHERE A DIFFERENT MODELLING CHOICE WOULD CHANGE THE ANSWER") # =========================================================================== print(""" Six knobs we set by hand, each re-run at c = 1 with four balanced groups, N = 240 and %d trials, changing one thing at a time. """ % TRIALS_GRID) sens_seeds = S_SENS.spawn(32) si = 0 print(" %-42s %10s %10s %18s" % ("variant", "rate", "SE", "vs base")) rule() base_r = grid[("random", 1, 1.0)] print(" %-42s %10.5f %10.5f %18s" % ("BASE: G=4, N=240, sigma=1, tau=1, b_w=1", base_r["p"], base_r["se"], "-")) variants = [ ("theta SD tau = 0.5 (weak y-loading)", dict(tau=0.5)), ("theta SD tau = 2.0 (strong y-loading)", dict(tau=2.0)), ("residual sigma = 0.25 (tight groups)", dict(sigma=0.25)), ("residual sigma = 4.0 (noisy groups)", dict(sigma=4.0)), ("within slope b_w = 0.25 (weak signal)", dict(b_w=0.25)), ("within slope b_w = 4.0 (strong signal)", dict(b_w=4.0)), ] sens = {} for lab, kw in variants: r = run_cell(1.0, bal, "random", TRIALS_GRID, sens_seeds[si], **kw) si += 1 sens[lab] = r dd4 = r["p"] - base_r["p"] sd = math.sqrt(r["se"] ** 2 + base_r["se"] ** 2) print(" %-42s %10.5f %10.5f %+9.5f (%+.1f SE)" % (lab, r["p"], r["se"], dd4, dd4 / sd if sd else float("nan"))) print("\n Number of groups, holding N = 240 and c = 1:") gtab = {} for G in (2, 3, 4, 6, 8, 12): sz = group_sizes(240, G, 1) r = run_cell(1.0, sz, "random", TRIALS_GRID, sens_seeds[si]) si += 1 gtab[G] = r print(" G = %2d (n_g = %3d): rate %.5f +/- %.5f, gated %.5f, mean w %.4f" % (G, int(sz[0]), r["p"], r["se"], r["p_sig"], r["mean_w"])) print("\n Total sample size, holding G = 4 and c = 1:") ntab = {} for Nv in (40, 80, 240, 800, 2400): sz = group_sizes(Nv, 4, 1) r = run_cell(1.0, sz, "random", TRIALS_GRID, sens_seeds[si]) si += 1 ntab[Nv] = r print(" N = %5d: rate %.5f +/- %.5f, gated %.5f +/- %.5f" % (Nv, r["p"], r["se"], r["p_sig"], r["se_sig"])) print(""" The point-estimate reversal rate is nearly flat in N, because it is set by the population geometry rather than by sampling error once N is moderate. The SIGNIFICANCE-GATED rate is not flat at all. More data does not protect you. It converts a reversal you might have dismissed as noise into a reversal with confidence intervals on opposite sides of zero.""") # =========================================================================== rule("11. FIGURE DATA") # =========================================================================== x, y, gi, d, within, pooled = constructed["flagship 5 groups"] print("\n FIG1: flagship constructed dataset, 5 groups of 15 points.") print(" within slope +1.000000, pooled slope %.6f" % pooled) print(" group means (x, y):") for g in range(5): m = gi == g print(" GM %d %.4f %.4f %d" % (g + 1, x[m].mean(), y[m].mean(), int(m.sum()))) print(" points, as 'PT group x y':") for g in range(5): for i in np.where(gi == g)[0]: print(" PT %d %.4f %.4f" % (g + 1, x[i], y[i])) print(" pooled fit: y = %.6f + %.6f x" % (ap_f, bp_f)) print(" between-group line slope b_b = %.6f" % d["b_between"]) print(" NOISY PANEL (sigma = 0.8), fitted pooled slope %.5f," % bpn) print(" fitted within slopes %s" % " ".join("%.4f" % v for v in wsl)) print(" noisy points, as 'NPT group x y':") for g in range(5): for i in np.where(gin == g)[0]: print(" NPT %d %.4f %.4f" % (g + 1, xn[i], yn[i])) print(" noisy pooled fit: y = %.6f + %.6f x" % (apn, bpn)) for g in range(5): m = gin == g bq, aq = ols(xn[m], yn[m])[0], ols(xn[m], yn[m])[1] print(" NFIT %d %.6f %.6f %.4f %.4f" % (g + 1, bq, aq, xn[m].min(), xn[m].max())) print("\n FIG3b: categorical allocation gap h against reversal probability.") print(" HGRID : %s" % " ".join("%.2f" % v for v in H_GRID)) print(" CATANA : %s" % " ".join("%.6f" % cat[h]["analytic"] for h in H_GRID)) print(" CATPOP : %s" % " ".join("%.6f" % cat[h]["p_pop"] for h in H_GRID)) print(" CATOBS : %s" % " ".join("%.6f" % cat[h]["p_obs"] for h in H_GRID)) print(" CATOBSSE : %s" % " ".join("%.6f" % cat[h]["se_obs"] for h in H_GRID)) print(" kidney-stone gap h = 0.523, analytic P = %.6f" % analytic_cat(0.523)) print("\n FIG2: the lever. Pooled slope as a function of the between-group") print(" weight 1-w, for b_w = +1 and three between-group slopes.") LEVER_X = [0.0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.85, 1.0] print(" LEVERX : %s" % " ".join("%.2f" % v for v in LEVER_X)) for bb in (-1.0, -3.0, -8.0): row = " LEVER %.1f :" % bb for wf in LEVER_X: row += " %.4f" % ((1 - wf) * 1.0 + wf * bb) print(row) print(" flagship dataset sits at 1-w = %.5f, b_b = %.4f" % (1 - d["w"], d["b_between"])) print("\n FIG3: reversal frequency against confounder strength.") print(" CGRID : %s" % " ".join("%.2f" % v for v in C_GRID)) for mode in ("random", "coupled"): for r in R_GRID: row = " CURVE %s %d:1 :" % (mode, r) for c in C_GRID: row += " %.5f" % grid[(mode, r, c)]["p"] print(row) row = " CURVE gated random 1:1 :" for c in C_GRID: row += " %.5f" % grid[("random", 1, c)]["p_sig"] print(row) row = " CURVE se random 1:1 :" for c in C_GRID: row += " %.5f" % grid[("random", 1, c)]["se"] print(row) print("\n FIG5: continuous case.") print(" CONTX : %s" % " ".join("%.2f" % v for v in S_LIST)) for key, lab in (("analytic", "ANA"), ("p_pop", "POP"), ("p_samp", "SAMP")): print(" CONT%-5s : %s" % (lab, " ".join("%.6f" % cont[s][key] for s in S_LIST))) print(" CONTSE : %s" % " ".join("%.6f" % cont[s]["se_samp"] for s in S_LIST)) S_FINE = [0.1 + 0.2 * i for i in range(40)] print(" CONTFINEX : %s" % " ".join("%.2f" % v for v in S_FINE)) print(" CONTFINEY : %s" % " ".join("%.6f" % analytic_reversal(v, nseg=20000) for v in S_FINE)) H_FINE = [0.0 + 0.02 * i for i in range(46)] print(" CATFINEX : %s" % " ".join("%.2f" % v for v in H_FINE)) print(" CATFINEY : %s" % " ".join("%.6f" % analytic_cat(v) for v in H_FINE)) # =========================================================================== rule("12. SUMMARY TABLE: club value beside the accepted value") # =========================================================================== print() print(" %-36s %18s %20s" % ("quantity", "club value", "reference")) print(" " + "-" * 74) rows = [ ("OLS slope, 5-point hand dataset", "%.12f" % cb, "3/5 exact"), ("OLS intercept, same dataset", "%.12f" % ca, "11/5 exact"), ("t_{0.975, df=10}", "%.4f" % t_crit(10), "2.228 (table)"), ("Weighted-average identity residual", "%.2e" % abs(b_pool - pred), "0 exact"), ("Constructed within slope (flagship)", "%.12f" % within[0], "+1 by design"), ("Constructed pooled slope (flagship)", "%.12f" % pooled, "-1 by design"), ("Worst deviation, 8 constructions", "%.2e" % max_dev_all, "0 exact"), ("Control, c=0 balanced", "%d in %d" % (n0, TRIALS_CTL), "0 exact"), ("Control, c=0 imbalance 30:1", "%d in %d" % (ctl_rows[1][1]["n_rev"], TRIALS_CTL), "0 exact"), ("Reversal rate, c=1 balanced", "%.5f" % hl["p"], "no published value"), ("Gated reversal rate, c=1 balanced", "%.5f" % hl["p_sig"], "no published value"), ("Reversal rate, c=8 balanced", "%.5f" % grid[("random", 1, 8.0)]["p"], "ceiling 0.5"), ("2x2x2 table, directional", "%.6f" % p_dir, "1/60 = %.6f" % PP), ("2x2x2 table, either direction", "%.6f" % p_any, "2/60 = %.6f" % (2 * PP)), ("Categorical, allocation gap h=0", "%.6f" % cat[0.0]["p_pop"], "0 exact"), ("Categorical, allocation gap h=0.5", "%.6f" % cat[0.5]["p_pop"], "%.6f closed form" % cat[0.5]["analytic"]), ("Continuous s=1, population MC", "%.6f" % cont[1.0]["p_pop"], "%.6f integral" % ana[1.0]), ("Continuous s=1, n=200 sample MC", "%.6f" % cont[1.0]["p_samp"], "no closed form"), ("Berkeley pooled gap (women-men)", "%+.1f pp" % (100 * (rw - rm)), "published table"), ("Berkeley standardised gap", "%+.1f pp" % (100 * (std_w - std_m)), "published table"), ("Kidney stones, pooled gap", "%+.1f pp" % (100 * (oa / ta - ob / tb)), "published table"), ("Kidney stones, small-stone gap", "%+.1f pp" % (100 * (81 / 87 - 234 / 270)), "published table"), ] for a_, b_, c_ in rows: print(" %-36s %18s %20s" % (a_, b_, c_)) rule() print("Wall clock: %.1f s" % (time.time() - T0)) print("Seed: %d. Rerunning this file reproduces every number above." % MASTER_SEED)