#!/usr/bin/env python3
"""Holm/Bonferroni/BH worked example (LAT-2297, facelift topic 3).

REAL DATA. UCI Wine Quality (red), 1,599 wines, 12 chemistry variables, all 66
pairwise Pearson correlations. No simulation, no construction. This is the same
corpus standard_multiple_comparisons is verified against, so the video, the Rmd
and the live tool all read one set of numbers.

WHAT IT SHOWS, and it is not what I expected to find:

    raw p<.05     55 of 66
    Bonferroni    43
    Holm          43      <- identical
    BH (FDR)      54

The chooser is NOT Bonferroni vs Holm. That choice moves nothing here. The
chooser is FWER vs FDR, and it moves eleven discoveries out of 66.

Holm is still the strictly better of the two FWER methods: same guarantee,
uniformly >= Bonferroni, never worse. Use it. Just do not expect it to rescue
an underpowered study -- it gains only when a p-value falls in the narrow band
between alpha/m and alpha/(m-k), and on real data that band is often empty.
Here it was nearly not: rank 44 sits at p=2.214e-03 against a Holm bar of
2.174e-03. It missed gaining one discovery by 4e-05.

Run:  python3 holm_example.py       (writes wine_correlation_pvalues.csv)
"""
import csv, io, itertools, urllib.request
import numpy as np
from scipy import stats

URL = ("https://archive.ics.uci.edu/ml/machine-learning-databases/"
       "wine-quality/winequality-red.csv")
ALPHA = 0.05

raw = urllib.request.urlopen(URL, timeout=60).read().decode()
rd = list(csv.reader(io.StringIO(raw), delimiter=";"))
hdr = [h.strip('"') for h in rd[0]]
data = np.array([[float(x) for x in r] for r in rd[1:] if r])

rows = []
for i, j in itertools.combinations(range(len(hdr)), 2):
    r, p = stats.pearsonr(data[:, i], data[:, j])
    rows.append((f"{hdr[i]} ~ {hdr[j]}", float(p), float(r)))
rows.sort(key=lambda t: t[1])

with open("wine_correlation_pvalues.csv", "w", newline="") as f:
    w = csv.writer(f)
    w.writerow(["comparison", "p_value"])
    for name, p, _ in rows:
        w.writerow([name, f"{p:.6g}"])

ps = np.array([p for _, p, _ in rows]); m = len(ps)
raw_n = int((ps < ALPHA).sum())
bonf = int((ps < ALPHA / m).sum())
holm = 0
for k, p in enumerate(ps):
    if p < ALPHA / (m - k): holm += 1
    else: break
bh = 0
for k, p in enumerate(ps, 1):
    if p <= k * ALPHA / m: bh = k

print(f"n wines {data.shape[0]}   variables {len(hdr)}   pairwise tests {m}")
print(f"  raw p<{ALPHA}   {raw_n:>3} of {m}")
print(f"  Bonferroni   {bonf:>3}   (flat bar {ALPHA/m:.4e})")
print(f"  Holm         {holm:>3}   ({holm-bonf:+d} vs Bonferroni)")
print(f"  BH (FDR)     {bh:>3}   ({bh-holm:+d} vs Holm)")
k = holm            # the first test Holm rejects
print(f"\n  where Holm stopped: rank {k+1}, p={ps[k]:.4e} vs its bar {ALPHA/(m-k):.4e}"
      f"  (missed by {ps[k]-ALPHA/(m-k):.1e})")