Lesson 10 · Statistics for data analysts

Multiple Comparisons & P-Value Corrections | Statistics #10

Video ten of the 15-part series: what happens to your real false-positive rate once you stop running one test and start running many. Real Sample Superstore…

⬇ Download notebook ⬇ Notebook + data (.zip)Open in Colab ↗

CourseStatistics for data analysts
Lesson10 of 15
Video13 min
FormatJupyter notebook · 8 code cells
Data1 dataset

statsmodels pandas NumPy SciPy itertools

What you'll learn

The Problem, in Theory and in Real Data
Running It for Real - 136 Sub-Category Comparisons
The Bonferroni Correction
Tukey's HSD and the Benjamini-Hochberg FDR
Choosing a Correction
Wrap-Up: What You Learned

Datasets used in this lesson

Save these next to the notebook. In Google Colab, upload them with the 📁 icon on the left first.

superstore_sales.csv715.5 KB

📓 Full notebook

Download .ipynb

Statistics for Data Analysts, Video 10: Multiple Comparisons#

Video ten of the 15-part series: what happens to your real false-positive rate once you stop running one test and start running many.
Real Sample Superstore order data throughout, including 136 real pairwise sub-category comparisons.
Let's get into it.

Before You Start#

Open a new Jupyter Notebook in VS Code and select your Python interpreter as the kernel.
Install statsmodels if you don't have it yet: pip install statsmodels.
Place superstore_sales.csv in the same folder as this notebook.

import pandas as pd
import numpy as np
from scipy import stats
from itertools import combinations
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.multicomp import pairwise_tukeyhsd

df = pd.read_csv('superstore_sales.csv')
print(df.shape)

(9994, 9)

Part 1: The Problem, in Theory and in Real Data#

for k in [1, 5, 10, 20, 50, 136]:
    fwer = 1 - (0.95) ** k
    print(f'{k} independent tests at alpha=0.05: at least one false positive with probability {fwer:.1%}')

1 independent tests at alpha=0.05: at least one false positive with probability 5.0%
5 independent tests at alpha=0.05: at least one false positive with probability 22.6%
10 independent tests at alpha=0.05: at least one false positive with probability 40.1%
20 independent tests at alpha=0.05: at least one false positive with probability 64.2%
50 independent tests at alpha=0.05: at least one false positive with probability 92.3%
136 independent tests at alpha=0.05: at least one false positive with probability 99.9%

office = df[df['Category'] == 'Office Supplies']['Sales'].values
rng = np.random.default_rng(seed=9)
shuffled = rng.permutation(office)
null_groups = np.array_split(shuffled, 20)
null_pairs = list(combinations(range(20), 2))
null_pvals = np.array([stats.ttest_ind(null_groups[i], null_groups[j])[1] for i, j in null_pairs])
print(f'Real null-world pairs tested: {len(null_pairs)} (all drawn from the same real Office Supplies population)')
print(f'Real false positives at raw alpha=0.05: {(null_pvals < 0.05).sum()}')

Real null-world pairs tested: 190 (all drawn from the same real Office Supplies population)
Real false positives at raw alpha=0.05: 2

Part 2: Running It for Real - 136 Sub-Category Comparisons#

subcats = df['Sub-Category'].unique()
pairs = list(combinations(subcats, 2))
pvals = np.array([stats.mannwhitneyu(df[df['Sub-Category'] == a]['Sales'], df[df['Sub-Category'] == b]['Sales']).pvalue for a, b in pairs])
print(f'Real sub-categories: {len(subcats)}')
print(f'Real pairwise comparisons: {len(pairs)}')
print(f'Real significant at raw alpha=0.05: {(pvals < 0.05).sum()} of {len(pairs)}')

Real sub-categories: 17
Real pairwise comparisons: 136
Real significant at raw alpha=0.05: 131 of 136

Part 3: The Bonferroni Correction#

bonferroni_alpha = 0.05 / len(pairs)
bonf_significant = pvals < bonferroni_alpha
print(f'Real Bonferroni-adjusted alpha: {bonferroni_alpha:.6f}')
print(f'Real significant after Bonferroni: {bonf_significant.sum()} of {len(pairs)}')

Real Bonferroni-adjusted alpha: 0.000368
Real significant after Bonferroni: 123 of 136

flipped = [(pairs[i], pvals[i]) for i in range(len(pairs)) if pvals[i] < 0.05 and pvals[i] >= bonferroni_alpha]
print(f'Real pairs significant raw but NOT after Bonferroni: {len(flipped)}')
for pair, p in sorted(flipped, key=lambda x: x[1])[:5]:
    print(f'  {pair}: raw p={p:.5f}')

Real pairs significant raw but NOT after Bonferroni: 8
  ('Furnishings', 'Envelopes'): raw p=0.00072
  ('Art', 'Binders'): raw p=0.00095
  ('Binders', 'Supplies'): raw p=0.00096
  ('Machines', 'Copiers'): raw p=0.00096
  ('Tables', 'Machines'): raw p=0.00216

Part 4: Tukey's HSD and the Benjamini-Hochberg FDR#

tukey_data = df[['Sales', 'Sub-Category']].dropna()
tukey_result = pairwise_tukeyhsd(tukey_data['Sales'], tukey_data['Sub-Category'], alpha=0.05)
print(tukey_result.summary().as_text()[:1200])

         Multiple Comparison of Means - Tukey HSD, FWER=0.05          
======================================================================
   group1      group2    meandiff  p-adj    lower      upper    reject
----------------------------------------------------------------------
Accessories  Appliances    14.7811    1.0   -98.3524   127.9146  False
Accessories         Art  -181.9058    0.0  -279.2992   -84.5123   True
Accessories     Binders    -82.414 0.0708  -167.5716     2.7435  False
Accessories   Bookcases    287.885    0.0   142.4794   433.2907   True
Accessories      Chairs   316.3578    0.0    212.228   420.4877   True
Accessories     Copiers   1982.967    0.0  1738.8728  2227.0613   True
Accessories   Envelopes  -151.1069 0.0188  -290.6438     -11.57   True
Accessories   Fasteners  -202.0378 0.0003  -350.2638   -53.8119   True
Accessories Furnishings  -120.1489  0.001  -213.4134   -26.8845   True
Accessories      Labels  -181.6715    0.0  -304.3051    -59.038   True
Accessories    Machines  1429.5787    0.0  1236.7178  1622.4397   True
Accessories       Paper  -158.6905    0.0  -245.4369   -71.9441   True
Accessories      Phones   155.2369    0.0    60.3898    250.084

reject_bh, pvals_bh, _, _ = multipletests(pvals, alpha=0.05, method='fdr_bh')
print(f'Real significant after Benjamini-Hochberg FDR correction: {reject_bh.sum()} of {len(pairs)}')
print(f'Real significant after Bonferroni: {bonf_significant.sum()} of {len(pairs)}')
print(f'Real significant, raw uncorrected: {(pvals < 0.05).sum()} of {len(pairs)}')

Real significant after Benjamini-Hochberg FDR correction: 131 of 136
Real significant after Bonferroni: 123 of 136
Real significant, raw uncorrected: 131 of 136

Part 5: Choosing a Correction#

Bonferroni: simplest, most conservative, controls the chance of even one false positive across the whole family of tests. Good default when a single false alarm would be costly.
Tukey's HSD: purpose-built for comparing every pair of groups after an ANOVA-style design; use it specifically for that all-pairs case.
Benjamini-Hochberg (FDR): less conservative, controls the expected proportion of false positives among the significant results rather than eliminating false positives entirely. Common choice when screening many real comparisons and some real true effects are expected, like this sub-category analysis.
Whichever correction you pick, decide on it before looking at the real results, exactly the same discipline as the pre-registration and sample-size planning from video eight.

Wrap-Up: What You Learned#

Why the family-wise error rate grows with the number of tests, in theory and confirmed with a real null-world simulation.
Running 136 real pairwise sub-category comparisons and seeing how the significant count shifts under different corrections.
The Bonferroni correction, and the specific real borderline pairs it flips from significant to not.
Tukey's HSD, purpose-built for all-pairs comparisons, and the Benjamini-Hochberg FDR correction as a middle ground.
A practical guide for choosing a correction method before looking at real results.
Video eleven shifts from comparing groups to modeling relationships directly, starting with simple linear regression on real mtcars data. Subscribe so it lands automatically see you there.

Found this useful?

All lessons, notebooks and datasets here are free. If they helped you, a coffee keeps new lessons coming.

☕ Buy me a coffee Other ways to help