ほくそ笑む

R言語と統計解析について

Benjamini-Hochberg法の原理

はじめに

FDR(False Discovery Rate)は多重検定を行う際の有意水準指標の一つです。
FDR を制御する最も簡単な方法は BH法(Benjamini-Hochberg法) です。
詳しくは下記をご覧ください。
Bonferroni法、Holm法、False Discovery Rate | 大阪大学腎臓内科
さて、BH 法は FDR の有意水準を 0.05 のようにあらかじめ決めてから行うのですが、実用上これは不便です。
なので、q-value という p値のように扱うことのできる FDR の指標をそれぞれの帰無仮説に対して求めるのが普通です。
BH 法を用いた q-value の求め方は以下の通りです。

  1. N個の帰無仮説を、p値の小さい順に並べます( p_1 \leq p_2 \leq p_3 \leq \dots \leq p_N)。
  2.  q_m = p_m \times N / m m = 1, \dots , N に対して求めます。
  3. i,j = 1, \dots , N かつ  i < j に対して、 q_i > q_j ならば  q_i = q_j とします。

q-value は、 p_i に対応する帰無仮説 H_i とすると、 H_1, \dots, H_m までの全てを棄却した場合の FDR が  q_m になります。
3 は  q_1 \leq q_2 \leq \dots \leq q_N を成り立たせるための補正です。
重要なのは 2 の  q_m = p_m \times N / m です。
なぜこれで FDR を計算できるのでしょうか?
今日はこの原理について説明しようと思います。

BH法の原理

FDR は一言で言うと「棄却された帰無仮説のうち、本当は帰無仮説が正しいものの割合」です。
もし、 H_1, \dots, H_m を棄却した場合、棄却された帰無仮説の数は m ですから、上記の  q_m = p_m \times N / m の分母については説明がつきます。
では、分子の  p_m \times N についてはどうでしょうか?
これは「棄却されたのに本当は帰無仮説が正しいもの」の数になるはずです。
この「棄却されたのに本当は帰無仮説が正しいもの」の数は、実際にはわからないので推定する必要があります。
これを推定するのに、BH法は「帰無仮説が全て正しい場合、p値は一様分布する」という原理を使用しています。
この原理は、非常に強力なもので、いかなる検定の p値に対しても成り立ちます。
少しシミュレーションで確認してみましょう。統計言語 R を使います。

N <- 10000
p.values <- double(N)
for(i in 1:N) {
  x <- rnorm(10)
  y <- rnorm(10)
  p.values[i] <- t.test(x,y)$p.value
}
hist(p.values)


有意差の無いデータに対して t検定を10000回繰り返して p値の分布を出力してみました。
みごとに一様分布しています。
このように、本当は有意差が無いデータでも、p値が棄却域に達することが一定の割合で起こります。
その分布が一様分布することから、上記で言及した「棄却されたのに本当は帰無仮説が正しいもの」の数を推定することができます。
上記のヒストグラムで具体的に求めると、p <= 0.05 となる帰無仮説の数は、一番左の棒だけですので、およそ500、p <= 0.1 となると、その右の棒も合わさるのでおよそ1000です。
一般に、 p \leq \alpha となる帰無仮説の数の期待値は  \alpha \times N となります。
BH法では、「棄却されたのに本当は帰無仮説が正しいもの」の数の推定値として、この期待値を使用します。
今、 H_1, \dots, H_m を棄却したとき、「棄却されたのに本当は帰無仮説が正しいもの」の数の推定値、すなわち、一様分布で  p \leq p_m となる帰無仮説の数の期待値は、上記より  p_m \times N です。
これは上記BH法 2 の  q_m = p_m \times N / m の分子と一致しました。
結局この式は、「棄却された帰無仮説のうち、本当は帰無仮説が正しいものの割合」、すなわち FDR を推定していることがわかりました。

BH法の改良

以上でBH法の原理についてはわかりました。
しかし、疑問が残ります。
実際に多重検定する場合、真に有意差のあるデータを含んでいると考えられるため、p値は一様分布しません。
実際は p=0 側に偏った次のようなグラフになると考えられます。

この分布は、「真に有意差のあるデータの分布」と「有意差の無いデータの分布(一様分布)」の混合分布と考えられます。
この場合、一様分布をなす帰無仮説の数が減少します。
「真に有意差のあるデータ」の数を M とすると、一様分布をなす帰無仮説の数は N - M になり、「棄却されたのに本当は帰無仮説が正しいもの」の数の推定値は  p_m \times (N-M) とするべきです。
まあ実際は M を正確に知ることはできないので、どうにか推定する必要がありそうです。
BH法のここらへんの改良が Storey らによって行われています。
(PDF)A direct approach to false discovery rates
Storey : The positive false discovery rate: a Bayesian interpretation and the q-value
Statistical significance for genomewide studies | PNAS
このへんの論文もそのうち読んでみようと思います。

おわりに

p値が一様分布で無い場合*1、BH法で推定された「棄却されたのに本当は帰無仮説が正しいもの」の数は過剰に見積もりされていることになります。
すなわち、BH法で求められる FDR は実際よりも大きめに算出される傾向にあります。*2
FDR はそもそも多重検定の有意水準指標としては「ゆるすぎる」という批判もあるようなので*3、これぐらい保守的でもあまり気にすることはなさそうです。
プログラム書く身としては、BH法は実装も簡単だし、気軽に使えるという点では便利な手法です。
あまりまとまってませんが、今日はこれまで。

*1:ほとんどの場合がそうですが

*2:保守的といいます

*3:http://d.hatena.ne.jp/what_a_dude/20111226/p1