Hatena::ブログ(Diary)

ほくそ笑む このページをアンテナに追加 RSSフィード

2012-04-13

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)

f:id:hoxo_m:20120413065004p:image

有意差の無いデータに対して 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 側に偏った次のようなグラフになると考えられます。

f:id:hoxo_m:20120413070123p:image

この分布は、「真に有意差のあるデータの分布」と「有意差の無いデータの分布(一様分布)」の混合分布と考えられます。

この場合、一様分布をなす帰無仮説の数が減少します。

「真に有意差のあるデータ」の数を 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

このへんの論文もそのうち読んでみようと思います。

おわりに

p値が一様分布で無い場合*1、BH法で推定された「棄却されたのに本当は帰無仮説が正しいもの」の数は過剰に見積もりされていることになります。

すなわち、BH法で求められる FDR は実際よりも大きめに算出される傾向にあります。*2

FDR はそもそも多重検定の有意水準指標としては「ゆるすぎる」という批判もあるようなので*3、これぐらい保守的でもあまり気にすることはなさそうです。

プログラム書く身としては、BH法は実装も簡単だし、気軽に使えるという点では便利な手法です。

あまりまとまってませんが、今日はこれまで。

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

*2:保守的といいます

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

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証