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$ より平均的に低い値でした。