거짓발견율과 재표본추출 다중검정¶

import numpy as np

# 검정 m 개가 모두 참 귀무가설이고 서로 독립이라면, 하나도 잘못 기각하지
# 않을 확률이 (1-a)^m 이다. 적어도 하나를 잘못 기각할 확률이 그 나머지다.
# m=100, a=0.05 면 0.994 — 거의 확실하게 거짓 양성이 하나는 나온다.
m_vals = np.arange(1, 501)
alphas = [0.05, 0.01, 0.001]
for a in alphas:
    fwer = 1 - (1 - a) ** m_vals
    print(f"alpha={a}, m=100: FWER={1 - (1-a)**100:.4f}")

# (2) FWER 을 0.05 로 맞추는 alpha.
m = 100
sidak = 1 - 0.95 ** (1 / m)
bonf = 0.05 / m
print(f"\n시닥   alpha = {sidak:.8f}   확인 FWER = {1 - (1 - sidak) ** m:.6f}")
print(f"본페로니 alpha = {bonf:.8f}   확인 FWER = {1 - (1 - bonf) ** m:.6f}")
print(f"비 시닥/본페로니 = {sidak / bonf:.4f}")

alpha=0.05, m=100: FWER=0.9941
alpha=0.01, m=100: FWER=0.6340
alpha=0.001, m=100: FWER=0.0952

시닥   alpha = 0.00051280   확인 FWER = 0.050000
본페로니 alpha = 0.00050000   확인 FWER = 0.048782
비 시닥/본페로니 = 1.0256

from scipy import stats
from statsmodels.stats.multitest import multipletests

np.random.seed(42)

n_tests = 2000
n_true_alt = 200
n_obs = 50
effect_size = 0.5

p_values = np.zeros(n_tests)
truth = np.zeros(n_tests, dtype=int)
truth[:n_true_alt] = 1      # 앞의 200개만 참 신호, 나머지 1800개는 귀무

for i in range(n_tests):
    # 모의실험이라 정답을 알고 있다. 그래서 TP와 FP를 직접 셀 수 있다.
    # 실제 자료에서는 이 정보가 없으므로 FDR을 추정해야 한다.
    mu = effect_size if i < n_true_alt else 0.0
    data = np.random.normal(mu, 1.0, n_obs)
    _, p_values[i] = stats.ttest_1samp(data, 0)

# 세 보정 방법을 같은 p-값 묶음에 적용한다.
_, p_bonf, _, _ = multipletests(p_values, method="bonferroni")
_, p_holm, _, _ = multipletests(p_values, method="holm")
_, p_bh, _, _   = multipletests(p_values, method="fdr_bh")

alpha = 0.05
for name, adj_p in [("Bonferroni", p_bonf), ("Holm", p_holm), ("BH", p_bh)]:
    rejected = adj_p < alpha
    tp = np.sum(rejected & (truth == 1))
    fp = np.sum(rejected & (truth == 0))
    fdr = fp / max(np.sum(rejected), 1)
    power = tp / np.sum(truth == 1)
    print(f"{name:12s}: TP={tp}, FP={fp}, FDR={fdr:.3f}, Power={power:.3f}")

def resampling_fdr(X_group1, X_group2, n_permutations=500):
    n1, n2 = X_group1.shape[0], X_group2.shape[0]
    n_features = X_group1.shape[1]
    X_combined = np.vstack([X_group1, X_group2])

    # 관측된 검정통계량
    t_obs = np.array([
        stats.ttest_ind(X_group1[:, j], X_group2[:, j]).statistic
        for j in range(n_features)
    ])

    # 집단 이름을 뒤섞어 만든 귀무분포. 모형을 가정하지 않는다.
    t_perm = np.zeros((n_permutations, n_features))
    for b in range(n_permutations):
        idx = np.random.permutation(n1 + n2)
        for j in range(n_features):
            t_perm[b, j] = stats.ttest_ind(
                X_combined[idx[:n1], j],
                X_combined[idx[n1:], j]
            ).statistic

    # 문턱마다 FDR을 추정한다.
    Rs, FDRs = [], []
    for thresh in np.sort(np.abs(t_obs)):
        R = np.sum(np.abs(t_obs) >= thresh)       # 실제 기각 수
        # 순열 자료에서 문턱을 넘은 총 개수를 순열 횟수로 나눈다.
        # 이것이 "H0가 참일 때 기대되는 기각 수", 즉 E[V]의 추정이다.
        V = np.sum(np.abs(t_perm) >= thresh) / n_permutations
        Rs.append(R)
        FDRs.append(V / max(R, 1))
    return np.array(Rs), np.array(FDRs)

np.random.seed(0)

m, m1, n1, n2, delta = 200, 40, 25, 25, 1.0
m0 = m - m1
truth = np.zeros(m, dtype=int)
truth[:m1] = 1                      # 앞의 40개 특징만 참 신호

X1 = np.random.normal(0.0, 1.0, (n1, m))
X2 = np.random.normal(0.0, 1.0, (n2, m))
X2[:, :m1] += delta                 # 신호는 집단 2 의 평균을 delta 만큼 옮긴다

t_obs = np.array([stats.ttest_ind(X1[:, j], X2[:, j]).statistic
                  for j in range(m)])
Rs, FDRs = resampling_fdr(X1, X2, n_permutations=200)

df = n1 + n2 - 2
print(f"m = {m},  m0 = {m0},  m/m0 = {m / m0:.4f},  자유도 = {df}")

print("\n순열 귀무분포가 t_48 과 맞는가")
print("    t     순열 P(|T|>=t)    2*sf(t, 48)")
# resampling_fdr 안에서 쓴 것과 같은 순열을 다시 만들 수는 없으므로
# 같은 방식으로 새로 200 번 섞어 비교한다.
Xc = np.vstack([X1, X2])
cnt = np.zeros(4)
grid = np.array([1.0, 2.0, 2.5, 3.0])
B = 200
for b in range(B):
    idx = np.random.permutation(n1 + n2)
    tp = np.array([stats.ttest_ind(Xc[idx[:n1], j], Xc[idx[n1:], j]).statistic
                   for j in range(m)])
    cnt += [(np.abs(tp) >= g).sum() for g in grid]
for g, c in zip(grid, cnt):
    print(f"{g:5.1f}      {c / (B * m):11.6f}    {2 * stats.t.sf(g, df):11.6f}")

print("\n   R    문턱    FDR 추정   m*p(t)/R   m0*p(t)/R   비   실현 FDP")
thr = np.sort(np.abs(t_obs))
for k in range(len(thr)):
    R = Rs[k]
    if R not in (200, 100, 60, 50, 40, 30, 20):
        continue
    t = thr[k]
    p_t = 2 * stats.t.sf(t, df)
    sel = np.abs(t_obs) >= t
    fdp = np.sum(sel & (truth == 0)) / max(sel.sum(), 1)
    print(f"{R:4d} {t:7.4f}  {FDRs[k]:9.4f}  {m * p_t / R:9.4f}"
          f"   {m0 * p_t / R:9.4f}  {m * p_t / (m0 * p_t):5.2f}  {fdp:9.4f}")

m = 200,  m0 = 160,  m/m0 = 1.2500,  자유도 = 48

순열 귀무분포가 t_48 과 맞는가
    t     순열 P(|T|>=t)    2*sf(t, 48)
  1.0         0.322125       0.322325
  2.0         0.052850       0.051176
  2.5         0.016450       0.015890
  3.0         0.004375       0.004272

   R    문턱    FDR 추정   m*p(t)/R   m0*p(t)/R   비   실현 FDP
 200  0.0249     0.9807     0.9802      0.7842   1.25     0.8000
 100  0.9308     0.7107     0.7133      0.5706   1.25     0.6000
  60  1.5739     0.3944     0.4069      0.3255   1.25     0.3500
  50  1.8511     0.2711     0.2813      0.2250   1.25     0.2200
  40  2.1825     0.1635     0.1700      0.1360   1.25     0.0500
  30  2.9122     0.0350     0.0362      0.0290   1.25     0.0000
  20  3.5186     0.0065     0.0096      0.0077   1.25     0.0000

for alpha in [0.05, 0.10]:
    _, p_bh, _, _ = multipletests(p_values, alpha=alpha, method="fdr_bh")
    rejected = p_bh < alpha
    tp = np.sum(rejected & (truth == 1))
    fp = np.sum(rejected & (truth == 0))
    fdr = fp / max(np.sum(rejected), 1)
    power = tp / np.sum(truth == 1)
    print(f"alpha={alpha}: Rejections={np.sum(rejected)}, "
          f"FDR={fdr:.3f}, Power={power:.3f}")

alpha=0.05: Rejections=142, FDR=0.063, Power=0.665
alpha=0.1: Rejections=179, FDR=0.128, Power=0.780

import numpy as np
from scipy import stats

rng = np.random.default_rng(202)
n, m, m1, B = 40, 500, 50, 500

rho = 0.6                                   # 결과변수 간 상관
S = rho * np.ones((m, m)) + (1 - rho) * np.eye(m)
L = np.linalg.cholesky(S)
g = np.repeat([0, 1], n // 2)
y = rng.standard_normal((n, m)) @ L.T
y[g == 1, :m1] += 1.0

def tstat(y, g):
    a, b = y[g == 0], y[g == 1]
    se = np.sqrt(a.var(0, ddof=1) / len(a) + b.var(0, ddof=1) / len(b))
    return np.abs(b.mean(0) - a.mean(0)) / se

t_obs = tstat(y, g)
null = np.empty((B, m))
for b in range(B):
    null[b] = tstat(y, rng.permutation(g))

grid = np.linspace(1.0, 6.0, 200)
best = None
for c in grid:
    R = (t_obs > c).sum()
    if R == 0:
        continue
    V = (null > c).sum() / B
    if V / R <= 0.05:
        best = (c, R, V / R)
        break
c, R, fdr_hat = best
true_V = (t_obs[m1:] > c).sum()
print(f"순열 FDR: 문턱 {c:.3f}, 기각 {R}개, 추정 FDR {fdr_hat:.4f}, "
      f"실제 거짓 {true_V}개 (실제 비율 {true_V / R:.4f})")

# 비교: 이론적 p-값 + BH
p = 2 * stats.t.sf(t_obs, n - 2)
o = np.argsort(p)
ok = np.where(p[o] <= 0.05 * np.arange(1, m + 1) / m)[0]
R2 = ok.max() + 1 if len(ok) else 0
V2 = (o[:R2] >= m1).sum()
print(f"이론 p + BH: 기각 {R2}개, 실제 거짓 {V2}개 "
      f"(실제 비율 {V2 / max(R2, 1):.4f})")

순열 FDR: 문턱 3.337, 기각 23개, 추정 FDR 0.0483, 실제 거짓 0개 (실제 비율 0.0000)
이론 p + BH: 기각 26개, 실제 거짓 0개 (실제 비율 0.0000)

import numpy as np
from scipy import stats

pi0 = 0.90
mu1, s1 = 3.0, 1.0                 # 대립분포 N(3,1)

def f0(z):
    return stats.norm.pdf(z)

def f1(z):
    return stats.norm.pdf(z, mu1, s1)

def fdr_local(z):
    return pi0 * f0(z) / (pi0 * f0(z) + (1 - pi0) * f1(z))

def fdr_tail(z):
    return (pi0 * stats.norm.sf(z)
            / (pi0 * stats.norm.sf(z) + (1 - pi0) * stats.norm.sf(z, mu1, s1)))

print(f"{'z':>5s} {'국소 fdr':>10s} {'꼬리 Fdr':>10s}")
for z in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0]:
    print(f"{z:5.1f} {fdr_local(z):10.4f} {fdr_tail(z):10.4f}")

    z    국소 fdr    꼬리 Fdr
  1.5     0.8438     0.7239
  2.0     0.7326     0.5688
  2.5     0.5786     0.3999
  3.0     0.4000     0.2519
  3.5     0.2338     0.1418
  4.0     0.1177     0.0714
  5.0     0.0184     0.0116

import numpy as np

print(f"{'B':>8s} {'p=0.05 의 SE':>13s} {'p=0.01 의 SE':>13s} "
      f"{'최소 p':>10s}")
for B in [99, 999, 9_999, 99_999]:
    print(f"{B:8d} {np.sqrt(0.05 * 0.95 / B):13.5f} "
          f"{np.sqrt(0.01 * 0.99 / B):13.5f} {1 / (B + 1):10.2e}")

       B   p=0.05 의 SE   p=0.01 의 SE       최소 p
      99       0.02190       0.00998   1.00e-02
     999       0.00689       0.00315   1.00e-03
    9999       0.00218       0.00099   1.00e-04
   99999       0.00069       0.00031   1.00e-05

import numpy as np
from scipy import stats

rng = np.random.default_rng(404)
m, m1 = 5_000, 250

for sigma0 in [1.0, 1.2, 1.5]:
    z = rng.normal(0, sigma0, m)          # 귀무 통계량이 퍼져 있음
    z[:m1] = rng.normal(3.5, 1.0, m1)     # 참 신호
    p = 2 * stats.norm.sf(np.abs(z))      # N(0,1) 을 가정한 p-값

    o = np.argsort(p)
    ok = np.where(p[o] <= 0.05 * np.arange(1, m + 1) / m)[0]
    R = ok.max() + 1 if len(ok) else 0
    V = (o[:R] >= m1).sum()
    print(f"σ0 = {sigma0}: BH 기각 {R:4d}개, 실제 거짓 {V:4d}개, "
          f"실제 FDR {V / max(R, 1):.4f}")

σ0 = 1.0: BH 기각  163개, 실제 거짓    8개, 실제 FDR 0.0491
σ0 = 1.2: BH 기각  225개, 실제 거짓   46개, 실제 FDR 0.2044
σ0 = 1.5: BH 기각  471개, 실제 거짓  289개, 실제 FDR 0.6136

from scipy.optimize import curve_fit

rng = np.random.default_rng(404)
z = rng.normal(0, 1.2, m)
z[:m1] = rng.normal(3.5, 1.0, m1)

mid = z[np.abs(z) < 2]                   # 중앙 부분
print(f"경험적 귀무: 중심 {mid.mean():.3f}, 척도 "
      f"{mid.std(ddof=1) / np.sqrt(1 - 2 * 2 * stats.norm.pdf(2)
                                   / (2 * stats.norm.cdf(2) - 1)):.3f}")

경험적 귀무: 중심 0.012, 척도 1.104