Skip to content

ResearchMethodsVolume 2, Issue 2, Winter 2026

Building a Dataset Where Every Group Trends Up and the Whole Trends Down

We wrote a generator that produces Simpson's paradox to order, then ran the experiment backwards to measure how often a sign reversal happens by accident. With a confounder as strong as the within-group spread, 9.70% of simulated datasets reversed and 6.04% reversed with significance on both sides.

Written jointly by the Science Journaling Club

Study
Computational study
Review
Peer-edited by the club review board
Length
12 sectionsabout 30 min

Abstract

Simpson's paradox is usually shown once, with one famous table, and then left alone as a curiosity. We treated it as a piece of engineering instead. The club wrote a generator that takes a target within-group slope and a target pooled slope of the opposite sign and returns a dataset hitting both, by solving for the required between-group slope rather than searching for one. Across eight specifications the worst deviation between the fitted slopes and the requested ones was 4.44×10−16. Nothing was measured in the world; the computation is the experiment. Then we ran the machine backwards and asked how often the paradox arrives uninvited. With four groups of sixty points, a within-group slope of +1+1, and a group-level confounder as strong as the within-group spread in xx, 9.70% ± 0.21% of 20,000 simulated datasets showed every group sloping one way and the pooled fit sloping the other, and 6.04% ± 0.17% did so with every slope involved clearing a two-sided tt test at α=0.05\alpha = 0.05. Turning the confounder off gave 0 reversals in 40,000 datasets, which it must, because the pooled slope is then equal to the within-group slope as an algebraic identity. Group-size imbalance turned out to work the opposite way from the folklore, cutting the rate from 9.70% to 4.86% at a 30:1 ratio, and coupling size to the confounder cut it further to 2.43%, a difference of 13.0 standard errors. In the continuous case with a single lurking variable the Monte Carlo matched a closed-form integral to within 1.60 standard errors at every strength. As an external check, our 3,000,000 random 2×2×2 tables reversed at 0.016734 ± 0.000074 against Pavlides and Perlman's exact 1/60 [10], a gap of 0.90 standard errors. Seed 20260214.

The Trick, Performed Once

Here is a dataset. Five groups, fifteen points each, seventy-five numbers in all. Inside every single group the relationship between xx and yy is a perfect straight line climbing at exactly one unit of yy per unit of xx. Not approximately. Exactly. The correlation inside each group is +1.000000+1.000000.

Fit one line to all seventy-five points and the slope is −1.000000-1.000000. The correlation is −0.797724-0.797724.

Every part of the data goes up. The data goes down.

+1.000000slope inside each of 5 groups
−1.000000slope through all 75 points
4.4e−16worst error over 8 such builds
0 in 40,000reversals with the confounder off

We did not hunt for this dataset. We did not generate ten thousand random tables and keep the one that misbehaved. We wrote down the two numbers we wanted, +1+1 and −1-1, and the program solved one linear equation and handed back the points. It would have handed back a dataset for any other pair we asked for, including pairs far more violent than this one, and section three shows the price list.

-6-3036xexact construction-6-3036xwith noise, sigma 0.8-8-4048y
  • within-group fit, slope +1
  • pooled fit, slope -1
Figure 1. The constructed dataset. Left, the exact build: five groups laid out along xx, each with an internal slope of exactly +1+1, and the single line fitted to all of them falling at exactly −1-1. Right, the same construction with independent Gaussian noise of σ=0.8\sigma = 0.8 added inside each group, which is what the thing looks like when it is pretending to be data. The noisy version fits within-group slopes of 0.92, 1.07, 0.55, 0.81 and 0.83, and a pooled slope of −1.009-1.009; the largest discrepancy from the design across all six fits is 1.05 standard errors.

The right-hand panel matters more than the left. A perfectly collinear group is obviously artificial and nobody would be fooled. Sprinkle noise on it and the five groups become five modest positive findings, none of them remarkable, each the sort of thing a person might report in a paragraph. The pooled line is still an emphatic negative. Everything about that panel is publishable.

How It Is Done

The mechanism is a weighted average, and once you have seen it you cannot unsee it.

Take grouped data. For group gg holding ngn_g points, write xˉg\bar{x}_g for its mean in xx and xˉ\bar{x} for the grand mean. Define two sums of squares:

Tw=∑g∑i(xgi−xˉg)2,Tb=∑gng(xˉg−xˉ)2.T_w = \sum_g \sum_i (x_{gi} - \bar{x}_g)^2, \qquad T_b = \sum_g n_g (\bar{x}_g - \bar{x})^2 .

TwT_w is how much xx varies inside groups. TbT_b is how much the groups are displaced from each other along xx. Their sum is the total variation in xx, because those two quantities partition it exactly. Now let w=Tw/(Tw+Tb)w = T_w/(T_w + T_b), let bwb_w be the common within-group slope, and let bbb_b be the least-squares slope of the GG group means (xˉg,yˉg)(\bar{x}_g, \bar{y}_g) weighted by ngn_g. Then the slope ordinary least squares returns for the pooled data is

bpooled  =  w bw  +  (1−w) bb.b_{\text{pooled}} \;=\; w\, b_w \;+\; (1 - w)\, b_b .

That is the whole apparatus. The pooled slope is an average of the within-group slope and the between-group slope, and the mixing weight is decided by nothing except geometry: how spread out the groups are in xx compared with how spread out the points are inside them. Simpson's paradox is what a weighted average does when the two things being averaged have opposite signs and the weight is large enough. It is not a paradox at all. It is a mean.

Which gives the constructor immediately. Choose bwb_w. Choose the pooled slope bpb_p you want. Lay out the group means, which fixes TbT_b, and choose the within-group spread, which fixes TwT_w. Then solve:

bb  =  bp−w bw1−w.b_b \;=\; \frac{b_p - w\, b_w}{1 - w} .

Place the group means on a line of that slope and you are finished. There is no search and no rejection sampling. Any pair of signs is reachable, for any GG, because the equation has a solution whenever w<1w < 1, and w<1w < 1 whenever the groups are displaced from one another at all.

We checked the identity to machine precision on a deliberately awkward dataset: five groups of sizes 37, 12, 61, 25 and 9, with different within-group spreads and different within-group slopes, none of which the identity requires. The fitted pooled slope was −0.325516544-0.325516544; the identity predicted −0.325516544-0.325516544; the difference was 1.11×10−161.11\times10^{-16}, a relative error of 3.4×10−163.4\times10^{-16}.

b_b = -1b_b = -3b_b = -8Figure 1 sits here0.00.20.40.60.81.0-8-6-4-201 - w, the share of the spread in x lying between groupspooled slope, within-group slope held at +1slope
Figure 2. The lever. With the within-group slope held at +1+1, the pooled slope is a straight line in 1−w1 - w, the share of the spread in xx that lies between groups rather than inside them. The three lines are three different between-group slopes. A steeper between-group line crosses zero sooner. The dataset in Figure 1 sits at 1−w=0.8751 - w = 0.875 with a between-group slope of −1.2857-1.2857, which is the single point marked. Nothing here was simulated; every value on this plot is one multiplication.

The Arithmetic

Work the flagship dataset by hand. The five groups sit at x=−4,−2,0,2,4x = -4, -2, 0, 2, 4 with fifteen points each, so xˉ=0\bar{x} = 0 and the between-group sum of squares is 15×(16+4+0+4+16)=60015 \times (16 + 4 + 0 + 4 + 16) = 600. The within-group offsets are fifteen evenly spaced values scaled so that each group's sum of squares comes to 85.714286/5=17.14285785.714286/5 = 17.142857, giving Tw=85.714286T_w = 85.714286. So

w=85.71428685.714286+600=0.125=18.w = \frac{85.714286}{85.714286 + 600} = 0.125 = \tfrac18 .

One eighth. The pooled fit gives seven eighths of its attention to where the groups sit and one eighth to what happens inside them. Now solve for the between-group slope that produces bp=−1b_p = -1:

bb=−1−18(1)1−18=−9/87/8=−97=−1.285714.b_b = \frac{-1 - \tfrac18 (1)}{1 - \tfrac18} = \frac{-9/8}{7/8} = -\frac97 = -1.285714 .

Check it forward. 18(1)+78(−9/7)=0.125−1.125=−1\tfrac18 (1) + \tfrac78 (-9/7) = 0.125 - 1.125 = -1. The group means then have to sit at y=−97xy = -\tfrac97 x, which for x=−4x = -4 gives y=36/7=5.142857y = 36/7 = 5.142857, and the program prints 5.14295.1429.

Notice the lever arm. It is 1/(1−w)=8/7=1.1431/(1-w) = 8/7 = 1.143. To move the pooled slope from +1+1 down to −1-1, a distance of 2, the between-group line has to move by 2×8/7=2.2862 \times 8/7 = 2.286, from +1+1 to −1.2857-1.2857. That is barely any effort at all. The layout in Figure 1 is not an extreme one. Groups two units apart with unit spread inside them is an ordinary-looking picture, and it already gives the between-group line seven eighths of the vote.

Push harder and the price stays cheap. Asking for a pooled slope of −5-5 needs a between-group slope of −5.857-5.857. Asking for −20-20 needs −23.0-23.0. The cost of an absurd answer is an absurd-looking layout, and the cost of a plausible answer is a layout nobody would look at twice.

Before any of that we checked the regression routine itself, because a study about arithmetic going wrong should not be running on arithmetic nobody checked. Five points: x=1,2,3,4,5x = 1,2,3,4,5 and y=2,4,5,4,5y = 2,4,5,4,5. By hand, xˉ=3\bar{x} = 3, yˉ=4\bar{y} = 4, Sxx=10S_{xx} = 10, Sxy=6S_{xy} = 6, so the slope is 6/10=3/56/10 = 3/5 and the intercept is 4−3×3/5=11/54 - 3\times 3/5 = 11/5. The club's routine returns 0.6000000000000.600000000000 and 2.2000000000002.200000000000, agreeing with the exact rational answer to the last bit. Rational arithmetic also gives the residual sum of squares as 12/512/5, from which the standard error of the slope is 0.28284271250.2828427125; the routine returns the same to 5.6×10−175.6\times10^{-17}. The Student tt quantiles, which the significance test later needs, come from our own continued-fraction incomplete beta function and match published tables at every degrees of freedom we use, the worst gap being 4.2×10−44.2\times10^{-4} at df=5\mathrm{df} = 5, where the published table is rounded to three decimals anyway.

Working Notes From the Club Table

Meeting 1
Somebody asked whether you could make the paradox happen on purpose. Twenty minutes of trying to do it by fiddling with numbers on the board. Got close twice. Gave up and did the algebra instead, which took four lines.

Meeting 1, later
Wrote construct(). Worked first time, which never happens. Then it did not work when noise was added, because the noise was shifting the group means and therefore moving TbT_b out from under the solved bbb_b. Fix: subtract the noise's own mean inside each group before adding it. Three characters.

Meeting 2
Argument about what counts as a reversal. Point estimates only, or significance on both sides? Settled it by reporting both. The gated number is smaller and more alarming.

Meeting 2
First frequency run had every reversal in the same direction, because the confounder's effect on yy was hard-coded negative. That is not a measurement of how often something happens by accident. Rewrote so the loading θ\theta is drawn fresh each dataset from N(0,1)N(0,1). The ceiling of the curve is now one half, by construction, and we say so.

Meeting 3
Imbalance result came back the wrong way round. Spent most of the meeting looking for the bug. There is no bug, and §7 is what we found instead. This is the result we would most like somebody to check.

Meeting 3
Categorical simulation sat 3.4 standard errors off its own closed form. Cause: we were clipping success probabilities at 1, which bit on about 3% of draws and quietly removed part of the treatment effect. Narrowed the draw range. Worst cell is now 2.1 SE. The closed form was right and we were wrong, which is the correct way round for that to happen.

Meeting 4
Seed fixed at 20260214 and not touched again. Every number in this article comes from one run of the file as it now stands, 80.9 seconds on a laptop.

The Control, Which Is the Dull Half

An argument of the form "this can happen" needs a companion of the form "and here is when it cannot". Otherwise the reader has no way to tell whether the simulation is measuring the phenomenon or measuring a bug.

So the frequency experiment is built as a two-by-two. The generator has exactly two ingredients that could plausibly cause a reversal: a confounder strength cc, which sets how far apart the groups sit along xx, and a group-size imbalance ratio. We switched each off and on. Four groups, 240 points, a true within-group slope of +1+1, residual standard deviation 1, and 40,000 datasets in each cell.

CellConfounderSizesReversalsRate
Controlc = 060, 60, 60, 6000.000000
Imbalance onlyc = 05, 17, 53, 16500.000000
Confounder onlyc = 160, 60, 60, 6037950.094875
Bothc = 15, 17, 53, 16517630.044075

Zero and zero. Not "approximately zero", not "a handful", zero, in 80,000 datasets between them. The reason is that the control is not a probabilistic statement. Set c=0c = 0 and all the group means in xx coincide, so Tb=0T_b = 0, so w=1w = 1, so bpooled=bwb_{\text{pooled}} = b_w identically. The pooled fit cannot disagree with the groups because it is the same number. Pooling has no power to reverse anything unless it is handed a between-group displacement to work with.

By the rule of three, zero events in 40,000 trials puts a one-sided 95% upper bound of 7.5×10−57.5\times10^{-5} on the true rate, which is as close to confirming an exact zero as counting can get. The categorical experiment in §9 has the same structure and the same control: set the allocation gap to zero and the reversal count across 200,000 simulated trials is also exactly zero.

How Often It Happens on Its Own

The reverse experiment. Instead of building the paradox, generate datasets with a confounder in them and count how often one turns up uninvited.

The generator: four groups, 240 points, a common within-group slope of +1+1, residual σ=1\sigma = 1, within-group spread in xx of 1. A group-level lurking variable zgz_g is drawn standard normal. It shifts group gg's mean in xx by czgc z_g and its level in yy by cθzgc \theta z_g, where θ∼N(0,1)\theta \sim N(0,1) is drawn fresh for each dataset. That last choice is the one carrying the most weight in this article and it deserves stating plainly: it means the confounder is exactly as likely to push with the within-group effect as against it. "By accident" has to mean something, and this is what we decided it means. A consequence is that the reversal curve has a ceiling at one half, and that ceiling is a design decision rather than a finding. Where the curve reaches it is the finding.

ccReversal rateSEGated rateSEMean ww
00.000000.000000.000000.000000.9875
0.250.000000.000000.000000.000000.9442
0.50.003400.000410.001000.000220.8442
1.00.096950.002090.060400.001680.6221
2.00.288500.003200.240900.003020.3415
4.00.426900.003500.395600.003460.1382
8.00.483600.003530.469150.003530.0457

Four balanced groups, 20,000 datasets per row. Read the c=1c = 1 row twice. A confounder whose between-group spread merely equals the within-group spread, which is not a dramatic confounder by any standard, produces a full sign reversal in 9.70% of datasets. One in ten. And 6.04% of datasets do it while every within-group slope and the pooled slope all clear a two-sided tt test at α=0.05\alpha = 0.05, which is to say that in one dataset in sixteen a careful analyst would find a significant positive effect in every group and a significant negative effect overall, with confidence intervals that do not touch zero on either side.

The gap between the two columns closes as cc grows. At c=8c = 8 the gated rate is 0.46915 against a plain rate of 0.48360, so 97% of the reversals present are also statistically convincing in both directions at once. A large confounder does not produce marginal, arguable reversals. It produces clean ones.

ceiling 0.501:13:110:130:1gated00.250.512480.00.10.20.30.40.5confounder strength c, in within-group SDs of xfraction of datasets with a sign reversalrateratios are largest group : smallest group, total held at 240error bars are 2 SE from 20,000 trials per point
Figure 3. Reversal frequency against confounder strength, for four group-size ratios, 20,000 datasets per point. The dashed curve is the significance-gated rate for balanced groups, which requires all four within-group slopes and the pooled slope to clear a two-sided tt test at α=0.05\alpha = 0.05. Error bars on the balanced curve are two standard errors. The ceiling at 0.50 is the modelling choice that the confounder's loading on yy is symmetric about zero, not a measured limit.

The Thing About Unequal Groups That Turned Out to Be Backwards

Everyone says unequal group sizes make Simpson's paradox more likely. We say it too, or we did. The Berkeley story is always told that way: the big departments were the hard ones, women applied to them, the marginal got dragged down. Imbalance sounds like the villain.

Our grid says otherwise, at least for the regression form. Holding c=1c = 1 and sliding the largest-to-smallest group ratio from 1:1 to 30:1, the reversal rate falls from 0.09695 to 0.04855. Every step down is far outside its standard error. At c=2c = 2 it falls from 0.28850 to 0.19285.

The reason is in the same identity as everything else. For two groups of sizes n1n_1 and n2n_2 separated in xx by a distance DD,

Tb=n1n2n1+n2D2,T_b = \frac{n_1 n_2}{n_1 + n_2} D^2 ,

which for a fixed total is largest when n1=n2n_1 = n_2. Between-group leverage is a product, and products of a fixed sum are maximised at balance. A dataset dominated by one enormous group has almost no between-group leverage to give away, so the pooled line is pinned to the big group's own slope and comes out agreeing with it. Watch the mean-ww column in the raw output climb from 0.6221 to 0.7253 as the ratio goes from 1:1 to 30:1. The pooled fit is listening to the groups more, not less.

We then tried to rescue the folklore by coupling size to the confounder, putting the biggest group at the most extreme value of the lurking variable, which is the Berkeley shape. It made things worse for the folklore, not better. At c=1c = 1 and 30:1 the coupled rate is 0.02425 against 0.04855 for random assignment, a difference of −0.0243±0.0019-0.0243 \pm 0.0019, which is 13.0 standard errors. The mechanism is the grand mean: park 165 of your 240 points at one extreme and xˉ\bar{x} follows them there, which collapses that group's own deviation (xˉg−xˉ)(\bar{x}_g - \bar{x}) toward zero and takes the between-group sum of squares with it.

So in the regression form, unequal group sizes are a protection. This is the one result in the study that we would most like somebody to try to break, and the code is a single file.

One Lurking Variable and No Groups at All

Groups are a convenience. The underlying problem has nothing to do with them, and the continuous version shows it with less machinery. Let ZZ be standard normal and write

X=aZ+u,Y=βX+bZ+v,X = aZ + u, \qquad Y = \beta X + bZ + v,

with uu and vv independent standard normal errors. The effect of XX on YY holding ZZ fixed is β\beta, and we set it to +0.5+0.5, positive by construction. The marginal slope, the one you get if nobody records ZZ, is

bmarg=Cov(X,Y)Var(X)=β+aba2+1,b_{\text{marg}} = \frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(X)} = \beta + \frac{ab}{a^2 + 1},

so the sign flips exactly when ab/(a2+1)<−βab/(a^2+1) < -\beta. Draw the loadings aa and bb independently from N(0,s2)N(0, s^2), with ss the lurking strength, and that probability becomes a one-dimensional integral with no simulation in it at all:

P(s)=2∫0∞ϕs(a)[1−Φ ⁣(β(a2+1)as)] da.P(s) = 2\int_0^\infty \phi_s(a)\left[1 - \Phi\!\left(\frac{\beta(a^2+1)}{as}\right)\right]\,da .

Which is useful for a reason beyond elegance. It gives the Monte Carlo something exact to be judged against. At s=1s = 1, 400,000 draws gave 0.097863 ±\pm 0.000470 against an integral value of 0.097321, a gap of 1.15 standard errors. Across all eight strengths the worst disagreement was 1.60 standard errors. The one cell with no events at all, s=0.25s = 0.25, is exactly what should happen: the integral predicts 0.0751 events in 400,000 draws, and the Poisson probability of observing none of them is 0.928.

The second run is the one worth staring at. Instead of evaluating the population slope we drew 200 actual observations per dataset and judged the reversal from the fitted marginal slope against the fitted slope adjusted for ZZ. Where the population rate is already large, the two agree inside a standard error, because a reversal that is going to happen happens whether or not the sample is noisy. Where the population rate is small they separate. At s=0.5s = 0.5 the population rate is 0.00565 and the sample rate is 0.00847, an excess of 5.2 standard errors. Finite samples manufacture reversals in precisely the region where confounding alone would produce none.

Rate Tables, and an Answer Somebody Already Knew

A simulation that has never been made to agree with an independently derived number is a rumour. So before reporting anything we did not have a check for, we went looking for something that was already known exactly.

Pavlides and Perlman [10], using a proof due to Hadjicostas, give the probability that a random 2×2×2 table shows Simpson's paradox when the eight cell probabilities are drawn uniformly from the simplex. The answer is exactly 1/601/60. It is a clean target, with no confidence interval attached to it and nothing to argue about.

We drew 3,000,000 such tables from a Dirichlet(1,…,1)(1,\ldots,1) and counted. Counting reversals in one nominated direction gives 25,113 events, a rate of 0.0083710 ±\pm 0.0000526. Counting them in either direction gives 50,201, a rate of 0.0167337 ±\pm 0.0000741. The published 1/601/60 is 0.0166667. Our either-direction rate sits 0.90 standard errors from it. Our one-directional rate sits 0.72 standard errors from 1/1201/120, which is exactly half, as the symmetry between the two treatment labels demands. So the published convention counts both directions, and our machinery reproduces it.

We print both numbers in the raw output, including the two comparisons that are meant to fail, because reporting only the one that agreed would have been the easy dishonest move and there would have been no way for a reader to catch it.

1/120 = 0.00833310^310^410^510^63e60.0040.0060.0080.0100.0120.0140.016tables drawnrunning estimate, band is +/- 1 SE1/120 is half of the published 1/60, which countsreversals in either direction
Figure 4. Convergence. The running Monte Carlo estimate of the one-directional 2×2×2 reversal rate as tables accumulate, with a band of one standard error, against its exact value of 1/1201/120. The estimate wanders by a factor of two before 10,000 tables and settles by about 100,000. At 3,000,000 it sits 0.72 standard errors from truth. 1/1201/120 is half of the published 1/601/60 [10], which counts reversals in either direction.

Which brings the second machine into view. The regression paradox and the rate-table paradox get taught as one thing. They are not one thing, and §7 is the proof: a risk factor for one of them is a protective factor for the other.

In a rate table there is no leverage weighting. An arm's pooled success rate is just its stratum-specific rates averaged with that arm's own allocation weights. So a reversal requires the two arms to be allocated differently across the strata, and nothing else will do it. Stratum size as such is irrelevant.

We simulated it. Two strata, two arms, 350 patients per arm to match the kidney-stone series. Arm B's stratum success probabilities q1,q2q_1, q_2 are drawn uniform on [0.05,0.85][0.05, 0.85], arm A gets an advantage δ\delta drawn uniform on [0,0.10][0, 0.10] applied in both strata so that A is strictly better everywhere, and the arms are allocated at fA=0.5+sh/2f_A = 0.5 + sh/2 and fB=0.5−sh/2f_B = 0.5 - sh/2 with s=±1s = \pm1 at random. Then hh is the allocation gap, and the population reversal condition is exactly sh(q1−q2)+δ<0sh(q_1 - q_2) + \delta < 0, which integrates to

P(h)=h2δmax⁡⋅L3[1−(L−UL)3],L=0.8,U=min⁡(δmax⁡/h, L).P(h) = \frac{h}{2\delta_{\max}}\cdot\frac{L}{3}\left[1 - \left(\frac{L-U}{L}\right)^{3}\right], \qquad L = 0.8,\quad U = \min(\delta_{\max}/h,\, L).

Another closed form, and another chance to be caught out. At h=0.5h = 0.5 the simulation gives 0.385755 against a predicted 0.385417. The worst cell over the whole sweep is 2.08 standard errors. At h=0h = 0 the count is exactly zero in 200,000 trials, which is the categorical control.

Now put the real numbers in. Charig and colleagues [7] reported 350 kidney-stone patients treated by open surgery and 350 by percutaneous nephrolithotomy. Open surgery won among small stones, 93.1% against 86.7%, and won among large stones, 73.0% against 68.8%. Pooled, it lost, 78.0% against 82.6%. The whole of that 4.6-point reversal comes from allocation: 75.1% of the open-surgery patients had large stones against 22.9% of the keyhole patients, an allocation gap of h=0.523h = 0.523. Read that gap off our curve and the probability of a reversal, for a treatment whose real advantage lies somewhere in nought to ten points, is 0.39. Julious and Mullee put those same figures in front of clinicians in 1994 under the heading of confounding rather than paradox [8], which is the better name for it.

The Berkeley figures behave the same way [5]. Across the six largest departments, women were admitted at the higher rate in four of the six. Pooled, the gap runs 14.2 percentage points against them, 30.3% against 44.5%. Standardise both sexes to the combined applicant pool and the gap becomes 4.2 points in their favour. Same numbers, two readings, opposite signs. This is the example everybody uses, and it is still the best one [4][6][11].

024680.00.10.20.30.4lurking strength scontinuous casekidney stones0.00.20.40.60.80.00.10.20.30.4allocation gap hcategorical caserate
  • closed form
  • population value
  • what a finite sample shows

left: 30,000 datasets of 200 observations, bars 2 SE. right: 350 per arm

Figure 5. Two versions of the same question. Left, the continuous case: reversal probability against lurking-variable strength, with the closed-form integral as the line, the population Monte Carlo as squares and 30,000 datasets of 200 observations each as circles with two-standard-error bars. Right, the categorical case: reversal probability against the allocation gap hh, with the closed form as the line, the population Monte Carlo as squares and what a 350-per-arm trial would actually observe as circles. The vertical mark is the kidney-stone series of Charig et al. [7] at h=0.523h = 0.523.

The Strongest Objection

We put the study to the room and asked for the hardest thing anybody could say against it. This was it, and it is good.

The objection

Your 9.70% is manufactured. You chose θ∼N(0,1)\theta \sim N(0,1), which says the confounder's pull on yy is on the same scale as the real effect you planted. Of course a coin-flip confounder of comparable size flips signs about a tenth of the time. You have measured your own prior and called it a frequency. In real research the confounder is usually weaker than the effect, and researchers usually know their field well enough to have measured the obvious confounders already. The number that matters is not how often this happens in your generator. It is how often it happens in journals, and your simulation cannot see a journal.

All of that is correct. The 9.70% is a property of our generator and carries no authority outside it. We can say how much it moves, because we measured that too, and the sensitivities are large. Halve the confounder's loading on yy to τ=0.5\tau = 0.5 and the rate drops to 0.02890. Double it to τ=2\tau = 2 and it rises to 0.19270. Both are more than 27 standard errors from the base case. A reader who believes confounders in their field are half our size should read 2.9%, not 9.7%, and we would not argue with them.

The ratio that actually governs the answer is not the confounder's absolute size but its size relative to the effect you are looking for. Hold the confounder fixed and weaken the real within-group slope from +1+1 to +0.25+0.25 and the rate climbs to 0.25680. Strengthen it to +4+4 and the rate collapses to 0.00405. That is a factor of sixty-three across a factor of sixteen in effect size, and it says something uncomfortable about where the risk concentrates. Fields hunting small effects are not slightly more exposed to this. They are the only ones seriously exposed to it.

Residual noise barely matters at all. Tightening the groups to σ=0.25\sigma = 0.25 moves the rate by 1.6 standard errors and loosening them to σ=4\sigma = 4 moves it by 2.4, against shifts of more than forty standard errors for the same relative change in effect size. So the phenomenon is not about how noisy the data are. It is about where the data sit.

The second half of the objection we concede without qualification. We cannot estimate how often this occurs in published work, because that quantity is a property of scientific practice rather than of probability, and no amount of simulation reaches it. What our numbers give is a floor under a conditional statement: if a confounder of this strength is present and unrecorded, then the reversal rate is this. Whether such confounders are present is a question for the people who know the field, which in our case is nobody.

One more thing on the objection's side, which nobody in the room raised and which we should have. Our model gives every group the same within-group slope. Real groups do not have the same slope, and unequal slopes make reversals easier rather than harder, because they give the pooled fit one more way to disagree with its parts. Our frequencies should be read as a floor in that direction too.

What the Numbers Do Not Decide

Sample size is the part that unsettles us most, and it is the opposite of reassuring. Holding the confounder at c=1c = 1 and growing the dataset from 40 points to 2,400, the plain reversal rate barely moves: 0.10340, 0.09750, 0.09245, 0.09215, 0.08980. It is set by the population geometry and sampling error has almost nothing to do with it. But the significance-gated rate climbs steadily, 0.00435, 0.03490, 0.05850, 0.07365, 0.07805. More data does not dissolve the contradiction. It sharpens it, and converts a discrepancy you could have written off as noise into two confident findings pointing opposite ways.

The number of groups pushes the same direction. At G=2G = 2 the rate is 0.05620 and at G=12G = 12 it is 0.13340, because more groups means less of the spread in xx sits inside them. Finer stratification of the same data raises the exposure.

Where a different modelling choice would have changed the answer, plainly. Drawing θ\theta from a distribution with a mean rather than centred at zero would have moved the whole curve up or down and removed the one-half ceiling; we chose symmetry because "by accident" seemed to require it, and a reader who thinks confounders in their field have a characteristic direction should expect a different number. Letting within-group slopes vary would have raised every frequency. Using heavy-tailed errors instead of Gaussian ones would have raised them again. Treating group membership as unobserved, which is the case that actually frightens people, would not have changed any frequency at all, because the arithmetic does not care whether you recorded the confounder. It only changes whether you can see what happened.

And the limit that no version of this study can pass. Nothing here tells you which answer is correct. The paradox is a statement about arithmetic and contains no information about causation whatsoever [16][12]. Whether to report the within-group number or the pooled number depends on what generated the data, which is not in the data. Stratify by something that happens after the exposure and the stratified answer is the wrong one. Stratify by a genuine confounder and the pooled answer is the wrong one. The same two columns of figures support both readings and the figures will never tell you which. Simpson said as much in 1951, in a paper that spends most of its four pages on exactly this point [1], and Pearson had already run into the same wall in 1899 while trying to understand correlations between mixed populations of racehorses [3]. It keeps coming back, in meta-analysis [13] as readily as in psychology [11], and it came back in 2020 when case fatality rates for one disease were higher in every age band in one country and lower overall [14].

What we can offer is a habit rather than a rule. When you read an aggregate number, ask what ww was. Ask how much of the variation in the thing being averaged sits between the groups rather than inside them, because that single fraction decides how much of the headline belongs to the effect and how much belongs to where the groups happened to be standing. If ww is small, the headline is mostly a map of the group layout. Our flagship dataset had w=0.125w = 0.125 and looked entirely ordinary.

Reproducing this

One file, one command. Python 3.12 and numpy, nothing else.

python simpsons-paradox.py > simpsons-paradox-output.txt

Expect about 80 seconds. Ours took 80.9 s on numpy 2.4.2 and Python 3.12.3. The master seed is 20260214, hard-coded at the top of the file, and each experiment draws an independent stream from it through numpy's SeedSequence.spawn, so the output is deterministic to the last digit. Two of the comparison lines in section 7 of the output are labelled MISMATCH on purpose; they are the conventions being ruled out. Any other number differing from ours by more than the printed standard errors means something is wrong, and we would like to hear about it. The interactive model runs both experiments in your browser, and its default settings reproduce the headline figures above.

QuantityClub valueReference valueAgreement
Machinery
OLS slope, five-point hand dataset0.6000000000003/5 exactexact
OLS intercept, same dataset2.20000000000011/5 exactexact
SE of that slope0.2828427125(12/5)/3/10\sqrt{(12/5)/3/10}5.6e−17
t0.975, df=10t_{0.975,\,\mathrm{df}=10}2.22812.228 (table)1.4e−04
Weighted-average identity−0.325516544−0.3255165441.1e−16
Construction
Flagship within-group slope1.000000000000+1 by design2.2e−16
Flagship pooled slope−1.000000000000−1 by design2.2e−16
Worst of eight constructions4.44e−160 exactexact
Noisy build, largest ∣z∣|z| over six fits1.05order 1within SE
Controls
Regression, c=0c = 0, balanced0 in 40,0000 exactexact
Regression, c=0c = 0, 30:1 imbalance0 in 40,0000 exactexact
Categorical, allocation gap h=0h = 00 in 200,0000 exactexact
Frequencies
Regression, c=1c = 1, balanced0.09695 ± 0.00209none published—
Same, significance-gated0.06040 ± 0.00168none published—
Regression, c=1c = 1, 30:1 imbalance0.04855 ± 0.00152none published—
Same, size coupled to confounder0.02425 ± 0.00109none published—
Categorical, h=0.5h = 0.50.3857550.385417 closed form0.31 SE
Continuous, s=1s = 1, population0.0978630.097321 integral1.15 SE
Continuous, s=1s = 1, n=200n = 2000.097800no closed form—
Random 2×2×2, either direction0.016734 ± 0.0000741/60 = 0.016667 [10]0.90 SE
Random 2×2×2, one direction0.008371 ± 0.0000531/120 = 0.0083330.72 SE
Published tables, recomputed
Berkeley, pooled gap (women − men)−14.2 ppBickel et al. [5]—
Berkeley, standardised gap+4.2 ppBickel et al. [5]—
Kidney stones, small-stone gap+6.4 ppCharig et al. [7]—
Kidney stones, large-stone gap+4.3 ppCharig et al. [7]—
Kidney stones, pooled gap−4.6 ppCharig et al. [7]—

References

Reproduce these numbers

From the analysis folder, run:

python simpsons-paradox.py > simpsons-paradox-output.txt
Seed
20260214
Analysis code
simpsons-paradox.py63 kB