거짓발견율 (Benjamini-Hochberg)¶

import numpy as np

def benjamini_hochberg(p_values, alpha=0.05):
    """벤자미니-호크버그 절차를 적용한다.

        p-값을 작은 것부터 늘어놓고 k번째를 k/m*alpha 와 견준다. 조건을 만족하는
        가장 큰 k 를 찾아 그 아래를 모두 기각한다. FWER 대신 FDR 을 통제하므로
        본페로니보다 훨씬 덜 보수적이다.

    Parameters
    ----------
    p_values : array-like
        Raw p-values from m hypothesis tests.
    alpha : float
        Target FDR level.

    Returns
    -------
    rejected : ndarray of bool
        True for hypotheses that are rejected.
    adjusted : ndarray of float
        BH-adjusted p-values.
    """
    p = np.asarray(p_values)
    m = len(p)
    order = np.argsort(p)
    sorted_p = p[order]

    # BH 문턱값: k번째로 작은 p-값을 k/m*alpha 와 견준다.
    thresholds = np.arange(1, m + 1) / m * alpha

    # 조건을 만족하는 가장 큰 k 를 찾는다. 그 아래는 모두 기각한다.
    below = sorted_p <= thresholds
    if not below.any():
        k = 0
    else:
        k = np.max(np.where(below)[0]) + 1

    # 기각 여부를 표시한다.
    rejected = np.zeros(m, dtype=bool)
    rejected[order[:k]] = True

    # 보정 p-값. 큰 것에서 작은 쪽으로 훑으며 단조성을 강제한다.
    # 이 되짚기가 없으면 원래 p-값의 순서가 뒤집히는 일이 생겨
    # "더 작은 p-값이 더 큰 보정 p-값을 갖는" 이상한 결과가 나온다.
    adjusted_sorted = np.minimum(1, sorted_p * m / np.arange(1, m + 1))
    for i in range(m - 2, -1, -1):
        adjusted_sorted[i] = min(adjusted_sorted[i], adjusted_sorted[i + 1])
    adjusted = np.empty(m)
    adjusted[order] = adjusted_sorted      # 원래 순서로 되돌린다

    return rejected, adjusted


# 아래 연습문제 1의 p-값으로 확인한다.
pvals = np.array([0.001, 0.008, 0.039, 0.041, 0.23, 0.76])
rej, adj = benjamini_hochberg(pvals, alpha=0.10)
print("rejected:", rej)
print("adjusted:", np.round(adj, 4))

# statsmodels와 대조
from statsmodels.stats.multitest import multipletests
rej_sm, adj_sm, _, _ = multipletests(pvals, alpha=0.10, method="fdr_bh")
print("statsmodels adjusted:", np.round(adj_sm, 4))

# (1) 되짚기 전후를 나란히 본다.
m = len(pvals)
raw = np.minimum(1, np.sort(pvals) * m / np.arange(1, m + 1))
print("\n되짚기 전:", np.round(raw, 4))
print("되짚기 후:", np.round(np.sort(adj), 4))
print("단조성이 깨진 자리:",
      [int(k) + 1 for k in np.where(np.diff(raw) < 0)[0]], "번째")

# (2) 조건을 만족하는 k 가 하나도 없는 경우.
none_p = np.array([0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
rej0, adj0 = benjamini_hochberg(none_p, alpha=0.10)
print("\n아무것도 기각되지 않는 경우:", rej0, np.round(adj0, 4))
print("statsmodels 도 같은가:",
      np.array_equal(rej0, multipletests(none_p, alpha=0.10,
                                         method="fdr_bh")[0]))

rejected: [ True  True  True  True False False]
adjusted: [0.006  0.024  0.0615 0.0615 0.276  0.76  ]
statsmodels adjusted: [0.006  0.024  0.0615 0.0615 0.276  0.76  ]

되짚기 전: [0.006  0.024  0.078  0.0615 0.276  0.76  ]
되짚기 후: [0.006  0.024  0.0615 0.0615 0.276  0.76  ]
단조성이 깨진 자리: [3] 번째

아무것도 기각되지 않는 경우: [False False False False False False] [0.9 0.9 0.9 0.9 0.9 0.9]
statsmodels 도 같은가: True

import numpy as np
from scipy import stats

rng = np.random.default_rng(31)
M, m = 400, 1_000

def bh_reject(p, q):
    m = len(p)
    o = np.argsort(p)
    ps = p[o]
    ok = np.where(ps <= q * np.arange(1, m + 1) / m)[0]
    r = np.zeros(m, bool)
    if len(ok):
        r[o[:ok.max() + 1]] = True
    return r

print(f"{'π1':>6s} {'실제 FDR':>10s} {'q·π0':>8s} {'검정력':>8s} "
      f"{'평균 기각 수':>12s}")
for pi1 in [0.01, 0.05, 0.10, 0.30]:
    fdrs, pows, nrej = [], [], []
    for _ in range(M):
        m1 = int(m * pi1)
        z = rng.standard_normal(m)
        z[:m1] += 3.5                        # 참 대립가설
        p = 2 * stats.norm.sf(np.abs(z))
        r = bh_reject(p, 0.05)
        R, V = r.sum(), r[m1:].sum()
        fdrs.append(V / max(R, 1))
        pows.append(r[:m1].sum() / m1)
        nrej.append(R)
    print(f"{pi1:6.2f} {np.mean(fdrs):10.4f} {0.05 * (1 - pi1):8.4f} "
          f"{np.mean(pows):8.4f} {np.mean(nrej):12.1f}")

    π1     실제 FDR     q·π0      검정력      평균 기각 수
  0.01     0.0501   0.0495   0.4325          4.6
  0.05     0.0486   0.0475   0.6409         33.7
  0.10     0.0476   0.0450   0.7320         76.9
  0.30     0.0341   0.0350   0.8449        262.5

def pi0_storey(p, lam=0.5):
    return min((p > lam).sum() / ((1 - lam) * len(p)), 1.0)

m1 = 300
z = rng.standard_normal(m)
z[:m1] += 3.5
p = 2 * stats.norm.sf(np.abs(z))
print(f"참 π0 = {1 - m1 / m:.3f},  스토리 추정 = {pi0_storey(p):.3f}")

참 π0 = 0.700,  스토리 추정 = 0.660

import numpy as np
from scipy import stats

rng = np.random.default_rng(31)
M, m = 400, 1_000

def bh_reject(p, q):
    mm = len(p)
    o = np.argsort(p)
    ok = np.where(p[o] <= q * np.arange(1, mm + 1) / mm)[0]
    r = np.zeros(mm, bool)
    if len(ok):
        r[o[:ok.max() + 1]] = True
    return r

print(f"{'π1':>6s} {'본페로니 검정력':>15s} {'BH 검정력':>11s} {'비':>6s}")
for pi1 in [0.01, 0.05, 0.10]:
    pb, pbh = [], []
    for _ in range(M):
        m1 = int(m * pi1)
        z = rng.standard_normal(m)
        z[:m1] += 3.5
        p = 2 * stats.norm.sf(np.abs(z))
        pb.append((p[:m1] < 0.05 / m).sum() / m1)
        pbh.append(bh_reject(p, 0.05)[:m1].sum() / m1)
    print(f"{pi1:6.2f} {np.mean(pb):15.4f} {np.mean(pbh):11.4f} "
          f"{np.mean(pbh) / np.mean(pb):6.2f}")

    π1    본페로니 검정력    BH 검정력      비
  0.01          0.2925      0.4455    1.52
  0.05          0.2868      0.6419    2.24
  0.10          0.2858      0.7252    2.54

import numpy as np
from scipy import stats

rng = np.random.default_rng(5)
m, m1 = 2_000, 200
z = rng.standard_normal(m)
z[:m1] += 3.0
p = 2 * stats.norm.sf(np.abs(z))

def pi0_storey(p, lam=0.5):
    return min((p > lam).sum() / ((1 - lam) * len(p)), 1.0)

def qvalues(p, pi0=1.0):
    mm = len(p)
    o = np.argsort(p)
    ps = p[o]
    q = pi0 * mm * ps / np.arange(1, mm + 1)
    q = np.minimum.accumulate(q[::-1])[::-1]
    out = np.empty(mm)
    out[o] = np.minimum(q, 1)
    return out

pi0 = pi0_storey(p)
q_cons = qvalues(p)                 # π0 = 1 (BH 와 동일)
q_ada = qvalues(p, pi0)             # π0 추정 반영
print(f"참 π0 = {1 - m1 / m:.3f},  스토리 추정 π0 = {pi0:.3f}")
for thr in [0.01, 0.05, 0.10]:
    print(f"q ≤ {thr:.2f}:  BH {(q_cons <= thr).sum():4d}개   "
          f"적응형 {(q_ada <= thr).sum():4d}개   "
          f"(참 신호 중 발견: {(q_ada[:m1] <= thr).sum():3d}/{m1})")

참 π0 = 0.900,  스토리 추정 π0 = 0.916
q ≤ 0.01:  BH   43개   적응형   43개   (참 신호 중 발견:  43/200)
q ≤ 0.05:  BH   83개   적응형   90개   (참 신호 중 발견:  87/200)
q ≤ 0.10:  BH  115개   적응형  122개   (참 신호 중 발견: 111/200)

import numpy as np
from scipy import stats

rng = np.random.default_rng(77)
M, m, m1 = 300, 500, 50

def bh_reject(p, q):
    mm = len(p)
    o = np.argsort(p)
    ok = np.where(p[o] <= q * np.arange(1, mm + 1) / mm)[0]
    r = np.zeros(mm, bool)
    if len(ok):
        r[o[:ok.max() + 1]] = True
    return r

print(f"{'상관 구조':>14s} {'BH FDR':>9s} {'BY FDR':>9s} "
      f"{'BH 검정력':>11s} {'BY 검정력':>11s}")
c_m = np.sum(1 / np.arange(1, m + 1))
for name, rho in [("독립", 0.0), ("양의 등상관 0.5", 0.5),
                  ("양의 등상관 0.9", 0.9)]:
    f1, f2, p1, p2 = [], [], [], []
    for _ in range(M):
        common = rng.standard_normal(1)
        z = np.sqrt(rho) * common + np.sqrt(1 - rho) * rng.standard_normal(m)
        z[:m1] += 3.5
        p = 2 * stats.norm.sf(np.abs(z))
        for q, fl, pl in [(0.05, f1, p1), (0.05 / c_m, f2, p2)]:
            r = bh_reject(p, q)
            fl.append(r[m1:].sum() / max(r.sum(), 1))
            pl.append(r[:m1].sum() / m1)
    print(f"{name:>14s} {np.mean(f1):9.4f} {np.mean(f2):9.4f} "
          f"{np.mean(p1):11.4f} {np.mean(p2):11.4f}")

         상관 구조    BH FDR    BY FDR      BH 검정력      BY 검정력
            독립    0.0475    0.0075      0.7312      0.4803
    양의 등상관 0.5    0.0396    0.0054      0.6847      0.4330
    양의 등상관 0.9    0.0185    0.0016      0.6963      0.4837

import numpy as np
from scipy import stats

rng = np.random.default_rng(13)
M, m, m1 = 2_000, 500, 25

def bh_reject(p, q):
    mm = len(p)
    o = np.argsort(p)
    ok = np.where(p[o] <= q * np.arange(1, mm + 1) / mm)[0]
    r = np.zeros(mm, bool)
    if len(ok):
        r[o[:ok.max() + 1]] = True
    return r

fdps = []
for _ in range(M):
    z = rng.standard_normal(m)
    z[:m1] += 3.0
    p = 2 * stats.norm.sf(np.abs(z))
    r = bh_reject(p, 0.05)
    fdps.append(r[m1:].sum() / max(r.sum(), 1))
fdps = np.array(fdps)
print(f"FDR(평균) = {fdps.mean():.4f}")
print(f"백분위 50/75/90/95/99: "
      + " ".join(f"{v:.3f}" for v in np.percentile(fdps, [50, 75, 90, 95, 99])))
print(f"실제 비율이 0.10 을 넘은 비율 = {np.mean(fdps > 0.10):.4f}")

FDR(평균) = 0.0478
백분위 50/75/90/95/99: 0.000 0.091 0.143 0.167 0.250
실제 비율이 0.10 을 넘은 비율 = 0.1910