偽発見率(FDR)の制御 — BH法の理論的基礎と実装

ゲノムワイド関連研究(GWAS)では、数百万の遺伝子座について「この遺伝子座は疾患と関連しているか」を同時に検定します。有意水準 $\alpha = 0.05$ で各検定を行うと、帰無仮説が正しい100万個の遺伝子座のうち約5万個が「偽陽性」として検出されてしまいます。真のシグナルが数十個しかない場合、検出結果のほぼすべてが偽陽性になります。

この多重検定問題(multiple testing problem)は、現代の大規模データ分析で最も重要な統計的課題の一つです。多重検定問題に対処する古典的な方法はボンフェローニ補正です。しかし、ボンフェローニ補正は非常に保守的で、真のシグナルまで見逃してしまうことが多くなります。検定数が百万を超えるGWASでは、ボンフェローニの閾値は $5 \times 10^{-8}$ 程度になり、多くの真の関連を見逃します。

この問題に対する革新的な解決策が、ベンジャミニとホッホベルグ(Benjamini and Hochberg, 1995)が提案した偽発見率(False Discovery Rate, FDR)の制御です。FDRは「発見(棄却)した中で偽陽性の割合」を制御するアプローチであり、ボンフェローニ補正よりもはるかに高い検出力を持ちます。

FDR制御を理解すると、以下のような応用が可能です。

  • ゲノミクス: 数千の遺伝子発現量の差異検定で、偽陽性を制御しつつ真のシグナルを検出します
  • 脳画像解析(fMRI): 数万のボクセル(立体画素)について活性化を検定する際に使われます
  • プロテオミクス: タンパク質の質量分析結果の統計的評価に標準的に利用されます
  • A/Bテスト: 複数の指標を同時に評価する場面で偽陽性を抑えつつシグナルを検出します

本記事の内容

  • 多重検定問題の本質と FWER vs FDR の違い
  • BH法(Benjamini-Hochberg法)の手順と理論的証明
  • q値(Storey の方法)の概念と計算
  • Pythonによる実装とシミュレーション

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

多重検定問題の本質

偽陽性の蓄積

$m$ 個の独立な検定を有意水準 $\alpha$ で行うとき、少なくとも1つの偽陽性が生じる確率は次のようになります。

$$ P(\text{少なくとも1つの偽陽性}) = 1 – (1 – \alpha)^{m_0} $$

ここで $m_0$ は帰無仮説が正しい検定の数です。$m_0 = 100$, $\alpha = 0.05$ のとき、この確率は $1 – 0.95^{100} \approx 0.994$ です。つまり、ほぼ確実に偽陽性が生じます。$m_0 = 1000$ では $1 – 0.95^{1000} \approx 1.0$ であり、偽陽性が生じないことはまず考えられません。

FWER(家族内過誤率)

多重検定への古典的なアプローチは、家族内過誤率(Family-Wise Error Rate, FWER)を制御することです。

$$ \text{FWER} = P(\text{少なくとも1つの偽陽性}) $$

ボンフェローニ補正は $\text{FWER} \leq \alpha$ を保証する最も単純な方法で、各検定の有意水準を $\alpha / m$ に設定します。しかし、$m$ が大きいと有意水準が極端に小さくなり、真のシグナルも検出できなくなります。

FDR(偽発見率)

FDRは、FWERよりも緩やかな基準です。

$$ \begin{equation} \text{FDR} = E\left[\frac{V}{R \vee 1}\right] = E\left[\frac{\text{偽陽性の数}}{\text{棄却した数}}\right] \end{equation} $$

ここで $V$ は偽陽性の数(帰無仮説が正しいのに棄却した数)、$R$ は棄却の総数、$R \vee 1 = \max(R, 1)$ です。

FDRは「発見(棄却されたもの)のうちどれだけの割合が偽物か」を制御します。$\text{FDR} = 0.05$ は「棄却した検定のうち、平均して $5\%$ が偽陽性である」ことを意味します。100個棄却したなら、そのうち約5個が偽陽性です。

FWERは「1つでも偽陽性が出る確率」を制御するのに対し、FDRは「偽陽性の割合」を制御します。多くの検定を同時に行い、一定割合の偽陽性を許容できる場面では、FDR制御の方がはるかに高い検出力を実現します。

FDRとFWERの関係

FDRとFWERの間には以下の関係が成り立ちます。

  1. FDR $\leq$ FWER: FDRはFWER以下です。FWERを制御する手法は自動的にFDRも制御します。したがってボンフェローニ補正はFDRも制御しますが、過度に保守的です
  2. すべての帰無仮説が正しい場合($m_0 = m$): FDR = FWER です。棄却が1つでもあればそれは必ず偽陽性なので、偽陽性の割合は0か1になり、その期待値がFWERに一致します
  3. FDRは弱い制御: FDRは「偽陽性の割合の期待値」であり、個々の実験での偽陽性の割合は確率的に変動します。「FDR $= 0.05$」は「平均して5%が偽陽性」であり、ある実験では10%、別の実験では0%かもしれません

具体例で理解するFDR

1000個の遺伝子について発現量の差を検定し、50個が真にdifferentially expressed(差がある)だとします。

ボンフェローニ補正($\alpha = 0.05$)の場合、各検定の閾値は $0.05/1000 = 5 \times 10^{-5}$ です。これは非常に厳しく、効果量が大きいシグナルしか検出できません。仮に20個のシグナルを検出できたとします。FWERは0.05以下が保証され、偽陽性はほぼゼロです。

BH法($\alpha = 0.05$)の場合、閾値は適応的に決まり、より多くのシグナルを検出できます。仮に40個のシグナルを検出したとすると、そのうち約2個($40 \times 0.05 = 2$)が偽陽性です。検出力は80%($40/50$)に達し、ボンフェローニの40%($20/50$)を大きく上回ります。

多くのゲノミクスや神経科学の研究では、偽陽性が2個混じっていても40個の真のシグナルを発見できる方が、確実に偽陽性ゼロだが20個しか発見できないより有益です。FDR制御が広く採用されている理由はここにあります。

FDRの概念を理解したところで、これを制御するBH法の手順を見ましょう。

BH法(Benjamini-Hochberg法)

手順

BH法は驚くほどシンプルな手順です。

ステップ1: $m$ 個の検定のp値を昇順にソートする: $p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}$

ステップ2: 次の条件を満たす最大の $k$ を見つける:

$$ \begin{equation} p_{(k)} \leq \frac{k}{m} \alpha \end{equation} $$

ステップ3: $p_{(1)}, p_{(2)}, \dots, p_{(k)}$ に対応する帰無仮説を棄却する。条件を満たす $k$ が存在しなければ、何も棄却しない。

BH法の幾何学的理解

BH法の手順は幾何学的に理解すると直感的です。横軸にランク $k$、縦軸にソートされたp値 $p_{(k)}$ をプロットします。BH法の閾値 $\frac{k}{m}\alpha$ は原点を通る直線(傾き $\alpha/m$)です。

このプロットでは、p値の点が直線の下にある限り棄却し、直線を超えた時点で棄却を打ち切ります。ボンフェローニ補正は水平線 $\alpha/m$($k$ に依存しない)であるのに対し、BH法の閾値は $k$ とともに増加するため、より多くの棄却を許容します。

BH法の適応的な閾値の直感は以下の通りです。もし多くの真のシグナルがあれば、小さいp値が多数存在し、ランク $k$ が大きくても直線の下に留まります。すなわち、BH法は「シグナルが豊富な場合は閾値を緩く、シグナルが少ない場合は閾値を厳しく」自動調整するのです。

BH法がFDRを制御する証明のスケッチ

ベンジャミニとホッホベルグの定理: p値が独立であるとき、BH法は $\text{FDR} \leq \frac{m_0}{m}\alpha \leq \alpha$ を保証します。

証明の核心的なアイデアは次の通りです。$m_0$ 個の真の帰無仮説のp値は $[0, 1]$ 上の一様分布に従います。BH法の閾値 $\frac{k}{m}\alpha$ は $k$(棄却数)が大きいほど緩くなります。真の帰無仮説のp値が閾値を超える確率は、一様分布の性質から $\frac{k}{m}\alpha$ 以下です。全体の帰無仮説の数 $m_0$ を掛けると、偽陽性の期待数は $\frac{m_0 k}{m}\alpha$ 以下になり、棄却数 $k$ で割るとFDRは $\frac{m_0}{m}\alpha$ 以下になります。

補正p値(adjusted p-value)

BH法の結果は補正p値(adjusted p-value, BH q-value)として表現できます。

$$ \begin{equation} \tilde{p}_{(k)} = \min_{j \geq k} \left\{\min\left(\frac{m}{j} p_{(j)},\, 1\right)\right\} \end{equation} $$

補正p値が $\alpha$ 以下の検定を棄却すれば、BH法と同じ結果が得られます。

BH法のFDR制御の証明

BH法がFDRを制御することの証明は、確率論の美しい結果です。ここでは、p値が独立な場合の証明を詳しく行います。

問題の設定

$m$ 個の検定のうち、$m_0$ 個は帰無仮説が真($H_{0i}$ は真)、$m_1 = m – m_0$ 個は対立仮説が真とします。帰無仮説が真の検定のp値は $U(0, 1)$ に従い、すべてのp値は独立です。

BH法の閾値は $\frac{k}{m}\alpha$ であり、棄却数は $R$、偽陽性の数は $V$ です。

証明の核心

帰無仮説が真の検定 $i$($i \in \mathcal{H}_0$)について、BH法でこの検定が棄却される確率を考えます。

検定 $i$ のp値を $p_i$ とし、他のすべての検定の結果を条件づけます。他の検定の結果が与えられたとき、BH法の閾値は「$p_i$ を除いた残りのp値」によって(ほぼ)決まります。

鍵となる不等式は以下です。帰無仮説が真の検定 $i$ が棄却される確率を、棄却数 $R$ で条件づけて評価します。

$$ E\left[\frac{V}{R \vee 1}\right] = E\left[\sum_{i \in \mathcal{H}_0} \frac{\mathbf{1}(i \text{が棄却})}{R \vee 1}\right] = \sum_{i \in \mathcal{H}_0} E\left[\frac{\mathbf{1}(i \text{が棄却})}{R \vee 1}\right] $$

各項について、$p_i \sim U(0, 1)$ であることと、棄却条件 $p_i \leq \frac{R}{m}\alpha$ を使うと、

$$ E\left[\frac{\mathbf{1}(p_i \leq R\alpha/m)}{R}\right] \leq E\left[\frac{1}{R} \cdot \frac{R\alpha}{m}\right] = \frac{\alpha}{m} $$

最後の等式は、$p_i$ が一様分布であるため $P(p_i \leq c) = c$ が成り立つことを使っています(ここでは証明の核心的なアイデアを示しており、条件付き期待値の厳密な扱いは省略しています)。

$m_0$ 個の帰無仮説について合計すると、

$$ \text{FDR} = E\left[\frac{V}{R \vee 1}\right] \leq m_0 \cdot \frac{\alpha}{m} = \frac{m_0}{m}\alpha \leq \alpha $$

$m_0 \leq m$ なので $\text{FDR} \leq \alpha$ が保証されます。

正の依存性の場合(BY法)

ベンジャミニとイェクティエリ(Benjamini and Yekutieli, 2001)は、p値が正の回帰依存性(PRDS: positive regression dependency on each subset)を持つ場合にもBH法のFDR制御が成り立つことを示しました。

正規分布に基づく片側検定や、多変量正規から生成されるp値(正の相関の場合)はPRDSを満たすため、BH法が適用できます。

一般の依存構造の場合には、BY法(Benjamini-Yekutieli法)が使えます。BY法の閾値は $\frac{k}{m \cdot c(m)}\alpha$ であり、$c(m) = \sum_{k=1}^m 1/k \approx \ln m + 0.577$ は調和級数です。BH法より保守的ですが、任意の依存構造でFDR制御を保証します。

q値 — StoreyのFDR推定

q値の概念

BH法の補正p値は「その検定を棄却する際に保証されるFDRの最小レベル」を表しますが、帰無仮説が真の割合 $\pi_0 = m_0/m$ を考慮していません。ストーリー(Storey, 2002, 2003)は、$\pi_0$ を推定してFDRの精度を向上させる方法を提案しました。

$\pi_0$ の推定のアイデアは次の通りです。p値のヒストグラムを描くと、帰無仮説が真のp値は $[0, 1]$ 上で一様に分布し、対立仮説のp値は0付近に集中します。したがって、p値が大きい領域(例えば $p > \lambda$ の領域)は帰無仮説のp値がほとんどを占めるため、

$$ \hat{\pi}_0(\lambda) = \frac{\#\{p_i > \lambda\}}{m(1 – \lambda)} $$

で $\pi_0$ を推定できます。$\lambda$ の選択にはブートストラップ法が使われます。

q値の定義

各検定のq値は、その検定を棄却に含めたときの推定FDRの最小値として定義されます。

$$ q(p_i) = \min_{t \geq p_i} \hat{\text{FDR}}(t) = \min_{t \geq p_i} \frac{\hat{\pi}_0 \cdot m \cdot t}{\#\{p_j \leq t\} \vee 1} $$

q値はBH法の補正p値に $\hat{\pi}_0$ の補正を加えたものと考えることができ、$\hat{\pi}_0 < 1$ のとき(つまり対立仮説が多いとき)、BH法よりも緩い閾値で検出力を上げられます。

ボンフェローニ補正との比較

理論的な比較

特性 ボンフェローニ BH法 BY法
制御する量 FWER FDR FDR
閾値 $\alpha / m$ $\frac{k}{m}\alpha$(適応的) $\frac{k}{mc(m)}\alpha$
仮定 なし 独立 or PRDS 任意の依存構造
検出力 低い($m$ が大きいと非常に保守的) 高い 中間的
偽陽性の扱い 1つも許さない 割合として許容 割合として許容

ホルム法とBH法の比較

FWERを制御する方法としてはボンフェローニよりもホルム法(Holm, 1979)の方が検出力が高く、実用上はボンフェローニの代わりにホルム法を使うべきです。ホルム法はp値をソートした後、ステップダウンで $p_{(k)} \leq \frac{\alpha}{m – k + 1}$ を判定します。ただし、それでもBH法の検出力には及びません。

検定数が異なる場面での比較

多重検定補正法の選択は、検定の数 $m$ に大きく依存します。

$m$ が小さい場合($m \leq 20$): ボンフェローニ補正やホルム法でもそれほど保守的にならず、FWERとFDRの差は小さいです。臨床試験の副次評価項目の比較など。

$m$ が中程度の場合($20 < m \leq 1000$): FDR制御の優位性が現れ始めます。遺伝子パネル検査、プロテオミクスのターゲット解析、多変量A/Bテストなど。

$m$ が大きい場合($m > 1000$): FDR制御がほぼ必須です。ゲノムワイド関連研究(GWAS)、RNA-seq、全脳fMRI解析など。ボンフェローニ補正では閾値が $\alpha/m < 10^{-7}$ のオーダーになり、効果量が非常に大きいシグナル以外は検出できません。

FDRの限界

FDR制御にもいくつかの限界があります。

  1. 偽発見率は期待値: $\text{FDR} = 0.05$ は平均であり、個々の実験では偽陽性の割合が大きく変動しうる
  2. 帰無仮説がすべて真の場合: 棄却が少数の場合にFDRの推定が不安定になる
  3. p値の品質に依存: p値が帰無仮説のもとで一様分布に従わない場合(例: 検定の前提条件の違反)、FDR制御は保証されない

Pythonによる実装

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

def bh_correction(p_values, alpha=0.05):
    """Benjamini-Hochberg法"""
    m = len(p_values)
    sorted_idx = np.argsort(p_values)
    sorted_p = p_values[sorted_idx]

    # BH閾値
    thresholds = np.arange(1, m + 1) / m * alpha

    # 棄却する最大のk
    rejected = sorted_p <= thresholds
    if np.any(rejected):
        k_max = np.max(np.where(rejected)[0]) + 1
        reject_mask = np.zeros(m, dtype=bool)
        reject_mask[sorted_idx[:k_max]] = True
    else:
        reject_mask = np.zeros(m, dtype=bool)

    # 補正p値
    adjusted_p = np.zeros(m)
    adjusted_p[sorted_idx[-1]] = min(sorted_p[-1] * m / m, 1.0)
    for i in range(m - 2, -1, -1):
        adjusted_p[sorted_idx[i]] = min(
            sorted_p[i] * m / (i + 1),
            adjusted_p[sorted_idx[i + 1]]
        )
    adjusted_p = np.minimum(adjusted_p, 1.0)

    return reject_mask, adjusted_p

def bonferroni_correction(p_values, alpha=0.05):
    """ボンフェローニ補正"""
    m = len(p_values)
    adjusted_p = np.minimum(p_values * m, 1.0)
    reject_mask = adjusted_p <= alpha
    return reject_mask, adjusted_p

# --- シミュレーション ---
np.random.seed(42)
m = 1000       # 検定の総数
m1 = 50        # 対立仮説が正しい数
m0 = m - m1    # 帰無仮説が正しい数
alpha = 0.05

# p値の生成
p_null = np.random.uniform(0, 1, m0)  # 帰無仮説: U(0,1)
effect_sizes = np.random.uniform(2, 4, m1)
p_alt = 2 * (1 - stats.norm.cdf(effect_sizes))  # 対立仮説: 小さいp値
p_values = np.concatenate([p_null, p_alt])
true_null = np.concatenate([np.ones(m0), np.zeros(m1)]).astype(bool)

# BH法
bh_reject, bh_adj = bh_correction(p_values, alpha)
# ボンフェローニ
bonf_reject, bonf_adj = bonferroni_correction(p_values, alpha)

# 結果
for name, reject in [("BH法", bh_reject), ("ボンフェローニ", bonf_reject)]:
    tp = np.sum(reject & ~true_null)
    fp = np.sum(reject & true_null)
    fn = np.sum(~reject & ~true_null)
    total_reject = np.sum(reject)
    fdr = fp / total_reject if total_reject > 0 else 0
    power = tp / m1
    print(f"{name:12}: 棄却={total_reject}, TP={tp}, FP={fp}, "
          f"FDR={fdr:.3f}, 検出力={power:.3f}")

# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

# p値のソートプロット + BH閾値
ax = axes[0]
sorted_p = np.sort(p_values)
k = np.arange(1, m + 1)
ax.scatter(k, sorted_p, s=3, alpha=0.5, color="steelblue", label="Sorted p-values")
ax.plot(k, k / m * alpha, "r-", linewidth=2, label=f"BH threshold (α={alpha})")
ax.axhline(alpha / m, color="green", linestyle="--", linewidth=1.5,
           label=f"Bonferroni (α/m={alpha/m:.1e})")
ax.set_xlabel("Rank (k)", fontsize=12)
ax.set_ylabel("p-value", fontsize=12)
ax.set_title("BH procedure", fontsize=13)
ax.legend(fontsize=9)
ax.set_xlim(0, 150)
ax.set_ylim(0, 0.15)
ax.grid(True, alpha=0.3)

# FDRと検出力のトレードオフ
ax = axes[1]
alpha_range = np.linspace(0.001, 0.3, 50)
fdr_bh, power_bh = [], []
fdr_bonf, power_bonf = [], []

for a in alpha_range:
    bh_r, _ = bh_correction(p_values, a)
    tp = np.sum(bh_r & ~true_null)
    fp = np.sum(bh_r & true_null)
    total_r = np.sum(bh_r)
    fdr_bh.append(fp / total_r if total_r > 0 else 0)
    power_bh.append(tp / m1)

    bonf_r, _ = bonferroni_correction(p_values, a)
    tp = np.sum(bonf_r & ~true_null)
    fp = np.sum(bonf_r & true_null)
    total_r = np.sum(bonf_r)
    fdr_bonf.append(fp / total_r if total_r > 0 else 0)
    power_bonf.append(tp / m1)

ax.plot(fdr_bh, power_bh, "b-o", markersize=3, linewidth=1.5, label="BH (FDR)")
ax.plot(fdr_bonf, power_bonf, "r-s", markersize=3, linewidth=1.5, label="Bonferroni (FWER)")
ax.set_xlabel("False Discovery Rate", fontsize=12)
ax.set_ylabel("Power (True Positive Rate)", fontsize=12)
ax.set_title("FDR vs Power tradeoff", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("fdr_bh_method.png", dpi=150, bbox_inches="tight")
plt.show()

このシミュレーション結果から、以下の重要な知見が得られます。

  1. BH法はボンフェローニ補正よりもはるかに多くの真のシグナルを検出する。FDR制御は偽陽性の「割合」を制御するため、多くの発見を許容しながらもその品質(偽陽性の割合)を保証します

  2. 左図のp値プロットはBH法の幾何学的な解釈を示している。ソートされたp値(青い点)がBH閾値の直線(赤線)を下回る範囲が棄却域です。ボンフェローニの閾値(緑の破線)は非常に低い位置にあり、少数の検定しか棄却しません

  3. 右図のトレードオフ曲線はBH法の優位性を示している。同じFDRに対してBH法の方が高い検出力を達成しています。特に、FDR $= 0.05$ の周辺でBH法の優位性が顕著です

FDR制御の実践的なガイドライン

いつFDR制御を使うべきか

FDR制御は以下の条件が揃ったときに特に有効です。

  1. 大量の仮説を同時に検定する場合: 数百〜数百万の検定。遺伝子発現、fMRIボクセル、A/Bテストの複数指標など
  2. 一定割合の偽陽性を許容できる場合: 検出された候補を次のステップ(実験的検証、追試験など)でフィルタリングできる場合
  3. FWERが過度に保守的な場合: 検定数が多いとボンフェローニ補正では何も検出できなくなる

逆に、以下の場合はFWER制御(ボンフェローニやホルム法)の方が適切です。

  • 少数の事前に定義された比較(例: 臨床試験の主要エンドポイント)
  • 1つの偽陽性も許容できない場合(安全性に関わる検定など)

p値ヒストグラムによる診断

FDR制御を適用する前に、p値のヒストグラムを確認することが重要な診断ステップです。

理想的なパターン: $[0, 1]$ 上の一様な基盤(帰無仮説のp値)の上に、0付近にスパイクがある(対立仮説のp値)。このパターンではBH法が有効に機能します。

問題のあるパターン1(右偏り): p値が1に近い方に偏っている場合、検定統計量の帰無分布が不正確である可能性があります。このパターンではBH法を適用する前に検定手法自体を見直すべきです。

問題のあるパターン2(U字型): p値が0と1の両方にピークを持つ場合、離散的なp値や、検定の前提条件の違反が疑われます。

局所的FDR

BH法のFDRは「棄却された検定全体」での偽陽性の割合を制御しますが、個々の検定について「この特定の検定が偽陽性である確率」を知りたいことがあります。これを局所的FDR(local FDR, fdr)と呼びます。

$$ \text{fdr}(p) = \frac{\pi_0 f_0(p)}{f(p)} $$

ここで $f_0(p)$ は帰無仮説のp値の密度(一様分布なので1)、$f(p)$ はp値の混合密度です。Efron(2004, 2008)による経験ベイズ的アプローチが広く用いられています。

まとめ

本記事では、偽発見率(FDR)の概念とBH法の理論を解説しました。

  • FDRは「発見(棄却)の中で偽陽性の割合の期待値」であり、FWERよりも緩やかな基準です
  • BH法はp値を昇順にソートして直線 $\frac{k}{m}\alpha$ と比較するシンプルな手順で、$\text{FDR} \leq \alpha$ を保証します
  • ボンフェローニ補正と比べてBH法は大幅に高い検出力を持つ。特に検定数が多く、真のシグナルが少数存在する場面で差が顕著です
  • q値は各検定に対するFDR制御のレベルを表し、帰無仮説の割合 $\pi_0$ を推定することでBH法よりも高い検出力を実現します
  • 局所的FDRは個々の検定が偽陽性である確率を推定する経験ベイズ的手法であり、各検定の信頼性を個別に評価できます
  • p値ヒストグラムの診断はFDR制御を適用する前に実施すべき重要なステップです。帰無分布の一様性が破れている場合、BH法の保証が成り立たない可能性があります

FDR制御の考え方は、多重検定問題を扱うあらゆる場面——ゲノミクス、脳画像、計量経済学、オンライン実験など——で現代統計学の基本ツールとなっています。

次のステップとして、以下の記事も参考にしてください。