============================================================================== SIMPSON'S PARADOX: CONSTRUCTION, FREQUENCY, AND THE CONTINUOUS CASE ============================================================================== Science Journaling Club, Volume 2 Issue 2, Winter 2026 Master seed: 20260214 | numpy 2.4.2 | Python 3.12.3 No data were collected. The computation is the experiment. ============================================================================== 1. VALIDATION: the regression routine against a closed-form solution ============================================================================== 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. exact rational xbar = 3, ybar = 4 exact rational Sxx = 10, Sxy = 6 exact rational slope = 3/5 = 0.59999999999999998 exact rational intercept = 11/5 = 2.2000000000000002 slope, club ols() vs exact rational club 0.6 ref 0.6 diff 0 [OK] intercept, club ols() vs exact rational club 2.2 ref 2.2 diff 0 [OK] slope, club ols() vs numpy.polyfit club 0.6 ref 0.6 diff 0 [OK] intercept, club ols() vs numpy.polyfit club 2.2 ref 2.2 diff 1.332268e-15 [OK] slope, club ols() vs numpy.linalg.lstsq club 0.6 ref 0.6 diff -5.551115e-16 [OK] residual sum of squares, exact: Syy - b*Sxy = 12/5 = 2.4 se(slope), club vs closed form club 0.2828427125 ref 0.2828427125 diff 5.551115e-17 [OK] club routine also reports t = 2.121320 on df = 3, r = 0.7745966692 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. t_{0.975, df=1} club 12.70620474 ref 12.706 diff 0.0002047361747 [OK] t_{0.975, df=2} club 4.30265273 ref 4.303 diff -0.0003472702505 [OK] t_{0.975, df=5} club 2.570581836 ref 2.571 diff -0.0004181643637 [OK] t_{0.975, df=10} club 2.228138852 ref 2.228 diff 0.0001388519863 [OK] t_{0.975, df=30} club 2.042272456 ref 2.042 diff 0.0002724563012 [OK] t_{0.975, df=60} club 2.000297822 ref 2 diff 0.0002978220143 [OK] t_{0.975, df=100} club 1.983971519 ref 1.984 diff -2.848147646e-05 [OK] t_{0.975, df=1000} club 1.962339081 ref 1.962 diff 0.0003390808264 [OK] Differences are at the rounding of the published tables (3 d.p.). ============================================================================== 2. VALIDATION: the weighted-average identity, to machine precision ============================================================================== 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. group sizes [37, 12, 61, 25, 9] (N = 144) T_w = 508.4379356650 T_b = 1560.1534981725 w = 0.245789442685 effective within slope b_w = 0.015471237425 between-group slope b_b = -0.436640945469 pooled slope, fitted vs identity club -0.325516544 ref -0.325516544 diff 1.110223e-16 [OK] relative difference: 3.411e-16 The identity is algebra, so the only thing this can catch is a coding error. It caught two while we were writing the file. ============================================================================== 3. THE CONSTRUCTOR ============================================================================== 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. dataset G b_w tgt b_p tgt w b_between max|dev| ------------------------------------------------------------------------------ flagship 5 groups 5 1.000 -1.000 0.12500 -1.2857 2.22e-16 mild reversal 4 0.500 -0.100 0.17373 -0.2262 2.78e-17 violent reversal 3 0.500 -2.000 0.28073 -2.9758 4.44e-16 downhill groups, uphill whole 3 -1.000 1.000 0.15122 1.3563 0.00e+00 nearly flat both ways 4 0.100 -0.100 0.05138 -0.1108 1.39e-17 wildly unequal sizes 4 2.000 -0.250 0.06052 -0.3949 4.44e-16 two groups only 2 1.000 -1.000 0.10227 -1.2278 0.00e+00 ten groups 10 0.800 -0.400 0.03457 -0.4430 2.22e-16 ------------------------------------------------------------------------------ worst deviation over all 8 constructions club 4.440892e-16 ref 0 diff 4.440892e-16 [OK] The construction is algebra, so this is a check on the code. The idea is checked by the fits themselves, printed next. The flagship dataset in full. Five groups, within-group slope +1.000 by design, pooled slope -1.000 by design. This is Figure 1. T_w = 85.714286 T_b = 600.000000 w = 0.12500000 required between-group slope b_b = -1.28571429 group | n | xbar | ybar | fitted within slope | target 1 | 15 | -4.000 | 5.143 | 1.000000000000 | +1.0 2 | 15 | -2.000 | 2.571 | 1.000000000000 | +1.0 3 | 15 | 0.000 | -0.000 | 1.000000000000 | +1.0 4 | 15 | 2.000 | -2.571 | 1.000000000000 | +1.0 5 | 15 | 4.000 | -5.143 | 1.000000000000 | +1.0 pooled| 75 | -0.000 | -0.000 | -1.000000000000 | -1.0 flagship within slope (all five identical) club 1 ref 1 diff 2.220446e-16 [OK] flagship pooled slope club -1 ref -1 diff -2.220446e-16 [OK] pooled correlation r = -0.797724; within-group correlation r = +1.000000 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. group | fitted | se | z 1 | 0.92432 | 0.28742 | -0.26 2 | 1.07493 | 0.18668 | +0.40 3 | 0.55266 | 0.42537 | -1.05 4 | 0.81004 | 0.27179 | -0.70 5 | 0.82793 | 0.19752 | -0.87 pooled| -1.00911 | 0.08309 | -0.11 largest |z| across the six fits: 1.05 (about 1 is expected; above 3 would mean the constructor is broken) all five within-group slopes positive: True pooled slope negative: True 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: b_p target | required b_between -0.10 | -0.2571 -0.50 | -0.7143 -1.00 | -1.2857 -5.00 | -5.8571 -20.00 | -23.0000 w = 0.12500 for this layout, so the lever arm 1/(1-w) = 1.143 ============================================================================== 4. THE CONTROL: remove the confounder, balance the groups ============================================================================== 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. balanced sizes : [np.int64(60), np.int64(60), np.int64(60), np.int64(60)] (sum 240) imbalanced sizes: [np.int64(5), np.int64(17), np.int64(53), np.int64(165)] (sum 240) cell reversals rate SE ------------------------------------------------------------------------------ no confounder (c=0), balanced [CONTROL] 0 0.000000 0.000000 no confounder (c=0), imbalance 30:1 0 0.000000 0.000000 confounder c=1, balanced 3795 0.094875 0.001465 confounder c=1, imbalance 30:1 1763 0.044075 0.001026 ------------------------------------------------------------------------------ The control cell produced 0 reversals in 40000 datasets. control reversal rate vs analytic 0 club 0 ref 0 diff 0 [OK] Zero events in 40000 trials. By the rule of three the one-sided 95% upper bound on the true rate is 3/40000 = 7.50e-05. 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. ============================================================================== 5. HOW OFTEN BY ACCIDENT: the confounder / imbalance grid ============================================================================== 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. mode c ratio rate SE gated SE mean w ------------------------------------------------------------------------------ random 0.00 1:1 0.00000 0.00000 0.00000 0.00000 0.9875 random 0.25 1:1 0.00000 0.00000 0.00000 0.00000 0.9442 random 0.50 1:1 0.00340 0.00041 0.00100 0.00022 0.8442 random 1.00 1:1 0.09695 0.00209 0.06040 0.00168 0.6221 random 2.00 1:1 0.28850 0.00320 0.24090 0.00302 0.3415 random 4.00 1:1 0.42690 0.00350 0.39560 0.00346 0.1382 random 8.00 1:1 0.48360 0.00353 0.46915 0.00353 0.0457 random 0.00 3:1 0.00000 0.00000 0.00000 0.00000 0.9874 random 0.25 3:1 0.00000 0.00000 0.00000 0.00000 0.9461 random 0.50 3:1 0.00455 0.00048 0.00140 0.00026 0.8504 random 1.00 3:1 0.08700 0.00199 0.05485 0.00161 0.6345 random 2.00 3:1 0.28105 0.00318 0.23430 0.00300 0.3573 random 4.00 3:1 0.42240 0.00349 0.38910 0.00345 0.1466 random 8.00 3:1 0.48410 0.00353 0.46725 0.00353 0.0490 random 0.00 10:1 0.00000 0.00000 0.00000 0.00000 0.9876 random 0.25 10:1 0.00000 0.00000 0.00000 0.00000 0.9528 random 0.50 10:1 0.00335 0.00041 0.00090 0.00021 0.8701 random 1.00 10:1 0.06795 0.00178 0.03430 0.00129 0.6770 random 2.00 10:1 0.24370 0.00304 0.16300 0.00261 0.4041 random 4.00 10:1 0.40310 0.00347 0.30445 0.00325 0.1802 random 8.00 10:1 0.46995 0.00353 0.37490 0.00342 0.0624 random 0.00 30:1 0.00000 0.00000 0.00000 0.00000 0.9873 random 0.25 30:1 0.00000 0.00000 0.00000 0.00000 0.9603 random 0.50 30:1 0.00175 0.00030 0.00010 0.00007 0.8925 random 1.00 30:1 0.04855 0.00152 0.00820 0.00064 0.7253 random 2.00 30:1 0.19285 0.00279 0.04245 0.00143 0.4668 random 4.00 30:1 0.34935 0.00337 0.08515 0.00197 0.2268 random 8.00 30:1 0.43295 0.00350 0.11345 0.00224 0.0854 coupled 0.00 1:1 0.00000 0.00000 0.00000 0.00000 0.9875 coupled 0.25 1:1 0.00000 0.00000 0.00000 0.00000 0.9444 coupled 0.50 1:1 0.00350 0.00042 0.00080 0.00020 0.8447 coupled 1.00 1:1 0.09395 0.00206 0.05970 0.00168 0.6216 coupled 2.00 1:1 0.29405 0.00322 0.24820 0.00305 0.3410 coupled 4.00 1:1 0.42845 0.00350 0.39910 0.00346 0.1381 coupled 8.00 1:1 0.47665 0.00353 0.46135 0.00352 0.0453 coupled 0.00 3:1 0.00000 0.00000 0.00000 0.00000 0.9874 coupled 0.25 3:1 0.00000 0.00000 0.00000 0.00000 0.9484 coupled 0.50 3:1 0.00325 0.00040 0.00105 0.00023 0.8546 coupled 1.00 3:1 0.08150 0.00193 0.04905 0.00153 0.6434 coupled 2.00 3:1 0.28045 0.00318 0.23135 0.00298 0.3614 coupled 4.00 3:1 0.41800 0.00349 0.38505 0.00344 0.1500 coupled 8.00 3:1 0.47305 0.00353 0.45540 0.00352 0.0513 coupled 0.00 10:1 0.00000 0.00000 0.00000 0.00000 0.9876 coupled 0.25 10:1 0.00000 0.00000 0.00000 0.00000 0.9585 coupled 0.50 10:1 0.00120 0.00024 0.00025 0.00011 0.8882 coupled 1.00 10:1 0.04950 0.00153 0.02435 0.00109 0.7080 coupled 2.00 10:1 0.22690 0.00296 0.14640 0.00250 0.4315 coupled 4.00 10:1 0.38235 0.00344 0.28700 0.00320 0.1965 coupled 8.00 10:1 0.47235 0.00353 0.37540 0.00342 0.0705 coupled 0.00 30:1 0.00000 0.00000 0.00000 0.00000 0.9874 coupled 0.25 30:1 0.00000 0.00000 0.00000 0.00000 0.9675 coupled 0.50 30:1 0.00040 0.00014 0.00000 0.00000 0.9148 coupled 1.00 30:1 0.02425 0.00109 0.00315 0.00040 0.7756 coupled 2.00 30:1 0.14905 0.00252 0.03000 0.00121 0.5278 coupled 4.00 30:1 0.31475 0.00328 0.07470 0.00186 0.2670 coupled 8.00 30:1 0.41925 0.00349 0.10725 0.00219 0.1033 ------------------------------------------------------------------------------ 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. c = 0.25 rate = 0.0000 +/- 0.0000 gated = 0.0000 +/- 0.0000 c = 0.50 rate = 0.0034 +/- 0.0004 gated = 0.0010 +/- 0.0002 c = 1.00 rate = 0.0969 +/- 0.0021 gated = 0.0604 +/- 0.0017 c = 2.00 rate = 0.2885 +/- 0.0032 gated = 0.2409 +/- 0.0030 c = 4.00 rate = 0.4269 +/- 0.0035 gated = 0.3956 +/- 0.0035 c = 8.00 rate = 0.4836 +/- 0.0035 gated = 0.4692 +/- 0.0035 The second reading surprised us. Compare the 'random' rows across imbalance ratios at fixed confounder strength: c = 0.50 : 1:1 -> 0.0034 3:1 -> 0.0046 10:1 -> 0.0034 30:1 -> 0.0018 c = 1.00 : 1:1 -> 0.0969 3:1 -> 0.0870 10:1 -> 0.0679 30:1 -> 0.0486 c = 2.00 : 1:1 -> 0.2885 3:1 -> 0.2811 10:1 -> 0.2437 30:1 -> 0.1928 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: c = 0.50 : 1:1 -> 0.0035 3:1 -> 0.0032 10:1 -> 0.0012 30:1 -> 0.0004 c = 1.00 : 1:1 -> 0.0940 3:1 -> 0.0815 10:1 -> 0.0495 30:1 -> 0.0243 c = 2.00 : 1:1 -> 0.2940 3:1 -> 0.2804 10:1 -> 0.2269 30:1 -> 0.1490 difference (coupled minus random), with the SE of the difference: c=0.50, 10:1 : -0.00215 +/- 0.00048 (-4.51 SE) c=0.50, 30:1 : -0.00135 +/- 0.00033 (-4.12 SE) c=1.00, 10:1 : -0.01845 +/- 0.00235 (-7.85 SE) c=1.00, 30:1 : -0.02430 +/- 0.00187 (-13.00 SE) c=2.00, 10:1 : -0.01680 +/- 0.00424 (-3.96 SE) c=2.00, 30:1 : -0.04380 +/- 0.00376 (-11.65 SE) 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. HEADLINE. With a confounder as strong as the within-group spread (c = 1), four balanced groups and 240 points, 9.70% of datasets showed a full sign reversal, and 6.04% 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. plain: 1939 reversals / 20000 trials, SE 0.00209 gated: 1208 reversals / 20000 trials, SE 0.00168 ============================================================================== 5b. THE CATEGORICAL FORM: where imbalance really is the culprit ============================================================================== 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. h MC pop SE closed form diff z observable SE ------------------------------------------------------------------------------ 0.00 0.000000 0.000000 0.000000 +0.000000 nan 0.000000 0.000000 0.10 0.134925 0.000764 0.133333 +0.001592 2.08 0.052785 0.000500 0.20 0.252830 0.000972 0.252604 +0.000226 0.23 0.134345 0.000763 0.30 0.320650 0.001044 0.320602 +0.000048 0.05 0.194405 0.000885 0.50 0.385755 0.001088 0.385417 +0.000338 0.31 0.247950 0.000966 0.70 0.414475 0.001102 0.416029 -0.001554 -1.41 0.265570 0.000988 0.90 0.434855 0.001109 0.433771 +0.001084 0.98 0.259335 0.000980 ------------------------------------------------------------------------------ largest |z|, categorical MC vs closed form club 2.083502495 ref 0 diff 2.083502495 [OK] h = 0 gives exactly 0 reversals in 200000 trials, which is the categorical control: equal allocation, no paradox, ever. 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 0.39 for a treatment whose real advantage is somewhere in 0 to 10 points. ============================================================================== 6. THE CONTINUOUS CASE: one lurking variable, no groups at all ============================================================================== 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. Analytic reversal probability, with a coarser Simpson grid as a convergence check on the quadrature itself: s P (200k seg) P (50k seg) diff 0.25 0.000000188 0.000000188 1.77e-19 0.50 0.005540270 0.005540270 -3.14e-16 1.00 0.097320986 0.097320986 -7.01e-15 1.50 0.181655914 0.181655914 -9.46e-15 2.00 0.232655251 0.232655251 -1.51e-14 3.00 0.284345338 0.284345338 -1.33e-14 5.00 0.321256463 0.321256463 -2.32e-14 8.00 0.337878621 0.337878621 -4.88e-12 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. s MC pop SE analytic diff z MC samp SE ------------------------------------------------------------------------------ 0.25 0.000000 0.000000 0.000000 -0.000000 nan 0.000000 0.000000 0.50 0.005647 0.000118 0.005540 +0.000107 0.90 0.008467 0.000529 1.00 0.097863 0.000470 0.097321 +0.000542 1.15 0.097800 0.001715 1.50 0.181820 0.000610 0.181656 +0.000164 0.27 0.179467 0.002216 2.00 0.233480 0.000669 0.232655 +0.000825 1.23 0.232867 0.002440 3.00 0.285490 0.000714 0.284345 +0.001145 1.60 0.285967 0.002609 5.00 0.321078 0.000738 0.321256 -0.000179 -0.24 0.323200 0.002700 8.00 0.337750 0.000748 0.337879 -0.000129 -0.17 0.339500 0.002734 ------------------------------------------------------------------------------ largest |z|, Monte Carlo vs analytic integral club 1.602904483 ref 0 diff 1.602904483 [OK] Every cell with a non-zero count agrees with the closed form inside 3 standard errors. The s = 0.25 cell produced zero events, so it has no usable normal standard error. The integral predicts 0.0751 events in 400000 draws, and the Poisson probability of seeing none of them is 0.9276, so zero is the expected outcome rather than a failure. 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. s = 0.50 : sample 0.00847, population 0.00565, excess +0.00282 (5.2 SE) s = 1.00 : sample 0.09780, population 0.09786, excess -0.00006 (-0.0 SE) s = 2.00 : sample 0.23287, population 0.23348, excess -0.00061 (-0.2 SE) s = 5.00 : sample 0.32320, population 0.32108, excess +0.00212 (0.8 SE) ============================================================================== 7. A PUBLISHED NUMBER TO CHECK AGAINST: random 2x2x2 tables ============================================================================== 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: 3000000 one-directional paradox (A better in both strata, worse pooled): count 25113, rate 0.0083710, SE 0.0000526 either direction: count 50201, rate 0.0167337, SE 0.0000741 directional rate vs published 1/60 club 0.008371 ref 0.01666666667 diff -0.008295666667 [MISMATCH] -157.71 standard errors from 1/60 either-direction rate vs published 1/60 club 0.01673366667 ref 0.01666666667 diff 6.7e-05 [OK] +0.90 standard errors from 1/60 either-direction rate vs 2/60 club 0.01673366667 ref 0.03333333333 diff -0.01659966667 [MISMATCH] -224.14 standard errors from 2/60 directional rate vs 1/120 club 0.008371 ref 0.008333333333 diff 3.766666667e-05 [OK] +0.72 standard errors from 1/120 Two of those four lines are MEANT to fail. They are printed so that the reader can see which convention is ruled out rather than being shown only the one that agrees. 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. ============================================================================== 8. CONVERGENCE ============================================================================== 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. The quantity tracked is the ONE-DIRECTIONAL rate, whose exact value is 1/120 = 0.0083333, half the published 1/60 because 1/60 counts both directions. The z column is against 1/120. trials estimate SE 1/120 z 1000 0.0080000 0.0028171 0.0083333 -0.12 3162 0.0110689 0.0018606 0.0083333 +1.47 10000 0.0103000 0.0010096 0.0083333 +1.95 31623 0.0086962 0.0005221 0.0083333 +0.69 100000 0.0082900 0.0002867 0.0083333 -0.15 316228 0.0084781 0.0001630 0.0083333 +0.89 1000000 0.0082670 0.0000905 0.0083333 -0.73 3000000 0.0083710 0.0000526 0.0083333 +0.72 FIGDATA convergence (trials, estimate, se): CONV 1000 0.00800000 0.00281709 CONV 1334 0.00974513 0.00268961 CONV 1778 0.01012373 0.00237408 CONV 2371 0.01096584 0.00213875 CONV 3162 0.01106894 0.00186061 CONV 4217 0.01138250 0.00163355 CONV 5623 0.01084830 0.00138143 CONV 7499 0.01080144 0.00119366 CONV 10000 0.01030000 0.00100965 CONV 13335 0.00989876 0.00085730 CONV 17783 0.00944722 0.00072542 CONV 23714 0.00902421 0.00061409 CONV 31623 0.00869620 0.00052212 CONV 42170 0.00844202 0.00044553 CONV 56234 0.00821567 0.00038065 CONV 74989 0.00813453 0.00032802 CONV 100000 0.00829000 0.00028673 CONV 133352 0.00831634 0.00024869 CONV 177828 0.00841825 0.00021666 CONV 237137 0.00848033 0.00018830 CONV 316228 0.00847806 0.00016304 CONV 421697 0.00846579 0.00014109 CONV 562341 0.00831880 0.00012112 CONV 749894 0.00829717 0.00010475 CONV 1000000 0.00826700 0.00009055 CONV 1333521 0.00834932 0.00007880 CONV 1778279 0.00836427 0.00006830 CONV 2371374 0.00835929 0.00005912 CONV 3000000 0.00837100 0.00005260 ============================================================================== 9. TWO REAL PUBLISHED TABLES, RECOMPUTED ============================================================================== 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. (a) Berkeley graduate admissions, autumn 1973, six largest departments. Transcribed from Bickel, Hammel & O'Connell (1975), Science 187, 398-404. dept men apps men % women apps women % women-men (pp) A 825 62 108 82 +20 B 560 63 25 68 +5 C 325 37 593 34 -3 D 417 33 375 35 +2 E 191 28 393 24 -4 F 373 6 341 7 +1 pooled: men 1198/2691 = 0.4452 (44.5%), women 557/1835 = 0.3033 (30.3%) pooled difference (women minus men): -14.2 percentage points departments admitting women at the higher rate: 4 of 6 directly standardised to the combined applicant pool: men 0.3876 (38.8%), women 0.4291 (42.9%), difference +4.2 pp SIGN FLIP: pooled -14.2 pp, standardised +4.2 pp. The campus-wide figures quoted in the same paper are 8442 male applicants at 44%, 4321 female applicants at 35%. (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). stratum open surgery PCNL difference small (<2 cm) 81/87 93.1% 234/270 86.7% +6.4 pp large (>=2 cm) 192/263 73.0% 55/80 68.8% +4.3 pp pooled 273/350 78.0% 289/350 82.6% -4.6 pp SIGN FLIP: open surgery wins both strata by +6.4 and +4.3 pp and loses the pooled comparison by 4.6 pp. The confounder is stone size: 75.1% of the open-surgery cases were large stones against 22.9% of the PCNL cases, a gap of 52.3 percentage points in exposure to the harder problem. ============================================================================== 10. WHERE A DIFFERENT MODELLING CHOICE WOULD CHANGE THE ANSWER ============================================================================== Six knobs we set by hand, each re-run at c = 1 with four balanced groups, N = 240 and 20000 trials, changing one thing at a time. variant rate SE vs base ------------------------------------------------------------------------------ BASE: G=4, N=240, sigma=1, tau=1, b_w=1 0.09695 0.00209 - theta SD tau = 0.5 (weak y-loading) 0.02890 0.00118 -0.06805 (-28.3 SE) theta SD tau = 2.0 (strong y-loading) 0.19270 0.00279 +0.09575 (+27.5 SE) residual sigma = 0.25 (tight groups) 0.09240 0.00205 -0.00455 (-1.6 SE) residual sigma = 4.0 (noisy groups) 0.09005 0.00202 -0.00690 (-2.4 SE) within slope b_w = 0.25 (weak signal) 0.25680 0.00309 +0.15985 (+42.8 SE) within slope b_w = 4.0 (strong signal) 0.00405 0.00045 -0.09290 (-43.4 SE) Number of groups, holding N = 240 and c = 1: G = 2 (n_g = 120): rate 0.05620 +/- 0.00163, gated 0.03835, mean w 0.7526 G = 3 (n_g = 80): rate 0.07720 +/- 0.00189, gated 0.04995, mean w 0.6670 G = 4 (n_g = 60): rate 0.09305 +/- 0.00205, gated 0.05995, mean w 0.6216 G = 6 (n_g = 40): rate 0.10735 +/- 0.00219, gated 0.06550, mean w 0.5745 G = 8 (n_g = 30): rate 0.11735 +/- 0.00228, gated 0.06850, mean w 0.5497 G = 12 (n_g = 20): rate 0.13340 +/- 0.00240, gated 0.04655, mean w 0.5181 Total sample size, holding G = 4 and c = 1: N = 40: rate 0.10340 +/- 0.00215, gated 0.00435 +/- 0.00047 N = 80: rate 0.09750 +/- 0.00210, gated 0.03490 +/- 0.00130 N = 240: rate 0.09245 +/- 0.00205, gated 0.05850 +/- 0.00166 N = 800: rate 0.09215 +/- 0.00205, gated 0.07365 +/- 0.00185 N = 2400: rate 0.08980 +/- 0.00202, gated 0.07805 +/- 0.00190 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. ============================================================================== 11. FIGURE DATA ============================================================================== FIG1: flagship constructed dataset, 5 groups of 15 points. within slope +1.000000, pooled slope -1.000000 group means (x, y): GM 1 -4.0000 5.1429 15 GM 2 -2.0000 2.5714 15 GM 3 0.0000 -0.0000 15 GM 4 2.0000 -2.5714 15 GM 5 4.0000 -5.1429 15 points, as 'PT group x y': PT 1 -5.7321 3.4108 PT 1 -5.4846 3.6582 PT 1 -5.2372 3.9057 PT 1 -4.9897 4.1531 PT 1 -4.7423 4.4005 PT 1 -4.4949 4.6480 PT 1 -4.2474 4.8954 PT 1 -4.0000 5.1429 PT 1 -3.7526 5.3903 PT 1 -3.5051 5.6377 PT 1 -3.2577 5.8852 PT 1 -3.0103 6.1326 PT 1 -2.7628 6.3800 PT 1 -2.5154 6.6275 PT 1 -2.2679 6.8749 PT 2 -3.7321 0.8394 PT 2 -3.4846 1.0868 PT 2 -3.2372 1.3342 PT 2 -2.9897 1.5817 PT 2 -2.7423 1.8291 PT 2 -2.4949 2.0766 PT 2 -2.2474 2.3240 PT 2 -2.0000 2.5714 PT 2 -1.7526 2.8189 PT 2 -1.5051 3.0663 PT 2 -1.2577 3.3137 PT 2 -1.0103 3.5612 PT 2 -0.7628 3.8086 PT 2 -0.5154 4.0560 PT 2 -0.2679 4.3035 PT 3 -1.7321 -1.7321 PT 3 -1.4846 -1.4846 PT 3 -1.2372 -1.2372 PT 3 -0.9897 -0.9897 PT 3 -0.7423 -0.7423 PT 3 -0.4949 -0.4949 PT 3 -0.2474 -0.2474 PT 3 0.0000 -0.0000 PT 3 0.2474 0.2474 PT 3 0.4949 0.4949 PT 3 0.7423 0.7423 PT 3 0.9897 0.9897 PT 3 1.2372 1.2372 PT 3 1.4846 1.4846 PT 3 1.7321 1.7321 PT 4 0.2679 -4.3035 PT 4 0.5154 -4.0560 PT 4 0.7628 -3.8086 PT 4 1.0103 -3.5612 PT 4 1.2577 -3.3137 PT 4 1.5051 -3.0663 PT 4 1.7526 -2.8189 PT 4 2.0000 -2.5714 PT 4 2.2474 -2.3240 PT 4 2.4949 -2.0766 PT 4 2.7423 -1.8291 PT 4 2.9897 -1.5817 PT 4 3.2372 -1.3342 PT 4 3.4846 -1.0868 PT 4 3.7321 -0.8394 PT 5 2.2679 -6.8749 PT 5 2.5154 -6.6275 PT 5 2.7628 -6.3800 PT 5 3.0103 -6.1326 PT 5 3.2577 -5.8852 PT 5 3.5051 -5.6377 PT 5 3.7526 -5.3903 PT 5 4.0000 -5.1429 PT 5 4.2474 -4.8954 PT 5 4.4949 -4.6480 PT 5 4.7423 -4.4005 PT 5 4.9897 -4.1531 PT 5 5.2372 -3.9057 PT 5 5.4846 -3.6582 PT 5 5.7321 -3.4108 pooled fit: y = -0.000000 + -1.000000 x between-group line slope b_b = -1.285714 NOISY PANEL (sigma = 0.8), fitted pooled slope -1.00911, fitted within slopes 0.9243 1.0749 0.5527 0.8100 0.8279 noisy points, as 'NPT group x y': NPT 1 -3.2420 5.1776 NPT 1 -4.5732 3.5406 NPT 1 -3.3942 4.3489 NPT 1 -3.4738 4.5215 NPT 1 -3.8714 3.3606 NPT 1 -3.6735 4.5064 NPT 1 -5.9373 3.2944 NPT 1 -3.5564 6.8797 NPT 1 -2.9335 7.4760 NPT 1 -3.4939 6.1837 NPT 1 -5.4931 4.1041 NPT 1 -5.1927 4.4092 NPT 1 -2.8081 6.4599 NPT 1 -4.2127 4.3004 NPT 1 -4.1444 5.2042 NPT 2 -4.6165 -0.5096 NPT 2 -3.1832 0.4795 NPT 2 -0.9086 3.2897 NPT 2 -2.1429 1.3341 NPT 2 -1.7102 3.4455 NPT 2 -2.0469 2.4070 NPT 2 -3.2672 2.7551 NPT 2 -1.7560 1.0185 NPT 2 -0.5578 3.7380 NPT 2 -1.6574 3.9595 NPT 2 -3.9149 0.6780 NPT 2 -0.2785 4.4868 NPT 2 -0.2910 5.2881 NPT 2 -1.9888 2.1875 NPT 2 -1.6803 2.3257 NPT 3 -0.1305 0.1897 NPT 3 -0.8730 0.2127 NPT 3 0.0593 0.7082 NPT 3 -0.0637 1.2023 NPT 3 -0.0688 -0.1829 NPT 3 0.5948 -0.1874 NPT 3 0.3655 0.0359 NPT 3 -0.0675 0.1945 NPT 3 0.1785 0.4421 NPT 3 -0.4754 -1.0020 NPT 3 0.6894 -0.3047 NPT 3 -0.1467 0.1077 NPT 3 -0.6242 -2.1285 NPT 3 -0.2388 0.1510 NPT 3 0.8010 0.5614 NPT 4 1.0199 -4.1963 NPT 4 2.3157 -1.8633 NPT 4 1.8074 -1.6783 NPT 4 2.7839 -1.4702 NPT 4 3.1633 -0.5119 NPT 4 2.2696 -1.9078 NPT 4 1.1098 -3.0582 NPT 4 2.5428 -2.9145 NPT 4 3.5784 -1.4409 NPT 4 2.2636 -2.4148 NPT 4 -0.0192 -4.6079 NPT 4 3.0325 -2.6612 NPT 4 1.2536 -4.9542 NPT 4 0.6748 -1.4088 NPT 4 2.2039 -1.7953 NPT 5 4.2763 -5.6523 NPT 5 4.4615 -4.7444 NPT 5 3.4117 -5.4412 NPT 5 4.8111 -3.6368 NPT 5 3.1402 -5.2850 NPT 5 4.2002 -6.1456 NPT 5 3.9965 -4.5919 NPT 5 4.3418 -6.2152 NPT 5 4.3617 -3.9095 NPT 5 3.1712 -5.3202 NPT 5 5.7208 -3.0430 NPT 5 2.8716 -5.1357 NPT 5 5.4307 -2.8183 NPT 5 1.3999 -6.9832 NPT 5 4.4049 -4.8448 noisy pooled fit: y = -0.000000 + -1.009110 x NFIT 1 0.924317 8.615065 -5.9373 -2.8081 NFIT 2 1.074928 4.608755 -4.6165 -0.2785 NFIT 3 0.552656 0.000000 -0.8730 0.8010 NFIT 4 0.810044 -4.078986 -0.0192 3.5784 NFIT 5 0.827933 -8.229531 1.3999 5.7208 FIG3b: categorical allocation gap h against reversal probability. HGRID : 0.00 0.10 0.20 0.30 0.50 0.70 0.90 CATANA : 0.000000 0.133333 0.252604 0.320602 0.385417 0.416029 0.433771 CATPOP : 0.000000 0.134925 0.252830 0.320650 0.385755 0.414475 0.434855 CATOBS : 0.000000 0.052785 0.134345 0.194405 0.247950 0.265570 0.259335 CATOBSSE : 0.000000 0.000500 0.000763 0.000885 0.000966 0.000988 0.000980 kidney-stone gap h = 0.523, analytic P = 0.390018 FIG2: the lever. Pooled slope as a function of the between-group weight 1-w, for b_w = +1 and three between-group slopes. LEVERX : 0.00 0.05 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.85 1.00 LEVER -1.0 : 1.0000 0.9000 0.8000 0.6000 0.4000 0.2000 0.0000 -0.2000 -0.4000 -0.7000 -1.0000 LEVER -3.0 : 1.0000 0.8000 0.6000 0.2000 -0.2000 -0.6000 -1.0000 -1.4000 -1.8000 -2.4000 -3.0000 LEVER -8.0 : 1.0000 0.5500 0.1000 -0.8000 -1.7000 -2.6000 -3.5000 -4.4000 -5.3000 -6.6500 -8.0000 flagship dataset sits at 1-w = 0.87500, b_b = -1.2857 FIG3: reversal frequency against confounder strength. CGRID : 0.00 0.25 0.50 1.00 2.00 4.00 8.00 CURVE random 1:1 : 0.00000 0.00000 0.00340 0.09695 0.28850 0.42690 0.48360 CURVE random 3:1 : 0.00000 0.00000 0.00455 0.08700 0.28105 0.42240 0.48410 CURVE random 10:1 : 0.00000 0.00000 0.00335 0.06795 0.24370 0.40310 0.46995 CURVE random 30:1 : 0.00000 0.00000 0.00175 0.04855 0.19285 0.34935 0.43295 CURVE coupled 1:1 : 0.00000 0.00000 0.00350 0.09395 0.29405 0.42845 0.47665 CURVE coupled 3:1 : 0.00000 0.00000 0.00325 0.08150 0.28045 0.41800 0.47305 CURVE coupled 10:1 : 0.00000 0.00000 0.00120 0.04950 0.22690 0.38235 0.47235 CURVE coupled 30:1 : 0.00000 0.00000 0.00040 0.02425 0.14905 0.31475 0.41925 CURVE gated random 1:1 : 0.00000 0.00000 0.00100 0.06040 0.24090 0.39560 0.46915 CURVE se random 1:1 : 0.00000 0.00000 0.00041 0.00209 0.00320 0.00350 0.00353 FIG5: continuous case. CONTX : 0.25 0.50 1.00 1.50 2.00 3.00 5.00 8.00 CONTANA : 0.000000 0.005540 0.097321 0.181656 0.232655 0.284345 0.321256 0.337879 CONTPOP : 0.000000 0.005647 0.097863 0.181820 0.233480 0.285490 0.321078 0.337750 CONTSAMP : 0.000000 0.008467 0.097800 0.179467 0.232867 0.285967 0.323200 0.339500 CONTSE : 0.000000 0.000529 0.001715 0.002216 0.002440 0.002609 0.002700 0.002734 CONTFINEX : 0.10 0.30 0.50 0.70 0.90 1.10 1.30 1.50 1.70 1.90 2.10 2.30 2.50 2.70 2.90 3.10 3.30 3.50 3.70 3.90 4.10 4.30 4.50 4.70 4.90 5.10 5.30 5.50 5.70 5.90 6.10 6.30 6.50 6.70 6.90 7.10 7.30 7.50 7.70 7.90 CONTFINEY : 0.000000 0.000012 0.005540 0.034186 0.076167 0.117315 0.152638 0.181656 0.205229 0.224412 0.240134 0.253133 0.263984 0.273124 0.280892 0.287548 0.293293 0.298287 0.302657 0.306502 0.309904 0.312930 0.315633 0.318059 0.320244 0.322220 0.324013 0.325646 0.327137 0.328502 0.329756 0.330910 0.331974 0.332959 0.333872 0.334720 0.335509 0.336244 0.336931 0.337573 CATFINEX : 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66 0.68 0.70 0.72 0.74 0.76 0.78 0.80 0.82 0.84 0.86 0.88 0.90 CATFINEY : 0.000000 0.026667 0.053333 0.080000 0.106667 0.133333 0.160000 0.186437 0.211100 0.233153 0.252604 0.269714 0.284795 0.298139 0.310002 0.320602 0.330119 0.338704 0.346483 0.353561 0.360026 0.365953 0.371406 0.376437 0.381094 0.385417 0.389438 0.393190 0.396697 0.399983 0.403067 0.405968 0.408702 0.411281 0.413720 0.416029 0.418218 0.420296 0.422272 0.424152 0.425944 0.427653 0.429286 0.430847 0.432340 0.433771 ============================================================================== 12. SUMMARY TABLE: club value beside the accepted value ============================================================================== quantity club value reference -------------------------------------------------------------------------- OLS slope, 5-point hand dataset 0.600000000000 3/5 exact OLS intercept, same dataset 2.200000000000 11/5 exact t_{0.975, df=10} 2.2281 2.228 (table) Weighted-average identity residual 1.11e-16 0 exact Constructed within slope (flagship) 1.000000000000 +1 by design Constructed pooled slope (flagship) -1.000000000000 -1 by design Worst deviation, 8 constructions 4.44e-16 0 exact Control, c=0 balanced 0 in 40000 0 exact Control, c=0 imbalance 30:1 0 in 40000 0 exact Reversal rate, c=1 balanced 0.09695 no published value Gated reversal rate, c=1 balanced 0.06040 no published value Reversal rate, c=8 balanced 0.48360 ceiling 0.5 2x2x2 table, directional 0.008371 1/60 = 0.016667 2x2x2 table, either direction 0.016734 2/60 = 0.033333 Categorical, allocation gap h=0 0.000000 0 exact Categorical, allocation gap h=0.5 0.385755 0.385417 closed form Continuous s=1, population MC 0.097863 0.097321 integral Continuous s=1, n=200 sample MC 0.097800 no closed form Berkeley pooled gap (women-men) -14.2 pp published table Berkeley standardised gap +4.2 pp published table Kidney stones, pooled gap -4.6 pp published table Kidney stones, small-stone gap +6.4 pp published table ------------------------------------------------------------------------------ Wall clock: 83.9 s Seed: 20260214. Rerunning this file reproduces every number above.