FDR制御法(ST)

Benjamini らが定義した BH 法 では全仮説が帰無仮説と仮定している。一方、Storey and Tibshirani が定義した FDR 制御法では、すべての仮説を帰無仮説と仮定せずに、帰無仮説が含まれる割合を線形予測してから q-value を計算する。

およそ次のような手順で q-value を計算する。

  1. p-value を昇順に並べる。
  2. λ=0, 0.01, 0.02, ..., 0.95 として、次の式に基づいて π0(λ) を計算する。
    \[ \hat{\pi_{0}} = \frac{\#\{p_{i} > \lambda \}}{m(1-\lambda)} \]
  3. π0(λ) の 3 次曲線を求め、f(λ) とする。そこで π0 推定量は次のように計算する。
    \[ \hat{\pi_{0}} = f(1) \]
  4. qm を計算する。
    \[ q_{m} = \hat{\pi_{0}}p_{m} \]
  5. i=m-1, m-2, ..., 2, 1 に対して、qi を計算する。
    \[ q_{i} = min\left\{ \frac{\hat{\pi_{0}}mp_{i}}{i} , q_{i+1} \right\} \]
  6. qi≤α ならば帰無仮説 Hi を棄却する。
ST法によるq-valueの計算手順を示した図

一方、統計量の分布がわからず、p-valueが求められない場合は、Storey らはブートストラップ法により再サンプリングし、 p-value を推定することを提唱した。論文中では以下の数式で書き表されている。

再サンプリング法によるp-valueを推定する数式

B は再サンプリングの回数、m は標本中の遺伝子数、tjb は b 回目にサンプリングした遺伝子 j の標本から計算される統計量を示す。また、ti は観測データから計算された遺伝子 i の統計量である。

R による q-value の計算

R では qvalue パッケージ中の qvalue 関数を利用して計算する。

library(qvalue)

# p値を用意
p.values <- runif(100)

# q値を計算
q.values <- qvalue(p.values, lambd = seq(0, 0.95, 0.05), pi0.method = "bootstrap")

# q-value
q.value$qvalues[1:10]
 [1] 0.6413621 0.6413621 0.8560256 0.7660976 0.6221614 0.6221614 0.7660976
 [8] 0.8206002 0.8771945 0.8771945

# Pi0の予測値
q.values$pi0
[1] 0.8857143

# Storeyらの論文ではPi0を推定する際に3次曲線を用いましたが、
# 3次曲線以外で近似を撮りたい場合は次のようにします。
q.values <- qvalue(p.values, smooth.df = 4)

References

  • Benjamini Y and Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. Roy. Statist. Soc. Ser. 1995, B57:289-300. JSTOR
  • Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci USA. 2003, 100(16):9440-5. PubMed Abstract