データ科学 多重検定 2 mul%ple test False Discovery Rate 藤博幸
前回の復習 1 多くの検定を繰り返す時には 単純に個々の検定を繰り返すだけでは不十分 5% 有意水準ということは, 1000 回検定を繰り返すと, 50 回くらいは帰無仮説が正しいのに 間違って棄却されてすまうじちがあるということ ex) 1 万個の遺伝子について 正常細胞とガン細胞で それぞれの遺伝子の発現に差があるかどうかを検定
RNA-seq で発現量を 5 回計測 正常細胞 ガン細胞 遺伝子 1 1.4 3.1 遺伝子正常細胞 2 4.5 ガン細胞 1.5 遺伝子 1 1.2 2.7 遺伝子正常細胞 2 4.5 ガン細胞. 1.5 遺伝子 1 1.8 2.3 遺伝子正常細胞 2 4.5 遺伝子ガン細胞. 3000 1.52.1 2.3 遺伝子 1 2.1 2.9 遺伝子 2 4.9 遺伝子. 3000 1.5 2.1 2.3 正常細胞ガン細胞 遺伝子 1 1.0 3.2 遺伝子. 3000 2.1 遺伝子 2 5.4 1.9 2.3 遺伝子 3000 2.1 2.3 遺伝子 3000 2.1 2.3
前回の復習 2 Family Wise Error Rate で制御する 1) Bonferroni 法 2) Holm 法
FWER: FP/(FP + TN) を制御 TN (True Nega*ve) 帰無仮説が正しく 棄却されない FP (False Posi%ve) 帰無仮説が正しいのに棄却される
FWER は厳しすぎる FP ( 帰無仮説が正しいのみ あやまって帰無仮説を棄却すること ; 正常細胞とガン細胞で 遺伝子 A の発現に差がないのに 棄却してしまうこと ) を抑えようとしているが そのため TP ( 帰無仮説が間違っていて 帰無仮説が正しく棄却されること ; 正常細胞とガン細胞で 遺伝子 A の発現に差があり 帰無仮説が棄却されること ) であっても棄却されにくい ----à 発現に差のある遺伝子を検出しにくい ( お金と時間がかかっているのに )
FDR (False Discovery Rate) FP が混じっても良いので TP を増やす ただし FP がどれだけ含めれているか (FDR) を推定して それを新たな基準とする TP をどの程度許可するかの基準を 有意水準の p- 値と区別するため q- 値とよぶ p- 値の分布に仮定を設けて FDR を計算
FDR = FP / (FP + TP) TP = true posi%ve: 対立仮説が正しく 帰無仮説が棄却 FP = false posi%ve: 帰無仮説が正しいのに 帰無仮説が棄却
Benjamini-Hochberg 法 1 m 回の多重検定の場合 (1) p-value を昇順に並べる (p 1 < p 2 < < p m ) (2) i = m とする (3) p i <= α i / m を満たす時 p 1 ~ p i を有意とする そうでなければ i を i 1 にして, 上の条件を確認する i = 1 になっても 条件を満たさない場合は有意なものはないとする FWER の Holm 法が p 1 から順次 p-value が大きい方に向かうのに対し BH 法は p m から出発して, p-value が小さい方に向かうことに注意
Benjamini-Hochberg 法 2 p 1 < p 2 < < p i < p m とする i 番目 (p i ) を検討する α を false posi%ve を含む割合とする. m α は m 回の検定中 false posi%ve の期待回数 i / m は, m 回の検定中 i 個の占める割合 すると, m α i / m = α i は i 回の検定の中で false posi%ve を起こす回数 (= 棄却されたのに本当は帰無仮説が正しいもの ) の期待値 FP の上限として α i を考える
Benjamini-Hochberg 法 3 今 p i 以下の p-value で棄却するすると posi%ve となる検定が i 個ある この i 個には true posi%ve も false posi%ve も含まれているとすると FP + TP = i FDR = FP/( FP + TP) = FP/i < α i / i = α このようにして FDR を α 以下になるよう制御できる
Benjamini-Hochberg 法 4 p i <= α i / m を満たす時 p 1 ~ p i を有意とする この式を書き直すと p i m / i <= α を満たす時 p 1 ~ p i を有意とする q i = p i m / i を p-value に対して q-value とよぶ
Benjamini-Hochberg 法 5 実際に行ってみよう p-value Gene 1 0.21 Gene 2 0.001 Gene 3 0.1 Gene 4 0.06 Gene 5 0.005 p- 値のソート p-value Gene 2 0.001 Gene 5 0.005 Gene 4 0.06 Gene 3 0.1 Gene 1 0.21 h_ps://www.slideshare.net/an%plas%cs/dna-21259335 より
Benjamini-Hochberg 法 6 p-value q-value q- 値への変換 Gene 2 0.001 Gene 5 0.005 Gene 4 0.06 Gene 3 0.1 Gene 1 0.21 0.001 x (5/1) 0.005 0.005 x (5/2) 0.00125 0.06 x (5/3) 0.1 0.1 x (5/4) 0.125 0.21 x 5/5 0.21
Benjamini-Hochberg 法 7 FDR の閾値を 0.05 として p-value の大きいものから順番に検討 p-value q-value Gene 2 0.001 Gene 5 0.005 Gene 4 0.06 Gene 3 0.1 Gene 1 0.21 0.001 x (5/1) 0.005 0.005 x (5/2) 0.00125 0.06 x (5/3) 0.1 0.1 x (5/4) 0.125 0.21 x 5/5 0.21 < 0.05 < 0.05 > 0.05 > 0.05 > 0.05 Gene 2 と Gene 5 が有意
R でやってみよう 1 Gene1~Gene5 の p-value をベクトルとして表現 > pv <- c(0.21, 0.001, 0.1, 0.06, 0.005) pv を昇順にソート > spv <- sort(pv) > spv [1] 0.001 0.005 0.060 0.100 0.210 q-value を記憶させる空ベクトルを作成 > qv <- c() for 文で q-value に変換 > for (i in 1:length(spv)) qv <- c(qv, spv[i]*length(spv)/i) > qv [1] 0.0050 0.0125 0.1000 0.1250 0.2100 前ページの表と同じ結果が得られたことを確認
R でやってみよう 2 p.adjust 関数を使ってみる > p.adjust(pv, method="bh", n=length(pv)) [1] 0.2100 0.0050 0.1250 0.1000 0.0125 > p.adjust(pv, method="bh", n=length(pv)) < 0.05 [1] FALSE TRUE FALSE FALSE TRUE 遺伝子の並び順のままで検定を行える
Benjamini-Hochberg 法の問題点 1 Benjamini-Hochberg 法のキモは q-value の計算で q i = p i m / i を p-value に対して q-value とよぶ p i を基準とした時 その時の false posi%ve の数が p i m と仮定している点にある p i が一様分布していればこの仮定は正しい 帰無仮説が全て正しい場合 (ex. 10000 個全ての遺伝子で正常細胞とガン細胞における発現量に差はない ) p-value は一様分布する
帰無仮説が全て正しい場合における p-value の一様分布の確認 前回の多重検定での false posi%ve 生成と同じシミュレーションを実施 > N <- 10000 # サンプル数を 10000 > p.values <- double(n) # p-value を格納するベクトル > for(i in 1:N) { + x <- rnorm(10, mean=0, sd=1) + y <- rnorm(10, mean=0, sd=1) + p.values[i] <- t.test(x,y)$p.value + } # 同じ正規分布から 10 こずつサンプルして平均値の差の検定 # 同じ分布からのサンプルなので差はないはず > hist(p.values) # p-value のヒストグラム作成
p-value がほぼ一様分布に従っていることがわかる
Benjamini-Hochberg 法の問題点 2 h_ps://sites.google.com/site/scriptocioinforma%cs/maikuroarei-guan-xi/fdr-zhi-yu-r
Benjamini-Hochberg 法の問題点 3 実際には 差があるものが含まれるので p-value は 0 に近い方に偏った分布となる > N <- 5000 > rm1 <- runif(n)*5 > rm2 <- runif(n)*5 > rs1 <- runif(n)*2 > rs2 <- runif(n)*2 > p.values <- double(n*2) > # 5000 個は同じ分布からサンプリング > for (i in 1:N) { + x <- rnorm(10, mean=0, sd=1) + y <- rnorm(10, mean=0, sd=1) + p.values[i] <- t.test(x,y)$p.value + } > # 残り 5000 個はランダムに生成した平均と標準偏差を持つ分布からサンプリング > for (i in 1:N) { + x <- rnorm(10, mean=rm1[i], sd=rs1[i]) + y <- rnorm(10, mean=rm2[i], sd=rs2[i]) + p.values[n+i] <- t.test(x,y)$p.value + } > hist(p.values)
p-value の一様性は成立していないこの点については Q-value 法や Local FDR 法などの改良があるが今回は説明を省く
帰無仮説の分布と対立仮説の分布の混合分布帰無仮説と対立仮説の密度比が各々 π 0 :1 - π 0 の比だったとする h_ps://sites.google.com/site/scriptocioinforma%cs/maikuroarei-guan-xi/fdr-zhi-yu-r h_ps://www.slideshare.net/yuifu/fdr-kashiwar-3 の記述だともう少し複雑
ヒストグラムのどの場所で高さをとるかという問題を考えたのが Q-value Q-value では 高さをとる場所を γ γ から推測される π 0 の値を π 0 (γ) としてプロットし このデータを自然スプライン関数で回帰 このスプライン関数を γ=0 まで外推した値を π 0 の予測値とする p-value のアナロジーとして使われている小文字の q-value とは違って これは Q value という手法名である点に注意 Q value はこの π 0 を BH 法の q-value にかける それ以外の計算手順は BH 法と同じなので 先にこの π 0 の予測値をかけておいた p -value を BH 法に適用すれば Q-vakue を求める事ができる R のパッケージ qvalue をインストールすると Q-value を計算できる
h_p://strimmerlab.org/notes/fdr.html
h_ps://www.slideshare.net/an%plas%cs/dna-21259335 h_ps://www.slideshare.net/yuifu/fdr-kashiwar-3 h_ps://sites.google.com/site/scriptocioinforma%cs/maikuroarei-guan-xi/fdr-zhi-yu-r h_p://d.hatena.ne.jp/hoxo_m/20120413/p1 h_p://www.med.osaka-u.ac.jp/pub/kid/clinicaljournalclub1.html