Common Statistical Mistakes

Most wrong conclusions in reports, news stories and even published research do not come from arithmetic errors. They come from a few classic traps: testing until something looks significant, adding up groups that should be kept apart, studying only the winners, and mistaking luck for skill. Each trap looks convincing from the inside.

In this lesson you will meet four famous traps, p-hacking, Simpson's paradox, survivorship bias and regression to the mean, reproduce each one with a short simulation or table so you can recognise it, and finish with a checklist of habits that keep your own work honest.

p-hacking and multiple comparisons

If you run one test at α = 0.05 on data with no real effect, you have a 5% chance of a false alarm. If you run 20 independent tests, the chance that at least one is "significant" is:

P(at least one false alarm) = 1 - 0.95^20 = 0.64

p-hacking means trying many analyses (different variables, subgroups, outlier rules, time windows) and reporting the ones that cross p < 0.05. Often it is not deliberate cheating, just curiosity without discipline, but the result is the same: findings that are pure noise.

import numpy as np

rng = np.random.default_rng(0)
studies = rng.random(size=(100_000, 20))       # under H0, p-values are uniform
print("Theory: 1 - 0.95^20 =", round(1 - 0.95 ** 20, 4))
print("Simulated:", (studies < 0.05).any(axis=1).mean())

Theory: 1 - 0.95^20 = 0.6415
Simulated: 0.64217

When H0 is true, p-values are spread evenly between 0 and 1 (you saw this in the hypothesis testing lesson), so random numbers stand in for them. In about 64% of "studies" that run 20 tests, at least one test is significant by chance alone.

Here is the same trap with data. An invented survey records 20 daily habits of 100 students, none of which has anything to do with their marks. Test each one:

import numpy as np
from scipy import stats

rng = np.random.default_rng(5)
marks = rng.normal(60, 12, size=100)            # invented marks of 100 students
habits = rng.normal(size=(100, 20))             # 20 "habits" unrelated to marks
pvals = np.array([stats.pearsonr(habits[:, j], marks).pvalue for j in range(20)])

print("Smallest p-values:", np.sort(pvals)[:3].round(4))
print("'Significant' at 0.05:", (pvals < 0.05).sum(), "of 20")
print("After Bonferroni (p < 0.05 / 20):", (pvals < 0.05 / 20).sum())
print("After Benjamini-Hochberg:", (stats.false_discovery_control(pvals) < 0.05).sum())

Smallest p-values: [0.0026 0.0156 0.1043]
'Significant' at 0.05: 2 of 20
After Bonferroni (p < 0.05 / 20): 0
After Benjamini-Hochberg: 0

Two habits look "significant", one with p = 0.0026, which would make an exciting headline. Both are noise. Two standard corrections remove them:

• Bonferroni: with m tests, require p < α / m (here 0.05 / 20 = 0.0025). Simple and strict.
• Benjamini-Hochberg: controls the false discovery rate, the expected share of false findings among your "discoveries". It is less strict than Bonferroni and popular when you test many things. stats.false_discovery_control(pvals) returns adjusted p-values to compare with α.

The best protection is to decide your main question and test before you look at the data, and to treat everything else as exploratory ideas to check with new data.

Simpson's paradox

Simpson's paradox happens when a trend that appears in every group reverses (or disappears) when the groups are combined. Two coaching centres report their pass rates (invented numbers):

import pandas as pd

data = pd.DataFrame({
    "centre": ["A", "A", "B", "B"],
    "group": ["Strong", "Weak", "Strong", "Weak"],
    "students": [20, 180, 170, 30],
    "passed": [19, 96, 153, 15],
})
data["rate"] = (data["passed"] / data["students"]).round(3)
print(data.pivot(index="group", columns="centre", values="rate"))
overall = data.groupby("centre")[["passed", "students"]].sum()
print((overall["passed"] / overall["students"]).round(3))

centre      A    B
group             
Strong  0.950  0.9
Weak    0.533  0.5
centre
A    0.575
B    0.840
dtype: float64

Centre A has the higher pass rate for strong students (95% vs 90%) and for weak students (53% vs 50%). Yet overall, B looks far better (84% vs 57.5%). The reason: A mostly teaches weak students (180 of its 200), while B mostly teaches strong ones. The overall rate mixes "which centre" with "which students", a confounder you met in the correlation lesson. The fair comparison is within groups: A does better. Whenever groups differ in their make-up, compare like with like, or adjust with multiple regression.

Survivorship bias

Survivorship bias means drawing conclusions only from the cases that "survived" some filter, while the failures are invisible. Famous examples: studying only successful startups to find the secret of success, or praising old buildings for being well built when the badly built ones fell down long ago.

During the Second World War, the statistician Abraham Wald is said to have pointed out that the bullet holes on returning bombers showed where a plane could be hit and still fly home; the armour belonged where the returning planes had no holes. Here is a modern version with invented investment funds that have no real skill at all:

import numpy as np

rng = np.random.default_rng(4)
funds = 1000
returns = rng.normal(0.0, 0.15, size=(funds, 5))     # 5 yearly returns, no real skill
value = (1 + returns).cumprod(axis=1)
closed = (value < 0.8).any(axis=1)                    # shut down after losing 20%
print(f"Funds closed: {closed.sum()} of {funds}")
print(f"Average yearly return, all funds:      {returns.mean():+.2%}")
print(f"Average yearly return, survivors only: {returns[~closed].mean():+.2%}")

Funds closed: 381 of 1000
Average yearly return, all funds:      +0.20%
Average yearly return, survivors only: +3.92%

Every fund's returns are pure luck, averaging about zero. But 381 funds that fell 20% below their starting value were closed and dropped from the record, and the survivors show an average of almost +4% a year. A database of "current funds" would make the whole industry look skilful. Always ask: who or what is missing from this data, and why?

Regression to the mean

When a measurement mixes skill and luck, extreme results are usually partly good (or bad) luck, and luck does not repeat. So the next measurement of an extreme case tends to be closer to the average. This is regression to the mean, and it fools people constantly.

import numpy as np

rng = np.random.default_rng(8)
players = 200
skill = rng.normal(30, 4, size=players)              # true batting average
season1 = skill + rng.normal(0, 8, size=players)     # skill plus luck
season2 = skill + rng.normal(0, 8, size=players)
top = np.argsort(season1)[-10:]
print(f"Top 10 of season 1: season 1 average {season1[top].mean():.1f}")
print(f"Same 10 players:    season 2 average {season2[top].mean():.1f}")
print(f"Their true skill:   {skill[top].mean():.1f}  (all players: {skill.mean():.1f})")

Top 10 of season 1: season 1 average 48.5
Same 10 players:    season 2 average 33.2
Their true skill:   32.6  (all players: 30.1)

The ten best batters of season 1 averaged 48.5, but their real skill is only 32.6: they were good and lucky. In season 2 their luck is new, and they average 33.2. Nothing went wrong with them. Yet commentators would invent a reason: pressure, a "jinx", complacency. The same effect makes punishments seem to work (a terrible result is usually followed by a better one anyway) and praise seem to fail. To show that a change really helps, compare with a control group that did not get it.

More traps to watch for

• Base rate neglect: forgetting how rare something is, as in the screening test of the Bayes lesson.
• Correlation read as causation, without asking about confounders or reverse causation.
• Tiny samples: a 100% success rate from 3 cases means very little; show the confidence interval.
• Misleading charts: bar charts that do not start at zero, cherry-picked time windows, or two y-axes chosen to make lines cross.
• Extrapolation: using a regression line far outside the data.
• Ecological fallacy: concluding about individuals from group averages. A district with higher average income and higher average marks does not show that richer students in it score higher.
• Confusing SD and SE, or percentage points and percent.

Common mistakes

• Running many tests and reporting the winners. Pre-register the main question, or correct for multiple comparisons.
• Comparing totals across groups with different make-up. Compare within strata.
• Only looking at survivors. Find out what was filtered out before the data reached you.
• Explaining away regression to the mean. Use a control group before claiming an effect.

Key points

• With 20 tests at α = 0.05, a false alarm somewhere is more likely than not; use Bonferroni or stats.false_discovery_control and decide tests in advance.
• Simpson's paradox: a trend inside every group can reverse when groups of different make-up are combined.
• Survivorship bias: studying only what passed a filter makes results look better than reality.
• Regression to the mean: extreme results mix skill and luck, so they tend to move towards the average next time.
• Good habits: plan, report everything you tested, show effect sizes and intervals, and look for missing data.

import numpy as np
from scipy import stats

rng = np.random.default_rng(11)
n = 600
tuition = rng.random(n) < 0.4
marks = rng.normal(60, 15, size=n)            # no real tuition effect anywhere
district = rng.integers(0, 12, size=n)

pvals = []
for d in range(12):
    here = district == d
    res = stats.ttest_ind(marks[here & tuition], marks[here & ~tuition], equal_var=False)
    pvals.append(res.pvalue)
pvals = np.array(pvals)

print("Smallest p-value:", round(pvals.min(), 4))
print("Districts with p < 0.05:", (pvals < 0.05).sum())
print("After Benjamini-Hochberg:", (stats.false_discovery_control(pvals) < 0.05).sum())

Smallest p-value: 0.0277
Districts with p < 0.05: 1
After Benjamini-Hochberg: 0