多重検定問題とBonferroni・FDR補正

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# シミュレーション設定
m = 1000  # 検定の数
m0 = 900  # 真の帰無仮説の数
m1 = m - m0  # 真に有意な検定の数
alpha = 0.05

# p値の生成
# H0が真: p ~ Uniform(0,1)
p_null = np.random.uniform(0, 1, m0)
# H0が偽: 効果ありのp値（小さい値に集中）
effect_size = 3.0
z_alt = np.random.normal(effect_size, 1, m1)
p_alt = 2 * (1 - stats.norm.cdf(np.abs(z_alt)))

p_values = np.concatenate([p_null, p_alt])
is_null = np.concatenate([np.ones(m0, dtype=bool), np.zeros(m1, dtype=bool)])

# シャッフル
idx = np.random.permutation(m)
p_values = p_values[idx]
is_null = is_null[idx]

# 補正方法の適用
# (1) 補正なし
reject_none = p_values < alpha

# (2) Bonferroni
reject_bonf = p_values < alpha / m

# (3) Holm
sorted_idx = np.argsort(p_values)
sorted_p = p_values[sorted_idx]
holm_thresholds = alpha / (m - np.arange(m))
reject_holm_sorted = np.zeros(m, dtype=bool)
for i in range(m):
    if sorted_p[i] <= holm_thresholds[i]:
        reject_holm_sorted[i] = True
    else:
        break
reject_holm = np.zeros(m, dtype=bool)
reject_holm[sorted_idx[reject_holm_sorted]] = True

# (4) Benjamini-Hochberg
bh_thresholds = np.arange(1, m+1) * alpha / m
reject_bh_sorted = sorted_p <= bh_thresholds
k_star = np.max(np.where(reject_bh_sorted)[0]) + 1 if np.any(reject_bh_sorted) else 0
reject_bh = np.zeros(m, dtype=bool)
reject_bh[sorted_idx[:k_star]] = True

# 結果の集計
def summarize(reject, is_null):
    V = np.sum(reject & is_null)   # 偽陽性
    S = np.sum(reject & ~is_null)  # 真陽性
    R = np.sum(reject)
    FDR = V / max(R, 1)
    Power = S / max(np.sum(~is_null), 1)
    return V, S, R, FDR, Power

results = {}
for name, reject in [('No correction', reject_none), ('Bonferroni', reject_bonf),
                      ('Holm', reject_holm), ('BH (FDR)', reject_bh)]:
    V, S, R, FDR, Power = summarize(reject, is_null)
    results[name] = {'V': V, 'S': S, 'R': R, 'FDR': FDR, 'Power': Power}

fig, axes = plt.subplots(1, 3, figsize=(16, 5.5))

# (a) p値と各閾値
ax = axes[0]
sorted_p_plot = np.sort(p_values)
rank = np.arange(1, m+1)

ax.scatter(rank, sorted_p_plot, s=2, alpha=0.3, color='gray')
ax.plot(rank, np.full(m, alpha), 'k--', linewidth=1.5, label=f'No correction ($\\alpha={alpha}$)')
ax.plot(rank, np.full(m, alpha/m), 'r-', linewidth=1.5, label=f'Bonferroni ($\\alpha/m$)')
ax.plot(rank, rank * alpha / m, 'b-', linewidth=2, label=f'BH ($i\\alpha/m$)')

# BHの棄却点を表示
if k_star > 0:
    ax.plot(k_star, sorted_p_plot[k_star-1], 'b*', markersize=15, zorder=5)
    ax.annotate(f'$k^* = {k_star}$', xy=(k_star, sorted_p_plot[k_star-1]),
                xytext=(k_star+50, sorted_p_plot[k_star-1]+0.05),
                fontsize=10, arrowprops=dict(arrowstyle='->', color='blue'))

ax.set_xlabel('Rank $i$', fontsize=12)
ax.set_ylabel('$p_{(i)}$', fontsize=12)
ax.set_title('Sorted p-values and Thresholds', fontsize=13)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 200)  # 最初の200個に拡大
ax.set_ylim(0, 0.15)

# (b) 結果の比較
ax = axes[1]
methods = list(results.keys())
x_pos = np.arange(len(methods))
width = 0.25

tp = [results[m]['S'] for m in methods]
fp = [results[m]['V'] for m in methods]
fn = [np.sum(~is_null) - results[m]['S'] for m in methods]

ax.bar(x_pos - width, tp, width, color='green', alpha=0.7, label='True Positive')
ax.bar(x_pos, fp, width, color='red', alpha=0.7, label='False Positive')
ax.bar(x_pos + width, fn, width, color='orange', alpha=0.7, label='False Negative')

for j, m_name in enumerate(methods):
    fdr = results[m_name]['FDR']
    ax.text(j, max(tp[j], fp[j], fn[j]) + 5, f'FDR={fdr:.3f}',
            ha='center', fontsize=8, fontweight='bold')

ax.set_ylabel('Count', fontsize=12)
ax.set_title('Detection Results by Method', fontsize=13)
ax.set_xticks(x_pos)
ax.set_xticklabels(methods, fontsize=8)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis='y')

# (c) FDR vs Power のトレードオフ
ax = axes[2]
alpha_range = np.linspace(0.001, 0.3, 50)
fdr_bonf = []
power_bonf = []
fdr_bh = []
power_bh = []

for a in alpha_range:
    # Bonferroni
    rej = p_values < a / m
    V, S, R, fdr, pwr = summarize(rej, is_null)
    fdr_bonf.append(fdr)
    power_bonf.append(pwr)

    # BH
    sorted_idx_a = np.argsort(p_values)
    sorted_p_a = p_values[sorted_idx_a]
    bh_thresh = np.arange(1, m+1) * a / m
    rej_sorted = sorted_p_a <= bh_thresh
    k = np.max(np.where(rej_sorted)[0]) + 1 if np.any(rej_sorted) else 0
    rej_bh = np.zeros(m, dtype=bool)
    rej_bh[sorted_idx_a[:k]] = True
    V, S, R, fdr, pwr = summarize(rej_bh, is_null)
    fdr_bh.append(fdr)
    power_bh.append(pwr)

ax.plot(fdr_bonf, power_bonf, 'r-', linewidth=2, label='Bonferroni')
ax.plot(fdr_bh, power_bh, 'b-', linewidth=2, label='BH (FDR)')
ax.plot([0, 1], [0, 1], 'k:', alpha=0.3)
ax.set_xlabel('False Discovery Rate', fontsize=12)
ax.set_ylabel('Power (True Positive Rate)', fontsize=12)
ax.set_title('FDR vs Power Trade-off', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 0.3)

plt.tight_layout()
plt.savefig('multiple_testing.png', dpi=150, bbox_inches='tight')
plt.show()