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

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 \(x\) and \(y\) is a perfect straight line climbing at exactly one unit of \(y\) per unit of \(x\). Not approximately. Exactly. The correlation inside each group is \(+1.000000\).

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

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

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\) and \(-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.

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 \(g\) holding \(n_g\) points, write \(\bar{x}_g\) for its mean in \(x\) and \(\bar{x}\) for the grand mean. Define two sums of squares:

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

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 \(x\) 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 \(b_w\). Choose the pooled slope \(b_p\) you want. Lay out the group means, which fixes \(T_b\), and choose the within-group spread, which fixes \(T_w\). Then solve:

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 \(G\), because the equation has a solution whenever \(w < 1\), and \(w < 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\); the identity predicted \(-0.325516544\); the difference was \(1.11\times10^{-16}\), a relative error of \(3.4\times10^{-16}\).

The Arithmetic

Work the flagship dataset by hand. The five groups sit at \(x = -4, -2, 0, 2, 4\) with fifteen points each, so \(\bar{x} = 0\) and the between-group sum of squares is \(15 \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.142857\), giving \(T_w = 85.714286\). So

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 \(b_p = -1\):

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

Notice the lever arm. It is \(1/(1-w) = 8/7 = 1.143\). To move the pooled slope from \(+1\) down to \(-1\), a distance of 2, the between-group line has to move by \(2 \times 8/7 = 2.286\), from \(+1\) to \(-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\) needs a between-group slope of \(-5.857\). Asking for \(-20\) needs \(-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,5\) and \(y = 2,4,5,4,5\). By hand, \(\bar{x} = 3\), \(\bar{y} = 4\), \(S_{xx} = 10\), \(S_{xy} = 6\), so the slope is \(6/10 = 3/5\) and the intercept is \(4 - 3\times 3/5 = 11/5\). The club's routine returns \(0.600000000000\) and \(2.200000000000\), agreeing with the exact rational answer to the last bit. Rational arithmetic also gives the residual sum of squares as \(12/5\), from which the standard error of the slope is \(0.2828427125\); the routine returns the same to \(5.6\times10^{-17}\). The Student \(t\) 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\times10^{-4}\) at \(\mathrm{df} = 5\), where the published table is rounded to three decimals anyway.

Working Notes From the Club Table

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 \(c\), which sets how far apart the groups sit along \(x\), and a group-size imbalance ratio. We switched each off and on. Four groups, 240 points, a true within-group slope of \(+1\), residual standard deviation 1, and 40,000 datasets in each cell.

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 = 0\) and all the group means in \(x\) coincide, so \(T_b = 0\), so \(w = 1\), so \(b_{\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\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\), residual \(\sigma = 1\), within-group spread in \(x\) of 1. A group-level lurking variable \(z_g\) is drawn standard normal. It shifts group \(g\)'s mean in \(x\) by \(c z_g\) and its level in \(y\) by \(c \theta z_g\), where \(\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.

Four balanced groups, 20,000 datasets per row. Read the \(c = 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 \(t\) test at \(\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 \(c\) grows. At \(c = 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.

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 = 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 = 2\) it falls from 0.28850 to 0.19285.

The reason is in the same identity as everything else. For two groups of sizes \(n_1\) and \(n_2\) separated in \(x\) by a distance \(D\),

which for a fixed total is largest when \(n_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-\(w\) 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 = 1\) and 30:1 the coupled rate is 0.02425 against 0.04855 for random assignment, a difference of \(-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 \(\bar{x}\) follows them there, which collapses that group's own deviation \((\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 \(Z\) be standard normal and write

with \(u\) and \(v\) independent standard normal errors. The effect of \(X\) on \(Y\) holding \(Z\) fixed is \(\beta\), and we set it to \(+0.5\), positive by construction. The marginal slope, the one you get if nobody records \(Z\), is

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

Which is useful for a reason beyond elegance. It gives the Monte Carlo something exact to be judged against. At \(s = 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.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 \(Z\). 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.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/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,\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/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/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.

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 \(q_1, q_2\) are drawn uniform on \([0.05, 0.85]\), arm A gets an advantage \(\delta\) drawn uniform on \([0, 0.10]\) applied in both strata so that A is strictly better everywhere, and the arms are allocated at \(f_A = 0.5 + sh/2\) and \(f_B = 0.5 - sh/2\) with \(s = \pm1\) at random. Then \(h\) is the allocation gap, and the population reversal condition is exactly \(sh(q_1 - q_2) + \delta < 0\), which integrates to

Another closed form, and another chance to be caught out. At \(h = 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 = 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.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].

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.

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 \(y\) to \(\tau = 0.5\) and the rate drops to 0.02890. Double it to \(\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\) to \(+0.25\) and the rate climbs to 0.25680. Strengthen it to \(+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 \(\sigma = 0.25\) moves the rate by 1.6 standard errors and loosening them to \(\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 = 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 = 2\) the rate is 0.05620 and at \(G = 12\) it is 0.13340, because more groups means less of the spread in \(x\) 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 \(w\) 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 \(w\) is small, the headline is mostly a map of the group layout. Our flagship dataset had \(w = 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.

References