FDR (False Discovery Rate) 制御は多重検定補正の一つです。FDR制御を行うための代表的なアルゴリズムがBenjamini-Hochberg法です。以下ではBH法と略します。

FDRは以下の式で表されます。TPはTrue Positive, FPはFalse Positiveです。FDR制御ではFDRの期待値を指定した値以下に抑えます。

$$FDR=\frac{FP}{TP+FP}$$

FDR制御とBH法は沢山のウェブサイトで紹介されていますので、詳細は割愛します。

このページでは、何故BH法でFDR制御できるのかを図で説明します。私の担当講義で話していますが、せっかくなのでQiitaにも載せることにしました。

FDR制御の説明

以下では、同じタイプの仮説検定を $N$ 回行うとします。例えば、$t$ 検定を$10^4$ 回行うとします。各仮説検定は互いに独立とします。簡単のため、仮説検定は片側検定とします。

図の上側の分布は、帰無仮説が正しい場合に検定統計量が従う分布です。例えば、$t$ 検定の場合、検定統計量は $t$ 値、その従う分布は $t$ 分布です。

図の下側の分布は、$N$ 回の検定によって得られた $N$ 個の検定統計量の実際の分布です。

閾値を図に示す値にした場合を考えます。実際に陽性と判定されるのは、下図のBの部分です。例えば、ここに $500$ 個いるとします。この中にはTPとFPが混ざっています。

このときのFPの個数の上限を上の図から見積もります。もし仮に、$N$ 個全ての仮説検定で帰無仮説が正しかった場合、$A$ の面積を $p$ として、FPの個数の期待値は $Np$ と表されます。例えば、$N=10^4$, $p=0.01$ の場合、$100$ 個になります。また、$N$ 個の中に対立仮説が正しい場合が含まれているときも、FPの個数の期待値は $Np$ より大きくなることはありません。

例えば、$500$ 個が出力されて、そのうち高々 $100$ 個が偽物だと予想されるなら、FDRの期待値の上限は $100/500=0.2$ になります。

あとは、指定したFDR (例えば $0.05$) を超えない範囲で、閾値をなるべく沢山の出力が得られるような値 (この図でいえば、なるべく左側の値) にすればおしまいです。

BH法との対応関係

ここまでは検定統計量の分布を使って説明をしてきました。一方、BH法は $p$ 値ベースのアルゴリズムです。対応関係を整理します。

$N$ 個の仮説検定で得られた $p$ 値を小さい順に並べます。簡単のため等号は含まないとします。

$$p_1<p_2 < \cdots < p_N$$

先ほどの図でいうと、下側の分布で右端の方から順に $p_1, p_2,\ldots$ に対応する検定統計量が並んでいます。

$1$ 番目から $i$ 番目までを陽性として出力する場合を考えます。出力される個数は $i$ です。これは先ほどの図のBの部分です。また、先ほどの図のAの面積は $p_i$ に相当します。従って、このときのFDRの期待値の上限 ($q$ 値) は以下のように表されます。

$$q_i=\frac{Np_i}{i}$$

あとは、この値が指定値を超えない範囲で最大の $i$ を見つければ良いので、BH法のアルゴリズムが出てきます。

Pythonのコード

以下に $q$ 値を計算するPythonのコードを載せます。単調増加となるように np.minimum.accumulate をかけています。

def calculate_q(p_seq):
  p_arr   = np.asarray(p_seq)
  N       = len(p_arr)
  idx_arr = np.argsort(p_arr)
  q_arr   = p_arr[idx_arr] * N / (np.arange(N) + 1)
  q_arr   = np.minimum.accumulate(q_arr[::-1])[::-1]
  q_arr[idx_arr] = q_arr.copy()
  return q_arr

数値実験

以下に数値実験結果を載せます。実験条件は以下のとおりです。

名前	設定値
検定の種類	ウェルチの $t$ 検定
各標本のサイズ $n$	$10$ から $100$ まで $10$ 刻み
仮説検定の数 $N$	$10^4$
真に陽性の数	$1000$
陽性の場合の効果量	$0.8$
FDRの閾値	$0.05$
試行回数	$10$

陽性として出力される数 (TP+FP) は、標本サイズに応じて増加した後に一定値に収束しました。収束値が真の陽性数 $1000$ を超えているのはFPを含むためと考えられます。

FDRは、どの標本サイズでも指定値 $0.05$ より平均的に低い値でした。