Simpson's paradox
Two treatments' success rates group by group and pooled, compared exactly by cross-multiplication: one better in every group and worse all together, why the weights decide it, the average of rates against the pooled rate, and a comparison that weighs the groups the same.
Every screen below was recorded under CPython. When this page was built, the EML interpreter replayed each session from the same input and printed the same bytes.
About
Two treatments, A and B, and their success rates group by group and pooled. One can do better in every group and worse with the groups put together. The program shows why the weights decide it, sets the average of the rates against the pooled rate, and gives a comparison that weighs the groups the same for both. It starts with the kidney-stone figures of Charig and others (1986): A wins on small stones and on large stones, B wins overall.
main.eml- the menu, the table, the explanations, and editing groupsrates.eml- comparing rates, percentages, totals, the verdicts, the mean of rates and standardized rates
How each part works:
- Every rate is a fraction of two whole numbers. As in the corpus case
simpsons-paradox, two rates are compared by cross-multiplication: a/b beats c/d exactly when a d beats c b. No rounding can cause a reversal; percentages are only for showing. - A paradox is every group going one way and the pooled rates going the other. The pooled rate is each group's rate weighted by the share of that arm's cases the group holds, so it depends on who got which treatment: A was given mostly large stones, the harder cases.
- The average of the group rates, as in the corpus case
average-of-averages, counts each group once, whatever its size. The pooled rate counts each case once. They agree only when the groups are the same size. - Weighing the groups the same: each arm's rate as if it had treated each group in proportion to that group's share of all the cases. Then a treatment that is at least as good in every group, and better in one, always comes out ahead.
What is checked: menu choices 1 to 7; a group name of up to 16 characters, new when adding and existing when changing or removing; up to 6 groups and at least 1; successes and cases as two whole numbers, at least 1 case and no more successes than cases. An empty answer cancels.
Sessions: sessions/basic.in gives: - the kidney-stone table: A 93.1% and 73.0% in the groups, 78.0% pooled; B 86.7% and 68.8%, 82.6% pooled; - why: A's cases are 75.1% large stones, B's 22.9%; weighing the groups the same gives A 83.3% against B 77.9%; - the average of the rates against the pooled rates; - a changed group, an added one and its removal, the paradox holding throughout.
sessions/bad-input.in gives: - menu choices 0 and x; - a group that does not exist, a name already used and one too long; - the counts "10", "5 3", "x 5" and "3 0", then a cancel; - removing a group that does not exist, then all but one, which has to stay; with one group there is no paradox.
Built on the verified corpus cases simpsons-paradox (rates compared by cross-multiplication, so the reversal is exact) and average-of-averages (the mean of ratios against the ratio of totals).
Recorded sessions
What the screen shows while someone uses the program. Each typed line appears after its prompt, the way a terminal shows it.
bad-input
interpreter: byte-equal== Simpson's paradox ==
Two treatments, A and B: better in every group, and worse all together.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 0
Pick a number from 1 to 7.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> x
Pick a number from 1 to 7.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 4
group name> tiny
No group is called tiny.
group name>
Cancelled.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 5
group name> small stones
There is a group called small stones already.
group name> a name far too long
A name is up to 16 characters.
group name> new
A: successes and cases> 10
Type the successes and the cases, like 81 87: whole numbers, at least 1 case, no more successes than cases.
A: successes and cases> 5 3
Type the successes and the cases, like 81 87: whole numbers, at least 1 case, no more successes than cases.
A: successes and cases> x 5
Type the successes and the cases, like 81 87: whole numbers, at least 1 case, no more successes than cases.
A: successes and cases> 3 0
Type the successes and the cases, like 81 87: whole numbers, at least 1 case, no more successes than cases.
A: successes and cases>
Cancelled.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 6
group name> nothing
No group is called nothing.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 6
group name> large stones
large stones removed.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 6
One group has to stay.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 1
group A B better
small stones 81/87 93.1% 234/270 86.7% A
all together 81/87 93.1% 234/270 86.7% A
No paradox here: the pooled rates do not turn against the groups.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 7
Bye.
What was typed (21 lines)
0
x
4
tiny
5
small stones
a name far too long
new
10
5 3
x 5
3 0
6
nothing
6
large stones
6
1
7
basic
interpreter: byte-equal== Simpson's paradox ==
Two treatments, A and B: better in every group, and worse all together.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 1
group A B better
small stones 81/87 93.1% 234/270 86.7% A
large stones 192/263 73.0% 55/80 68.8% A
all together 273/350 78.0% 289/350 82.6% B
Simpson's paradox: A does better in every group, and B does better with the groups put together.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 2
Where each treatment's cases came from:
small stones A 24.9% of its cases B 77.1% of its cases
large stones A 75.1% of its cases B 22.9% of its cases
A pooled rate is each group's rate weighted by how many of that arm's cases the group holds - so it depends on who got which treatment.
Weighing the groups the same for both (each by its share of all 700 cases): A 83.3%, B 77.9% - better: A.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 3
Average of the group rates: A 83.1%, B 77.7%.
Pooled rate, every case counted once: A 78.0%, B 82.6%.
The average counts every group once, whatever its size; the pooled rate counts every case once. They agree only when the groups are the same size.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 4
group name> large stones
A: successes and cases> 192 263
B: successes and cases> 80 120
large stones changed.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 1
group A B better
small stones 81/87 93.1% 234/270 86.7% A
large stones 192/263 73.0% 80/120 66.7% A
all together 273/350 78.0% 314/390 80.5% B
Simpson's paradox: A does better in every group, and B does better with the groups put together.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 5
group name> medium stones
A: successes and cases> 50 60
B: successes and cases> 40 50
medium stones added.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 1
group A B better
small stones 81/87 93.1% 234/270 86.7% A
large stones 192/263 73.0% 80/120 66.7% A
medium stones 50/60 83.3% 40/50 80.0% A
all together 323/410 78.8% 354/440 80.5% B
Simpson's paradox: A does better in every group, and B does better with the groups put together.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 6
group name> medium stones
medium stones removed.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 1
group A B better
small stones 81/87 93.1% 234/270 86.7% A
large stones 192/263 73.0% 80/120 66.7% A
all together 273/350 78.0% 314/390 80.5% B
Simpson's paradox: A does better in every group, and B does better with the groups put together.
1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit
choice> 7
Bye.
What was typed (17 lines)
1
2
3
4
large stones
192 263
80 120
1
5
medium stones
50 60
40 50
1
6
medium stones
1
7
Modules
The program as written, entry module first. Each module transpiles to its own Python file, which is what eml project run executes.
main.eml(entry)
eml# P086 Simpson's paradox: two treatments, their success rates group by
# group and pooled, why the pooled comparison can turn round, the average
# of the rates against the pooled rate, and a comparison that weighs the
# groups the same for both.
import rates
def trim(s):
0 => i
len(s) => j
while i < j and s[i] == " ":
i + 1 => i
while j > i and s[j - 1] == " ":
j - 1 => j
return s[i:j]
def number(s):
if len(s) == 0 or len(s) > 6:
return -1
0 => n
for c in s:
if not (c >= "0" and c <= "9"):
return -1
n * 10 + int(c) => n
return n
def pad(s, n):
while len(s) < n:
s + " " => s
return s
def lpad(s, n):
while len(s) < n:
" " + s => s
return s
def who(v):
if v == 1:
return "A"
if v == -1:
return "B"
return "tie"
def cell(ok, n):
return lpad(str(ok) + "/" + str(n), 11) + " " + lpad(rates.percent(ok, n), 6)
def table(groups):
rates.verdicts(groups) => v
(" " + pad("group", 16) + lpad("A", 18) + lpad("B", 18) + " better") ^0
for i in [0:len(groups) - 1]:
groups[i] => g
(" " + pad(g[0], 16) + cell(g[1], g[2]) + cell(g[3], g[4]) + " " + who(v[0][i])) ^0
rates.totals(groups) => t
(" " + pad("all together", 16) + cell(t[0], t[1]) + cell(t[2], t[3]) + " " + who(v[1])) ^0
if v[2]:
("Simpson's paradox: " + who(0 - v[1]) + " does better in every group, and " + who(v[1]) + " does better with the groups put together.") ^0
else:
"No paradox here: the pooled rates do not turn against the groups." ^0
def why(groups):
rates.totals(groups) => t
"Where each treatment's cases came from:" ^0
for g in groups:
(" " + pad(g[0], 16) + "A " + lpad(rates.percent(g[2], t[1]), 6) + " of its cases B " + lpad(rates.percent(g[4], t[3]), 6) + " of its cases") ^0
"A pooled rate is each group's rate weighted by how many of that arm's cases the group holds - so it depends on who got which treatment." ^0
rates.standardized(groups, 1) => sa
rates.standardized(groups, 3) => sb
rates.better(sa[0], sa[1], sb[0], sb[1]) => v
("Weighing the groups the same for both (each by its share of all " + str(t[1] + t[3]) + " cases): A " + rates.percent(sa[0], sa[1]) + ", B " + rates.percent(sb[0], sb[1]) + " - better: " + who(v) + ".") ^0
def means(groups):
rates.totals(groups) => t
rates.mean_of_rates(groups, 1) => ma
rates.mean_of_rates(groups, 3) => mb
("Average of the group rates: A " + rates.percent(ma[0], ma[1]) + ", B " + rates.percent(mb[0], mb[1]) + ".") ^0
("Pooled rate, every case counted once: A " + rates.percent(t[0], t[1]) + ", B " + rates.percent(t[2], t[3]) + ".") ^0
"The average counts every group once, whatever its size; the pooled rate counts every case once. They agree only when the groups are the same size." ^0
def ask_counts(prompt):
# [successes, cases] from "81 87", or [] if cancelled.
while True:
trim(input(prompt)) => s
if s == "":
return []
0 => k
while k < len(s) and s[k] != " ":
k + 1 => k
number(trim(s[0:k])) => ok
-1 => n
if k < len(s):
number(trim(s[k + 1:len(s)])) => n
if ok >= 0 and n >= 1 and ok <= n:
return [ok, n]
"Type the successes and the cases, like 81 87: whole numbers, at least 1 case, no more successes than cases." ^0
def find(groups, name):
for i in [0:len(groups) - 1]:
if groups[i][0] == name:
return i
return -1
def edit(groups, adding):
while True:
trim(input("group name> ")) => name
if name == "":
"Cancelled." ^0
return groups
find(groups, name) => i
if len(name) > 16:
"A name is up to 16 characters." ^0
elif adding and i >= 0:
("There is a group called " + name + " already.") ^0
elif not adding and i == -1:
("No group is called " + name + ".") ^0
elif adding and len(groups) == 6:
"There are 6 groups already." ^0
return groups
else:
ask_counts("A: successes and cases> ") => a
if len(a) == 0:
"Cancelled." ^0
return groups
ask_counts("B: successes and cases> ") => b
if len(b) == 0:
"Cancelled." ^0
return groups
[name, a[0], a[1], b[0], b[1]] => g
if adding:
(name + " added.") ^0
return groups + [g]
[] => out
for h in groups:
if h[0] == name:
out + [g] => out
else:
out + [h] => out
(name + " changed.") ^0
return out
def remove(groups):
if len(groups) == 1:
"One group has to stay." ^0
return groups
trim(input("group name> ")) => name
if name == "":
"Cancelled." ^0
return groups
[] => out
for h in groups:
if h[0] != name:
out + [h] => out
if len(out) == len(groups):
("No group is called " + name + ".") ^0
return groups
(name + " removed.") ^0
return out
[["small stones", 81, 87, 234, 270], ["large stones", 192, 263, 55, 80]] => groups
"== Simpson's paradox ==" ^0
"Two treatments, A and B: better in every group, and worse all together." ^0
True => running
while running:
"" ^0
"1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit" ^0
trim(input("choice> ")) => choice
if choice == "1":
table(groups)
elif choice == "2":
why(groups)
elif choice == "3":
means(groups)
elif choice == "4":
edit(groups, False) => groups
elif choice == "5":
edit(groups, True) => groups
elif choice == "6":
remove(groups) => groups
elif choice == "7":
False => running
else:
"Pick a number from 1 to 7." ^0
"Bye." ^0
Python projection (main.py)
import rates
def trim(s):
i = 0
j = len(s)
while i < j and s[i] == " ":
i = i + 1
while j > i and s[j - 1] == " ":
j = j - 1
return s[i:j]
def number(s):
if len(s) == 0 or len(s) > 6:
return -1
n = 0
for c in s:
if not (c >= "0" and c <= "9"):
return -1
n = n * 10 + int(c)
return n
def pad(s, n):
while len(s) < n:
s = s + " "
return s
def lpad(s, n):
while len(s) < n:
s = " " + s
return s
def who(v):
if v == 1:
return "A"
if v == -1:
return "B"
return "tie"
def cell(ok, n):
return lpad(str(ok) + "/" + str(n), 11) + " " + lpad(rates.percent(ok, n), 6)
def table(groups):
v = rates.verdicts(groups)
print(" " + pad("group", 16) + lpad("A", 18) + lpad("B", 18) + " better")
for i in range(0, len(groups)):
g = groups[i]
print(" " + pad(g[0], 16) + cell(g[1], g[2]) + cell(g[3], g[4]) + " " + who(v[0][i]))
t = rates.totals(groups)
print(" " + pad("all together", 16) + cell(t[0], t[1]) + cell(t[2], t[3]) + " " + who(v[1]))
if v[2]:
print("Simpson's paradox: " + who(0 - v[1]) + " does better in every group, and " + who(v[1]) + " does better with the groups put together.")
else:
print("No paradox here: the pooled rates do not turn against the groups.")
def why(groups):
t = rates.totals(groups)
print("Where each treatment's cases came from:")
for g in groups:
print(" " + pad(g[0], 16) + "A " + lpad(rates.percent(g[2], t[1]), 6) + " of its cases B " + lpad(rates.percent(g[4], t[3]), 6) + " of its cases")
print("A pooled rate is each group's rate weighted by how many of that arm's cases the group holds - so it depends on who got which treatment.")
sa = rates.standardized(groups, 1)
sb = rates.standardized(groups, 3)
v = rates.better(sa[0], sa[1], sb[0], sb[1])
print("Weighing the groups the same for both (each by its share of all " + str(t[1] + t[3]) + " cases): A " + rates.percent(sa[0], sa[1]) + ", B " + rates.percent(sb[0], sb[1]) + " - better: " + who(v) + ".")
def means(groups):
t = rates.totals(groups)
ma = rates.mean_of_rates(groups, 1)
mb = rates.mean_of_rates(groups, 3)
print("Average of the group rates: A " + rates.percent(ma[0], ma[1]) + ", B " + rates.percent(mb[0], mb[1]) + ".")
print("Pooled rate, every case counted once: A " + rates.percent(t[0], t[1]) + ", B " + rates.percent(t[2], t[3]) + ".")
print("The average counts every group once, whatever its size; the pooled rate counts every case once. They agree only when the groups are the same size.")
def ask_counts(prompt):
while True:
s = trim(input(prompt))
if s == "":
return []
k = 0
while k < len(s) and s[k] != " ":
k = k + 1
ok = number(trim(s[0:k]))
n = -1
if k < len(s):
n = number(trim(s[k + 1:len(s)]))
if ok >= 0 and n >= 1 and ok <= n:
return [ok, n]
print("Type the successes and the cases, like 81 87: whole numbers, at least 1 case, no more successes than cases.")
def find(groups, name):
for i in range(0, len(groups)):
if groups[i][0] == name:
return i
return -1
def edit(groups, adding):
while True:
name = trim(input("group name> "))
if name == "":
print("Cancelled.")
return groups
i = find(groups, name)
if len(name) > 16:
print("A name is up to 16 characters.")
elif adding and i >= 0:
print("There is a group called " + name + " already.")
elif not adding and i == -1:
print("No group is called " + name + ".")
elif adding and len(groups) == 6:
print("There are 6 groups already.")
return groups
else:
a = ask_counts("A: successes and cases> ")
if len(a) == 0:
print("Cancelled.")
return groups
b = ask_counts("B: successes and cases> ")
if len(b) == 0:
print("Cancelled.")
return groups
g = [name, a[0], a[1], b[0], b[1]]
if adding:
print(name + " added.")
return groups + [g]
out = []
for h in groups:
if h[0] == name:
out = out + [g]
else:
out = out + [h]
print(name + " changed.")
return out
def remove(groups):
if len(groups) == 1:
print("One group has to stay.")
return groups
name = trim(input("group name> "))
if name == "":
print("Cancelled.")
return groups
out = []
for h in groups:
if h[0] != name:
out = out + [h]
if len(out) == len(groups):
print("No group is called " + name + ".")
return groups
print(name + " removed.")
return out
groups = [["small stones", 81, 87, 234, 270], ["large stones", 192, 263, 55, 80]]
print("== Simpson's paradox ==")
print("Two treatments, A and B: better in every group, and worse all together.")
running = True
while running:
print("")
print("1) table 2) why 3) average of rates 4) change a group 5) add a group 6) remove a group 7) quit")
choice = trim(input("choice> "))
if choice == "1":
table(groups)
elif choice == "2":
why(groups)
elif choice == "3":
means(groups)
elif choice == "4":
groups = edit(groups, False)
elif choice == "5":
groups = edit(groups, True)
elif choice == "6":
groups = remove(groups)
elif choice == "7":
running = False
else:
print("Pick a number from 1 to 7.")
print("Bye.")
rates.eml
eml# P086 Simpson's paradox - success rates of two treatments, A and B, group
# by group and pooled. Every rate is a fraction of two whole numbers and
# rates are compared by cross-multiplication, as in the corpus case
# simpsons-paradox, so no rounding can cause a reversal. A group is
# [name, A successes, A cases, B successes, B cases].
def better(a_ok, a_n, b_ok, b_n):
# 1 if A's rate beats B's, -1 if B's beats A's, 0 if equal:
# a_ok/a_n > b_ok/b_n exactly when a_ok b_n > b_ok a_n.
a_ok * b_n => l
b_ok * a_n => r
if l > r:
return 1
if l < r:
return -1
return 0
def percent(num, den):
# num/den as a percentage to one decimal, rounded half up, by long
# division in whole numbers.
2000 * num + den => n
2 * den => d
0 => q
while n >= d:
d => m
1 => k
while m + m <= n:
m + m => m
k + k => k
n - m => n
q + k => q
0 => w
q => r
while r >= 10:
r - 10 => r
w + 1 => w
return str(w) + "." + str(r) + "%"
def totals(groups):
# [A successes, A cases, B successes, B cases] over all groups.
[0, 0, 0, 0] => t
for g in groups:
[t[0] + g[1], t[1] + g[2], t[2] + g[3], t[3] + g[4]] => t
return t
def verdicts(groups):
# [winner in each group (1 A, -1 B, 0 tie), winner pooled, paradox]:
# a paradox is every group going one way (ties aside, at least one
# group won) and the pooled rates going the other way.
[] => each
0 => a_wins
0 => b_wins
for g in groups:
better(g[1], g[2], g[3], g[4]) => v
each + [v] => each
if v == 1:
a_wins + 1 => a_wins
elif v == -1:
b_wins + 1 => b_wins
totals(groups) => t
better(t[0], t[1], t[2], t[3]) => pooled
((a_wins > 0 and b_wins == 0 and pooled == -1) or (b_wins > 0 and a_wins == 0 and pooled == 1)) => paradox
return [each, pooled, paradox]
def mean_of_rates(groups, arm):
# The plain average of one arm's group rates (arm 1 is A, 3 is B), as a
# fraction: each group counts once, whatever its size - the corpus case
# average-of-averages.
[0, 1] => s
for g in groups:
[s[0] * g[arm + 1] + g[arm] * s[1], s[1] * g[arm + 1]] => s
return [s[0], s[1] * len(groups)]
def standardized(groups, arm):
# One arm's rate if it had treated each group in the same proportion -
# the share of all cases (both arms) that group has: sum over groups of
# rate x group size, over all cases. Comparing these compares like with
# like.
0 => all_cases
for g in groups:
all_cases + g[2] + g[4] => all_cases
[0, 1] => s
for g in groups:
g[2] + g[4] => w
[s[0] * g[arm + 1] + g[arm] * w * s[1], s[1] * g[arm + 1]] => s
return [s[0], s[1] * all_cases]
Python projection (rates.py)
def better(a_ok, a_n, b_ok, b_n):
l = a_ok * b_n
r = b_ok * a_n
if l > r:
return 1
if l < r:
return -1
return 0
def percent(num, den):
n = 2000 * num + den
d = 2 * den
q = 0
while n >= d:
m = d
k = 1
while m + m <= n:
m = m + m
k = k + k
n = n - m
q = q + k
w = 0
r = q
while r >= 10:
r = r - 10
w = w + 1
return str(w) + "." + str(r) + "%"
def totals(groups):
t = [0, 0, 0, 0]
for g in groups:
t = [t[0] + g[1], t[1] + g[2], t[2] + g[3], t[3] + g[4]]
return t
def verdicts(groups):
each = []
a_wins = 0
b_wins = 0
for g in groups:
v = better(g[1], g[2], g[3], g[4])
each = each + [v]
if v == 1:
a_wins = a_wins + 1
elif v == -1:
b_wins = b_wins + 1
t = totals(groups)
pooled = better(t[0], t[1], t[2], t[3])
paradox = a_wins > 0 and b_wins == 0 and pooled == -1 or b_wins > 0 and a_wins == 0 and pooled == 1
return [each, pooled, paradox]
def mean_of_rates(groups, arm):
s = [0, 1]
for g in groups:
s = [s[0] * g[arm + 1] + g[arm] * s[1], s[1] * g[arm + 1]]
return [s[0], s[1] * len(groups)]
def standardized(groups, arm):
all_cases = 0
for g in groups:
all_cases = all_cases + g[2] + g[4]
s = [0, 1]
for g in groups:
w = g[2] + g[4]
s = [s[0] * g[arm + 1] + g[arm] * w * s[1], s[1] * g[arm + 1]]
return [s[0], s[1] * all_cases]