ほくそ笑む

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

もっと保守的なBH法思いついた

前回、FDR 算出方法の一つ、BH法(Benjamini-Hochberg法)の原理について書きました。
そこで、

一般に、 p \leq \alpha となる帰無仮説の数の期待値は  \alpha \times N となります。
BH法では、「棄却されたのに本当は帰無仮説が正しいもの」の数の推定値として、この期待値を使用します。

と書きました。
ところで、この  p \leq \alpha となる帰無仮説の数はどういう分布をしているのでしょうか?
統計言語 R でシミュレーションしてみましょう。

counts <- integer(1000)
for(i in 1:1000) {
  p.values <- double(100)
  for(j in 1:100) {
    x <- rnorm(10)
    y <- rnorm(10)
    p.values[j] = t.test(x,y)$p.value
  }
  counts[i] <- sum(p.values <= 0.05)
}
hist(counts)

「有意差の無いデータに対して、t検定を 100回行い、そのうち p <= 0.05 となる数をカウントする」という試行を 1000回やってみます。
分布は次のようになりました。

このグラフを見て気付いたんですが、帰無仮説が正しい場合、一回の検定というのは、成功確率 0.05 のベルヌーイ試行になります。
つまり、この分布は二項分布  Bi(100, 0.05) になるはずです。

もっと保守的なBH法

 p \leq \alpha となる帰無仮説の数が二項分布  Bi(N, \alpha) に従うことがわかりました。
二項分布  Bi(N, \alpha) の期待値は  \alpha \times N ですので、前回と整合性がとれます。
通常の BH法では、「棄却されたのに本当は帰無仮説が正しいものの数」の推定値として、この期待値  \alpha \times N を使います。
しかし、それでよいのでしょうか?
例えば上記のヒストグラム \alpha = 0.05, N = 100 のときの分布です。
このとき、 \alpha \times N = 5 であり、「棄却されたのに本当は帰無仮説が正しいものの数」の推定値として 5 が使われます。
しかし、ヒストグラムを見ると、5 より大きい値となる可能性が十分あるではありませんか。
これでは FDR を低く見積もってしまう可能性が見逃せません。
我々はもっと保守的になるべきではないでしょうか?
というわけで、二項分布が取り得る値を十分カバーするように、「棄却されたのに本当は帰無仮説が正しいものの数」の推定値として、二項分布の 99パーセンタイル点を使います。
統計言語 R には二項分布のパーセンタイル点を計算する関数があるので、それで計算してみましょう。

perc99 <- qbinom(p = 0.99, size = 100, prob = 0.05)
perc99

[1] 11

これにより、上記の場合、「棄却されたのに本当は帰無仮説が正しいものの数」の推定値として 11 を取れば、99% は低く見積もる可能性を排除できるということがわかりました。
非常に保守的なBH法の完成です。

Rで実装

早速 R で実装してみましょう。

p.adjust.conservativeBH <- function(p.values, perc = 0.99) {
  n <- length(p.values)
  if (n <= 1) 
    return(p.values)
  i <- n:1L
  o <- order(p.values, decreasing = TRUE)
  ro <- order(o)
  pmin(1, cummin(pmax(qbinom(p = perc, size = n, prob = p.values[o]), p.values[o] * n) / i))[ro]
}

二項分布の 99パーセンタイル点は 0 になる可能性もあるので、そのときは期待値を取るように pmax() 関数で大きい方を選択しています。
適当にでっちあげた p値に対して適用した結果が下記の図です(左が全体、右はpが小さい部分の拡大図)。

ちゃんと保守的に FDR を算出できているようです。
以上です。