統計的仮説検定を1回だけ行う場合、第1種の過誤の確率は有意水準 $\alpha$ で制御されています。しかし、同じデータに対して複数の検定を同時に行うと、少なくとも1つの検定で偽陽性（第1種の過誤）が生じる確率が $\alpha$ を大幅に超えてしまいます。これが 多重比較問題（multiple comparison problem） です。
ゲノム解析では数万の遺伝子を同時に検定し、脳画像解析では数万のボクセルを同時に検定し、A/Bテストでは複数の指標を同時に比較します。これらの場面で多重比較問題を無視すると、大量の偽陽性が生じ、誤った結論に至ります。
本記事の内容
- 多重比較問題とFWER（族ごとの過誤率）の増大の数学的導出
- ボンフェローニ補正の理論と導出
- Holm法（ステップダウン法）のアルゴリズムと理論
- FDR（False Discovery Rate）の概念
- Benjamini-Hochberg法の導出と手順
- Pythonでのシミュレーションと各補正法の比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
多重比較問題とは
1回の検定の場合
有意水準 $\alpha = 0.05$ で1回の検定を行うとき、帰無仮説が真であれば第1種の過誤を犯す確率は $\alpha = 0.05$ です。つまり、正しく帰無仮説を棄却しない確率は $1 – \alpha = 0.95$ です。
複数回の検定の場合
$m$ 個の独立な検定を同時に行い、すべての帰無仮説が真であるとしましょう。各検定で帰無仮説を棄却しない確率は $1 – \alpha$ ですから、すべての検定で帰無仮説を棄却しない確率は次のようになります。
$$ P(\text{すべて正しく棄却しない}) = (1 – \alpha)^m $$
ここで各検定の独立性を仮定しました。
したがって、少なくとも1つの検定で誤って棄却してしまう確率（族ごとの過誤率: FWER）は次のようになります。
$$ \text{FWER} = 1 – (1 – \alpha)^m $$
FWERの急激な増大
具体的な値を見てみましょう。
| 検定数 $m$ | FWER（$\alpha = 0.05$） |
|---|---|
| 1 | 0.050 |
| 5 | 0.226 |
| 10 | 0.401 |
| 20 | 0.642 |
| 50 | 0.923 |
| 100 | 0.994 |
$m = 20$ の場合、すべての帰無仮説が真であっても、64%以上の確率で少なくとも1つは「有意」と判定されてしまいます。$m = 100$ では、ほぼ確実に偽陽性が生じます。
FWERの近似
$\alpha$ が小さく $m$ がそれほど大きくない場合、FWER は次のように近似できます。
$$ \text{FWER} = 1 – (1 – \alpha)^m \approx m\alpha $$
この近似は、$(1-\alpha)^m \approx 1 – m\alpha$（二項近似、$m\alpha \ll 1$ のとき）に基づいています。より厳密には、ブールの不等式（union bound）を用いると次が成り立ちます。
$$ \text{FWER} = P\left(\bigcup_{i=1}^{m} A_i\right) \leq \sum_{i=1}^{m} P(A_i) = m\alpha $$
ここで $A_i$ は「$i$ 番目の検定で第1種の過誤を犯す」事象です。この不等式は、検定間の独立性を仮定せずに成り立つ点が重要です。
ボンフェローニ補正
基本的なアイデア
FWERを $\alpha$ 以下に制御するための最もシンプルな方法が ボンフェローニ補正（Bonferroni correction） です。
アイデアは単純です。$m$ 個の検定を行うとき、各検定の有意水準を $\alpha/m$ に引き下げます。
$$ \text{補正後の有意水準}: \alpha^* = \frac{\alpha}{m} $$
数学的導出
ボンフェローニ補正がFWERを $\alpha$ 以下に制御することを証明します。
$A_i$ を「$i$ 番目の検定で第1種の過誤を犯す」事象とします。FWERは次のように表されます。
$$ \text{FWER} = P\left(\bigcup_{i=1}^{m} A_i\right) $$
ブールの不等式（確率の劣加法性）を適用すると、
$$ P\left(\bigcup_{i=1}^{m} A_i\right) \leq \sum_{i=1}^{m} P(A_i) $$
各検定の有意水準を $\alpha/m$ に設定すれば $P(A_i) \leq \alpha/m$ ですから、
$$ \text{FWER} \leq \sum_{i=1}^{m} \frac{\alpha}{m} = m \cdot \frac{\alpha}{m} = \alpha $$
これでFWERが $\alpha$ 以下に制御されることが示されました。
等価な表現: p値の補正
実際の適用では、有意水準を下げる代わりに、p値を $m$ 倍して通常の有意水準 $\alpha$ と比較する方法も等価です。
$$ p_i^{\text{adj}} = \min(m \cdot p_i, 1) $$
$p_i^{\text{adj}} \leq \alpha$ のとき、$i$ 番目の帰無仮説を棄却します。
ボンフェローニ補正の限界
ボンフェローニ補正は保守的です。つまり、FWERを確実に制御する代わりに、検出力が大きく低下します。
有意水準が $\alpha/m$ に下がるため、検出すべき真の効果も見逃しやすくなります（第2種の過誤が増える）。$m$ が大きいほどこの問題は深刻になります。
Holm法（ステップダウン法）
動機
ボンフェローニ補正の保守性を緩和する方法として、Holm法（Holm-Bonferroni法） が提案されました（Holm, 1979）。Holm法はボンフェローニ補正と同じくFWERを $\alpha$ 以下に制御しますが、検出力が向上します。
アルゴリズム
Holm法のアルゴリズムは以下の通りです。
ステップ1: $m$ 個のp値を小さい順にソートします。
$$ p_{(1)} \leq p_{(2)} \leq \dots \leq p_{(m)} $$
ステップ2: 小さい方から順に、以下の基準で比較します。
- $p_{(1)}$ を $\frac{\alpha}{m}$ と比較
- $p_{(2)}$ を $\frac{\alpha}{m-1}$ と比較
- $p_{(3)}$ を $\frac{\alpha}{m-2}$ と比較
- …一般に、$p_{(k)}$ を $\frac{\alpha}{m-k+1}$ と比較
ステップ3: $p_{(k)} > \frac{\alpha}{m-k+1}$ となる最小の $k$ を見つけます。$p_{(k)}, p_{(k+1)}, \dots, p_{(m)}$ に対応する帰無仮説は棄却せず、$p_{(1)}, \dots, p_{(k-1)}$ に対応する帰無仮説を棄却します。
Holm法がFWERを制御する理由
Holm法がFWERを制御することの証明を示します。
$m_0$ 個の真の帰無仮説があるとします（$m_0 \leq m$）。FWER制御に失敗する（真の帰無仮説を少なくとも1つ棄却する）ためには、ソートされたp値列のどこかで、真の帰無仮説に対応するp値がその段階の閾値 $\alpha/(m – k + 1)$ 以下である必要があります。
真の帰無仮説 $H_{0,i}$ のp値 $p_i$ は $[0, 1]$ 上の一様分布に従います（帰無仮説が真のとき）。Holm法で真の帰無仮説が棄却されるためには、少なくとも次のいずれかが成り立つ必要があります。
$$ p_i \leq \frac{\alpha}{m – k + 1} \leq \frac{\alpha}{m_0} $$
最後の不等式は、$k \leq m – m_0 + 1$（真でない帰無仮説がすべて棄却された後でも、$m – k + 1 \geq m_0$）から従います。
ブールの不等式を適用すると、
$$ \text{FWER} \leq \sum_{i \in \mathcal{H}_0} P\left(p_i \leq \frac{\alpha}{m_0}\right) = m_0 \cdot \frac{\alpha}{m_0} = \alpha $$
ここで $\mathcal{H}_0$ は真の帰無仮説の集合です。
補正p値
Holm法の補正p値は次のように計算されます。
$$ p_{(k)}^{\text{adj}} = \max_{j \leq k}\left\{(m – j + 1) \cdot p_{(j)}\right\} $$
最大値を取るのは、補正p値の単調性（$p_{(1)}^{\text{adj}} \leq p_{(2)}^{\text{adj}} \leq \dots$）を保証するためです。
FDR（False Discovery Rate）
FWERの限界
大規模な検定（$m$ が数千〜数万）では、FWERの制御は非常に保守的になります。ゲノム解析で2万の遺伝子を同時に検定する場合、ボンフェローニ補正では有意水準が $0.05 / 20000 = 2.5 \times 10^{-6}$ となり、真の効果をほとんど検出できなくなります。
FDRの定義
FDR（False Discovery Rate） は、棄却した帰無仮説のうち、真の帰無仮説（偽陽性）の割合の期待値です。
$$ \text{FDR} = E\left[\frac{V}{R}\right] $$
ここで、
- $V$: 偽陽性の数（真の帰無仮説のうち棄却されたもの）
- $R$: 棄却した帰無仮説の総数（$R = 0$ の場合は $V/R = 0$ と定義）
これを表にまとめます。
| 棄却しない | 棄却する | 合計 | |
|---|---|---|---|
| $H_0$ が真 | $U$ | $V$ | $m_0$ |
| $H_0$ が偽 | $T$ | $S$ | $m – m_0$ |
| 合計 | $m – R$ | $R$ | $m$ |
FDRは「発見」（棄却）のうち偽の発見の割合をコントロールします。
FWERとFDRの関係
すべての帰無仮説が真（$m_0 = m$）の場合、$R > 0$ ならば必ず $V = R$ なので $V/R = 1$ となります。このとき、
$$ \text{FDR} = E\left[\frac{V}{R}\right] = P(R > 0) = \text{FWER} $$
一般には、FDRはFWER以下になります。
$$ \text{FDR} \leq \text{FWER} $$
FDRを制御すれば、FWERよりも緩やかな基準になるため、より多くの真の効果を検出できます。
Benjamini-Hochberg法（BH法）
アルゴリズム
Benjamini-Hochberg法（BH法）（Benjamini & Hochberg, 1995）は、FDRを指定レベル $q$ 以下に制御する手続きです。
ステップ1: $m$ 個のp値を小さい順にソートします。
$$ p_{(1)} \leq p_{(2)} \leq \dots \leq p_{(m)} $$
ステップ2: 次の条件を満たす最大の $k$ を見つけます。
$$ p_{(k)} \leq \frac{k}{m} \cdot q $$
ステップ3: $p_{(1)}, p_{(2)}, \dots, p_{(k)}$ に対応する帰無仮説を棄却します。条件を満たす $k$ が存在しなければ、何も棄却しません。
BH法がFDRを制御する証明のスケッチ
BH法がFDR $\leq q$ を保証することの証明の核心を示します（Benjamini & Hochberg, 1995）。
p値が独立で、真の帰無仮説のp値が一様分布 $U(0, 1)$ に従うとします。$m_0$ 個の真の帰無仮説に対して、
$$ \text{FDR} = E\left[\frac{V}{R} \mid R > 0\right] P(R > 0) \leq \frac{m_0}{m} q \leq q $$
直感的には、BH法の閾値 $\frac{k}{m}q$ が $k$ に対して線形に増加するため、偽の帰無仮説（真の効果がある）のp値は小さい値に集中し、真の帰無仮説のp値は一様分布に従うことから、両者をうまく分離できるという仕組みです。
補正p値（adjusted p-value / q-value）
BH法の補正p値は次のように計算されます。
$$ p_{(k)}^{\text{adj}} = \min_{j \geq k}\left\{\frac{m}{j} \cdot p_{(j)}\right\} $$
最小値を取るのは、補正p値の単調性を保証するためです。逆方向（大きい方から小さい方）に処理する点がHolm法と対照的です。
具体的なアルゴリズムとしては、以下の手順で計算します。
$$ p_{(m)}^{\text{adj}} = p_{(m)} $$
$$ p_{(k)}^{\text{adj}} = \min\left(\frac{m}{k} \cdot p_{(k)}, \ p_{(k+1)}^{\text{adj}}\right) \quad (k = m-1, m-2, \dots, 1) $$
具体例
数値例
5つの検定を同時に行い、以下のp値を得たとします。
$$ p_1 = 0.005, \quad p_2 = 0.013, \quad p_3 = 0.031, \quad p_4 = 0.045, \quad p_5 = 0.250 $$
$\alpha = 0.05$, $m = 5$ です。
補正なし: $p_1, p_2, p_3, p_4$ が有意（4つ棄却）
ボンフェローニ補正: 補正有意水準 $\alpha^* = 0.05/5 = 0.01$ – $p_1 = 0.005 < 0.01$ → 棄却 - 他はすべて $> 0.01$ → 棄却しない（1つ棄却）
Holm法: p値をソートして順に比較 – $p_{(1)} = 0.005$ vs $\alpha/5 = 0.010$ → $0.005 < 0.010$ → 棄却 - $p_{(2)} = 0.013$ vs $\alpha/4 = 0.0125$ → $0.013 > 0.0125$ → 棄却しない → ここで停止
結果: $p_{(1)}$ のみ棄却（1つ棄却）
BH法: $q = 0.05$ として、$p_{(k)} \leq \frac{k}{5} \times 0.05$ を満たす最大の $k$ を探す
| $k$ | $p_{(k)}$ | $\frac{k}{5} \times 0.05$ | 比較 |
|---|---|---|---|
| 1 | 0.005 | 0.010 | $0.005 \leq 0.010$ |
| 2 | 0.013 | 0.020 | $0.013 \leq 0.020$ |
| 3 | 0.031 | 0.030 | $0.031 > 0.030$ |
| 4 | 0.045 | 0.040 | $0.045 > 0.040$ |
| 5 | 0.250 | 0.050 | $0.250 > 0.050$ |
条件を満たす最大の $k = 2$ なので、$p_{(1)}, p_{(2)}$ を棄却（2つ棄却）。
Pythonでの実装
各補正法のスクラッチ実装
import numpy as np
def bonferroni_correction(p_values, alpha=0.05):
"""
ボンフェローニ補正
Parameters
----------
p_values : array-like
p値の配列
alpha : float
有意水準
Returns
-------
adjusted_p : ndarray
補正p値
rejected : ndarray
棄却されたかどうか (True/False)
"""
p = np.array(p_values)
m = len(p)
adjusted_p = np.minimum(p * m, 1.0)
rejected = adjusted_p <= alpha
return adjusted_p, rejected
def holm_correction(p_values, alpha=0.05):
"""
Holm法（ステップダウン法）
Parameters
----------
p_values : array-like
p値の配列
alpha : float
有意水準
Returns
-------
adjusted_p : ndarray
補正p値（元の順序）
rejected : ndarray
棄却されたかどうか (True/False)
"""
p = np.array(p_values)
m = len(p)
# ソートされたインデックス
sorted_idx = np.argsort(p)
sorted_p = p[sorted_idx]
# 補正p値の計算
adjusted_sorted = np.zeros(m)
for k in range(m):
adjusted_sorted[k] = (m - k) * sorted_p[k]
# 単調性の保証（累積最大値）
for k in range(1, m):
adjusted_sorted[k] = max(adjusted_sorted[k], adjusted_sorted[k-1])
# 1を超えないようにクリップ
adjusted_sorted = np.minimum(adjusted_sorted, 1.0)
# 元の順序に戻す
adjusted_p = np.zeros(m)
adjusted_p[sorted_idx] = adjusted_sorted
rejected = adjusted_p <= alpha
return adjusted_p, rejected
def benjamini_hochberg(p_values, q=0.05):
"""
Benjamini-Hochberg法（FDR制御）
Parameters
----------
p_values : array-like
p値の配列
q : float
FDRの目標水準
Returns
-------
adjusted_p : ndarray
補正p値（元の順序）
rejected : ndarray
棄却されたかどうか (True/False)
"""
p = np.array(p_values)
m = len(p)
# ソートされたインデックス
sorted_idx = np.argsort(p)
sorted_p = p[sorted_idx]
# 補正p値の計算（逆方向）
adjusted_sorted = np.zeros(m)
adjusted_sorted[m-1] = sorted_p[m-1]
for k in range(m-2, -1, -1):
adjusted_sorted[k] = min(m / (k + 1) * sorted_p[k], adjusted_sorted[k+1])
# 1を超えないようにクリップ
adjusted_sorted = np.minimum(adjusted_sorted, 1.0)
# 元の順序に戻す
adjusted_p = np.zeros(m)
adjusted_p[sorted_idx] = adjusted_sorted
rejected = adjusted_p <= q
return adjusted_p, rejected
# --- 具体例 ---
p_values = [0.005, 0.013, 0.031, 0.045, 0.250]
print("元のp値:", p_values)
print("=" * 60)
# 補正なし
print("\n【補正なし】")
for i, p in enumerate(p_values):
print(f" p{i+1} = {p:.3f} → {'棄却' if p <= 0.05 else '棄却しない'}")
# ボンフェローニ
adj_p, rej = bonferroni_correction(p_values)
print("\n【ボンフェローニ補正】")
for i in range(len(p_values)):
print(f" p{i+1} = {p_values[i]:.3f} → adj_p = {adj_p[i]:.3f} → {'棄却' if rej[i] else '棄却しない'}")
# Holm法
adj_p, rej = holm_correction(p_values)
print("\n【Holm法】")
for i in range(len(p_values)):
print(f" p{i+1} = {p_values[i]:.3f} → adj_p = {adj_p[i]:.3f} → {'棄却' if rej[i] else '棄却しない'}")
# BH法
adj_p, rej = benjamini_hochberg(p_values)
print("\n【Benjamini-Hochberg法】")
for i in range(len(p_values)):
print(f" p{i+1} = {p_values[i]:.3f} → adj_p = {adj_p[i]:.3f} → {'棄却' if rej[i] else '棄却しない'}")
実行結果は次のようになります。
元のp値: [0.005, 0.013, 0.031, 0.045, 0.25]
============================================================
【補正なし】
p1 = 0.005 → 棄却
p2 = 0.013 → 棄却
p3 = 0.031 → 棄却
p4 = 0.045 → 棄却
p5 = 0.250 → 棄却しない
【ボンフェローニ補正】
p1 = 0.005 → adj_p = 0.025 → 棄却
p2 = 0.013 → adj_p = 0.065 → 棄却しない
p3 = 0.031 → adj_p = 0.155 → 棄却しない
p4 = 0.045 → adj_p = 0.225 → 棄却しない
p5 = 0.250 → adj_p = 1.000 → 棄却しない
【Holm法】
p1 = 0.005 → adj_p = 0.025 → 棄却
p2 = 0.013 → adj_p = 0.052 → 棄却しない
p3 = 0.031 → adj_p = 0.093 → 棄却しない
p4 = 0.045 → adj_p = 0.093 → 棄却しない
p5 = 0.250 → adj_p = 0.250 → 棄却しない
【Benjamini-Hochberg法】
p1 = 0.005 → adj_p = 0.025 → 棄却
p2 = 0.013 → adj_p = 0.033 → 棄却
p3 = 0.031 → adj_p = 0.052 → 棄却しない
p4 = 0.045 → adj_p = 0.056 → 棄却しない
p5 = 0.250 → adj_p = 0.250 → 棄却しない
FWERの増大をシミュレーションで確認
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n_simulations = 10000
alpha = 0.05
m_values = [1, 2, 5, 10, 20, 50, 100]
# シミュレーション: すべてのH0が真のケース
fwer_empirical = []
fwer_theoretical = []
for m in m_values:
false_positive_count = 0
for _ in range(n_simulations):
# m個の検定（すべてH0が真 → p値は一様分布）
p_vals = np.random.uniform(0, 1, m)
# 少なくとも1つがα以下なら偽陽性
if np.any(p_vals <= alpha):
false_positive_count += 1
fwer_empirical.append(false_positive_count / n_simulations)
fwer_theoretical.append(1 - (1 - alpha)**m)
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(m_values, fwer_theoretical, 'b-o', linewidth=2, markersize=8,
label='Theoretical: $1 - (1-\\alpha)^m$')
ax.plot(m_values, fwer_empirical, 'r--s', linewidth=2, markersize=8,
label='Simulated FWER')
ax.axhline(y=alpha, color='gray', linestyle=':', linewidth=1.5,
label=f'$\\alpha$ = {alpha}')
ax.set_xlabel('Number of Tests $m$', fontsize=13)
ax.set_ylabel('FWER', fontsize=13)
ax.set_title('Family-Wise Error Rate vs Number of Tests', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)
plt.tight_layout()
plt.show()
各補正法のFWER・FDR・検出力の比較シミュレーション
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
def simulate_multiple_testing(m, m1, effect_size, n_obs, n_sim, alpha=0.05):
"""
多重検定のシミュレーション
Parameters
----------
m : int
検定の総数
m1 : int
対立仮説が真の検定の数
effect_size : float
効果量（Cohen's d）
n_obs : int
各検定の標本サイズ
n_sim : int
シミュレーション回数
alpha : float
有意水準
Returns
-------
results : dict
各補正法ごとのFWER, FDR, Power
"""
m0 = m - m1 # 帰無仮説が真の数
results = {method: {'V': [], 'S': [], 'R': []} for method in ['none', 'bonferroni', 'holm', 'bh']}
for _ in range(n_sim):
p_values = np.zeros(m)
# 帰無仮説が真の検定: p値は一様分布
for i in range(m0):
# 1標本t検定を実際にシミュレーション
sample = np.random.normal(0, 1, n_obs)
_, p_values[i] = stats.ttest_1samp(sample, 0)
# 対立仮説が真の検定: 効果ありのデータ
for i in range(m0, m):
sample = np.random.normal(effect_size, 1, n_obs)
_, p_values[i] = stats.ttest_1samp(sample, 0)
# 各補正法を適用
# 補正なし
rej_none = p_values <= alpha
# ボンフェローニ
adj_bonf, rej_bonf = bonferroni_correction(p_values, alpha)
# Holm法
adj_holm, rej_holm = holm_correction(p_values, alpha)
# BH法
adj_bh, rej_bh = benjamini_hochberg(p_values, alpha)
for method, rej in [('none', rej_none), ('bonferroni', rej_bonf),
('holm', rej_holm), ('bh', rej_bh)]:
V = np.sum(rej[:m0]) # 偽陽性の数
S = np.sum(rej[m0:]) # 真陽性の数
R = np.sum(rej) # 棄却数の合計
results[method]['V'].append(V)
results[method]['S'].append(S)
results[method]['R'].append(R)
# 統計量の計算
summary = {}
for method in results:
V = np.array(results[method]['V'])
S = np.array(results[method]['S'])
R = np.array(results[method]['R'])
fwer = np.mean(V > 0)
fdr_vals = np.where(R > 0, V / R, 0)
fdr = np.mean(fdr_vals)
power = np.mean(S) / m1 if m1 > 0 else 0
summary[method] = {'FWER': fwer, 'FDR': fdr, 'Power': power}
return summary
# シミュレーション実行
m = 20 # 検定の総数
m1 = 5 # 対立仮説が真の数
effect_size = 0.5
n_obs = 30 # 各検定の標本サイズ
n_sim = 5000
summary = simulate_multiple_testing(m, m1, effect_size, n_obs, n_sim)
# 結果の表示
print(f"シミュレーション設定: m = {m}, m₁ = {m1}, d = {effect_size}, n = {n_obs}")
print(f"シミュレーション回数: {n_sim}")
print("=" * 60)
print(f"{'Method':<15} {'FWER':>8} {'FDR':>8} {'Power':>8}")
print("-" * 45)
method_names = {'none': 'No correction', 'bonferroni': 'Bonferroni',
'holm': 'Holm', 'bh': 'BH (FDR)'}
for method in ['none', 'bonferroni', 'holm', 'bh']:
s = summary[method]
print(f"{method_names[method]:<15} {s['FWER']:>8.3f} {s['FDR']:>8.3f} {s['Power']:>8.3f}")
# 可視化
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
methods = ['none', 'bonferroni', 'holm', 'bh']
labels = ['No correction', 'Bonferroni', 'Holm', 'BH (FDR)']
colors = ['#e74c3c', '#3498db', '#2ecc71', '#9b59b6']
# FWER
fwer_vals = [summary[m]['FWER'] for m in methods]
axes[0].bar(labels, fwer_vals, color=colors, alpha=0.8)
axes[0].axhline(y=0.05, color='black', linestyle='--', linewidth=1.5, label='$\\alpha$ = 0.05')
axes[0].set_ylabel('FWER', fontsize=12)
axes[0].set_title('Family-Wise Error Rate', fontsize=13)
axes[0].legend()
axes[0].set_ylim(0, max(fwer_vals) * 1.2)
axes[0].tick_params(axis='x', rotation=30)
# FDR
fdr_vals = [summary[m]['FDR'] for m in methods]
axes[1].bar(labels, fdr_vals, color=colors, alpha=0.8)
axes[1].axhline(y=0.05, color='black', linestyle='--', linewidth=1.5, label='$q$ = 0.05')
axes[1].set_ylabel('FDR', fontsize=12)
axes[1].set_title('False Discovery Rate', fontsize=13)
axes[1].legend()
axes[1].set_ylim(0, max(fdr_vals) * 1.2)
axes[1].tick_params(axis='x', rotation=30)
# Power
power_vals = [summary[m]['Power'] for m in methods]
axes[2].bar(labels, power_vals, color=colors, alpha=0.8)
axes[2].set_ylabel('Power', fontsize=12)
axes[2].set_title('Average Power', fontsize=13)
axes[2].set_ylim(0, 1.05)
axes[2].tick_params(axis='x', rotation=30)
plt.tight_layout()
plt.show()
p値分布の可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
m = 1000 # 検定の総数
m1 = 100 # 対立仮説が真の数
m0 = m - m1 # 帰無仮説が真の数
effect_size = 0.5
n_obs = 30
# p値を生成
p_values = np.zeros(m)
# 帰無仮説が真: p値は一様分布
for i in range(m0):
sample = np.random.normal(0, 1, n_obs)
_, p_values[i] = stats.ttest_1samp(sample, 0)
# 対立仮説が真: p値は小さい値に偏る
for i in range(m0, m):
sample = np.random.normal(effect_size, 1, n_obs)
_, p_values[i] = stats.ttest_1samp(sample, 0)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# (1) すべてのp値のヒストグラム
axes[0].hist(p_values, bins=50, density=True, color='steelblue', alpha=0.7, edgecolor='white')
axes[0].axhline(y=1.0, color='red', linestyle='--', linewidth=1.5, label='Uniform(0,1)')
axes[0].set_xlabel('p-value', fontsize=12)
axes[0].set_ylabel('Density', fontsize=12)
axes[0].set_title(f'All p-values (m = {m})', fontsize=13)
axes[0].legend(fontsize=10)
# (2) 帰無仮説が真のp値
axes[1].hist(p_values[:m0], bins=50, density=True, color='green', alpha=0.7, edgecolor='white')
axes[1].axhline(y=1.0, color='red', linestyle='--', linewidth=1.5, label='Uniform(0,1)')
axes[1].set_xlabel('p-value', fontsize=12)
axes[1].set_ylabel('Density', fontsize=12)
axes[1].set_title(f'$H_0$ True (m₀ = {m0})', fontsize=13)
axes[1].legend(fontsize=10)
# (3) 対立仮説が真のp値
axes[2].hist(p_values[m0:], bins=50, density=True, color='orange', alpha=0.7, edgecolor='white')
axes[2].set_xlabel('p-value', fontsize=12)
axes[2].set_ylabel('Density', fontsize=12)
axes[2].set_title(f'$H_1$ True (m₁ = {m1})', fontsize=13)
plt.tight_layout()
plt.show()
# BH法の閾値と棄却の可視化
sorted_idx = np.argsort(p_values)
sorted_p = p_values[sorted_idx]
bh_threshold = np.arange(1, m+1) / m * 0.05
fig, ax = plt.subplots(figsize=(10, 6))
# 帰無仮説が真のp値と対立仮説が真のp値を色分け
is_null = sorted_idx < m0 # 元のインデックスがm0未満なら帰無仮説が真
ax.scatter(np.arange(1, m+1)[is_null], sorted_p[is_null],
s=10, color='green', alpha=0.5, label='$H_0$ True')
ax.scatter(np.arange(1, m+1)[~is_null], sorted_p[~is_null],
s=10, color='orange', alpha=0.5, label='$H_1$ True')
# BH法の閾値線
ax.plot(np.arange(1, m+1), bh_threshold, 'r-', linewidth=2, label='BH threshold ($q$ = 0.05)')
# 棄却される点
adj_bh, rej_bh = benjamini_hochberg(p_values, q=0.05)
n_rejected = np.sum(rej_bh)
ax.axvline(x=n_rejected, color='blue', linestyle='--', alpha=0.7,
label=f'Rejected: {n_rejected}')
ax.set_xlabel('Rank $k$', fontsize=13)
ax.set_ylabel('Sorted p-value $p_{(k)}$', fontsize=13)
ax.set_title('Benjamini-Hochberg Procedure', fontsize=14)
ax.legend(fontsize=10)
ax.set_xlim(0, m * 0.3)
ax.set_ylim(0, 0.1)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 結果の集計
V = np.sum(rej_bh[:m0]) # 偽陽性
S = np.sum(rej_bh[m0:]) # 真陽性
R = np.sum(rej_bh)
print(f"BH法の結果（q = 0.05）:")
print(f" 棄却数: {R}")
print(f" 真陽性: {S}")
print(f" 偽陽性: {V}")
print(f" FDR (実測): {V/R:.3f}" if R > 0 else " FDR: N/A")
print(f" 検出力: {S/m1:.3f}")
statsmodelsとの比較
import numpy as np
from statsmodels.stats.multitest import multipletests
# 先ほどの5つのp値で各補正法を適用
p_values = [0.005, 0.013, 0.031, 0.045, 0.250]
methods_sm = [
('bonferroni', 'Bonferroni'),
('holm', 'Holm'),
('fdr_bh', 'Benjamini-Hochberg'),
]
print("statsmodelsを使った多重比較補正")
print("=" * 60)
print(f"元のp値: {p_values}")
print()
for method_code, method_name in methods_sm:
rejected, adj_p, _, _ = multipletests(p_values, alpha=0.05, method=method_code)
print(f"【{method_name}】")
for i in range(len(p_values)):
print(f" p{i+1} = {p_values[i]:.3f} → adj_p = {adj_p[i]:.4f} → {'棄却' if rejected[i] else '棄却しない'}")
print()
まとめ
本記事では、多重比較問題とその補正法について解説しました。
- 多重比較問題: $m$ 個の検定を同時に行うとFWERが $1-(1-\alpha)^m$ に増大し、偽陽性が大量に発生する
- ボンフェローニ補正: 各検定の有意水準を $\alpha/m$ に下げる。簡単だが保守的
- Holm法: ステップダウン方式でFWERを制御。ボンフェローニより検出力が高い
- FDR: 棄却した中の偽陽性割合の期待値。FWERより緩やかな基準
- BH法: FDRを $q$ 以下に制御する手続き。大規模検定で有効
次のステップとして、以下の記事も参考にしてください。