マルチプルテスティングとFDR〜オミックス統計学入門2014
- 2014年度講義資料
- FDRとマルチプルテスティング
- 90分で1-2回分相当
- Rmdファイルです。html化、epub化できます(やり方はこちら)
- html化、epub化が面倒くさければ、kindleで1米ドルでも(こちら)
マルチプルテスティングとFDR?オミックス統計学2014_4? 統計遺伝学
- 作者: ryamada
- 発売日: 2014/02/07
- メディア: Kindle版
- この商品を含むブログ (1件) を見る
# マルチプルテスティングとFalse Discovery Rate # 検定とp値と、どれくらい信じるか(帰無仮説の棄却) ## 一様分布とp値 一様分布というのは、どんな値も同じ確率で出現する分布である。 横軸が0から1、縦軸も0から1の正方形になっている。 x軸の値によらず、出現確率が同じ高さである。 検定で用いられるp値はこのような性質を持っている。 小さいp値も大きいp値も中くらいのp値も同じ確率で出現する。 ```{r} x <- seq(from=0,to=1,length=10000) d <- dunif(x) plot(x,d,type="l",ylim=c(0,1.1),col=2) abline(v=c(0,1)) abline(h=0) ``` たとえば0.05より小さいp値は次の図の赤い部分であって、全体の0.05を占める。 ```{r} x <- seq(from=0,to=1,length=10000) col <- rep(1,length(x)) col[which(x<0.05)] <- 2 plot(x,d,type="h",ylim=c(0,1.1),col=col) abline(v=c(0,1)) abline(h=0) ``` 従って、p値が0.05であるということは、「p値が0.05以下になる確率が0.05である」ということを意味している。 これを別の方法で示す。 ```{r} p <- punif(x) plot(x,p,type="p",cex=0.01,ylim=c(0,1.1),col=col) abline(v=c(0,0.05,1)) abline(h=c(0,0.05,1)) ``` p値を小さい方から累積してプロットすると、グラフは$y=x$の直線となり、原点と(1,1)の２点を結んだものとなる。 p値が0.05であるということは、この直線の左側0.05の部分を登ったところであることを示している。 ## 帰無仮説を棄却する ある試料を一斉に集め、適当に分配して、Ａ・Ｂ２社に外注したところ、Ａ社とＢ社とのそれぞれからの納品濃度の分布が次のようになっていたとする。 Ａ社とＢ社とで納品濃度に違いがないと考えるかどうかを試してみる。 ２社の納品濃度の分布が同じであると仮定してデータを作成し、t検定をしてみる。 同じことを仮定しているので、帰無仮説は棄却されないはずである。 p値は0.05より小さくならないだろうと予想されると言うこともできる。 ```{r} ma <- 120 sa <- 10 mb <- 120 sb <- 10 na <- 20 nb <- 15 A <- rnorm(na,ma,sa) B <- rnorm(nb,mb,sb) h <- hist(c(A,B),plot=FALSE) hist(A,breaks=h$breaks,ylim=h$ylim,density=10) hist(B,breaks=h$breaks,density=7,col=2,add=TRUE) t.test(A,B)$p.value ``` 小さいp値が出ないだろうという予想は、高い確率で裏切られませんが、 低い確率では裏切られる。 それを試してみる。 ```{r} n.iter <- 10000 ps <- rep(0,n.iter) for(i in 1:n.iter){ A <- rnorm(na,ma,sa) B <- rnorm(nb,mb,sb) ps[i] <- t.test(A,B)$p.value } hist(ps) plot(sort(ps)) ``` たしかに、p値の分布は一様分布らしくなっており、0.05より小さいp値は確率0.05で出現する。 ## たくさんの検定〜GWAS〜 ゲノムワイドに50万個のSNPと誕生月が偶数か奇数かとの関係を調べるとする。 無関係な50万の検定が繰り返される。 ```{r} n.snp <- 500000 p <- runif(n.snp) plot(p,cex=0.1) ``` 点の数が多すぎて何のことだかわからないので、対数表示にしてやる。 こうすると、いわゆるGWASのマンハッタンプロットになる ```{r} plot(-log10(p),cex=0.1,ylim=c(0,10)) ``` 一番小さいp値はいくつだろうか。 プロット上で強調しておく。 ```{r} min(p) plot(-log10(p),cex=0.1,ylim=c(0,10)) minp <- which(p==min(p)) points(minp,-log10(p[minp]),col=2,pch=20) ``` 今、p値はすべて帰無仮説が成り立つ状態で発生させたのであるから、 この小さいp値を見ても、 > 帰無仮説が成り立っているときの、ごく普通の値であって、珍しくはない と感じるべきである。 ではどれくらい小さいp値を見たら、帰無仮説を疑うほど、小さいことになるだろか。 珍しさは、「何度も繰り返してもめったに起きない」ことで決まるので、50万SNPの検定を何度も繰り返し、そのもっとも小さいp値がどのような分布をするかを調べてみればよい。 その分布に照らしてめったに起きないほど小さければ、そのとき初めて帰無仮説を疑えばよい。 ```{r} n.iter <- 1000 minp <- apply(matrix(runif(n.iter*n.snp),nrow=n.iter),1,min) hist(minp) ``` 図に示すように、大きい方に裾を引いた分布である。 ```{r} summary(minp) quantile(minp,c(0.01,0.05,0.1)) ``` 50万個のSNPによるGWASを1000回実施したときの各GWASでの最小p値の分布である。 1000回のGWASのそれぞれは、全く遺伝因子が関与していないと思われるようなフェノタイプに対して実施してある。 小さい方から１パーセンタイル、５パーセンタイル、１０パーセンタイルを示したが、この値に比べて、小さい最小p値を観測したときにはじめて、そのp値をもたらしたSNPとフェノタイプとの間では帰無仮説が成り立たないのではないかと考えることとなる。 この数値は大まかには$0.01/500000=2\times 10^{-8}$でとして算出できる。 珍しさ0.01を検定の個数50万で割った値である。 このようにして、珍しさの基準を決める方法をボンフェロニ補正法と呼ぶ。 最も小さいp値の大小を問題にしたが、すべてのp値についても情報を得ることにする。 たくさんのp値のすべてが帰無仮説に由来しているのかどうかを図で示すには、小さい順にならべて、$y=x$の直線に乗るかどうかでわかるのだった。 実際にやってみる。 横軸が0から1となるように工夫しておく。 ```{r} plot(ppoints(n.snp,a=0),sort(p),cex=0.01) ``` ただし、これでは、p値が小さい部分が、原点付近の狭い範囲に固まっていて、よくわからないので、その部分を拡大してみるために、対数をとれば ```{r} plot(log10(ppoints(n.snp,a=0)),log10(sort(p)),cex=0.01) ``` また、場合によって(論文ではそちらの方が多いが)、p値ではなく、統計量($\chi^2$統計量)の観測値とその期待値とで、直線になっているかを示す場合もある。 ```{r} plot(qchisq(ppoints(n.snp,a=0),1),qchisq(sort(p),1),cex=0.01) ``` # 対立仮説と推定 ## たくさんの治療法報告 ２種類の治療薬の効果を比較した研究が次から次へと発表される。 小さなp値とともに有効だ、とする報告もあれば、あと一歩で帰無仮説が棄却できるレベルだったけれど、ちょっと足らなかったが、でも、惜しかった、という感じで報告されることもある。 単純に考えるためにp値が0.05であったら、必ず発表され、0.05から0.07の間のp値の場合には、10%が発表されるとする。 そんな報告をたくさん集めると、p値はどんな分布になるか。 ```{r} n.t <- 100000 p <- runif(n.t) ok <- which(p<0.05) sub.ok <-which(p>=0.05 & p<0.07) published.sub.ok <- sample(sub.ok,length(sub.ok)*0.1) hist(p[c(ok,published.sub.ok)]) ``` 0.05より小さい部分と0.05から0.07の部分とで段差はできるが、それぞれの範囲で一様分布になる。 ## 検出力 実際の治療法の報告のp値の分布がこうなっているかというとそんなことはない。 カットオフの閾値0.05を下回ったあたりにたくさんの報告が集中する。 それはなぜかというと、 > 効果があるならこれくらい、という見込みの下、必要な標本数を計算して、スタディをデザインする からである。 一般的な標本数の計算をしてみる。 治療法Aの有効率が50%、治療法Bのそれが60%であると見積もられるとき、p値の閾値として0.05で両者に違いがないという帰無仮説を棄却することにして、デザインする。 きちんとp値が0.05を下回る確率が80%となるには、検体数を2群に何人ずつとればよいかを考える。 ```{r} pw.out <- power.prop.test(p1 = .50, p2 = .60, power = .80) pw.out$n n. <- round(pw.out$n) n. ``` で得られる数値が、２群に同数の検体数を取る場合に確保するべき検体数である。 実際、この条件で、このスタディを繰り返してp値がどうなるかを見てみる。 ```{r} n.iter <- 1000 ps <- rep(0,n.iter) pa <- 0.5 pb <- 0.6 for(i in 1:n.iter){ A <- sample(0:1,n.,replace=TRUE,prob=c(1-pa,pa)) B <- sample(0:1,n.,replace=TRUE,prob=c(1-pb,pb)) tbl <- matrix(c(table(A),table(B)),ncol=2) ps[i] <- chisq.test(tbl,correct=FALSE)$p.value } hist(ps) ``` 実際p値が0.05を下回った割合は ```{r} ok <- which(ps<0.05) length(ok)/n.iter ``` と、80% 程度になっている。 では、このようなスタディがたくさんあって、p値が0.05を下回ったら報告され、 若干、上回ったときには、その10%が「惜しい」として報告されるとすると、p値はどんな分布になるのだろうか。 ```{r} sub.ok <- which(ps>=0.05 & ps<0.07) published.sub.ok <- sample(sub.ok,length(sub.ok)*0.1) hist(ps[c(ok,published.sub.ok)]) ``` 0.05ぎりぎりの報告は少なく、かなり小さいp値の報告が目白押しになるはずであることがわかる。 実際の報告のp値の分布は、もっと、0.05ぎりぎりのものが多い印象がないだろうか。 ## 色々な治療効果のスタディの混成 ここで示した分布は、2群の成功率が0.5,0.6と「きちんとした差がある」場合のp値の分布であり、 実際の報告の集まりの中には、このように違いのはっきりしたスタディの報告もあれば、まったく2群に違いがないのにもかかわらず、確率0.05で生じてしまう小さなp値のせいで報告されるものも含まれること、また、若干の差はあるが、2群の成功率が0.5と0.55だったり、0.5と0.51だったりするものも含まれるからである。 では、A群の成功率は0.5であって、B群の成功率は0.5の場合と0.6の場合とが半々だとどういうことになるだろうか。 ```{r} n.iter <- 1000 ps <- rep(0,n.iter) pa <- 0.5 pb <- 0.6 for(i in 1:n.iter){ A <- sample(0:1,n.,replace=TRUE,prob=c(1-pa,pa)) if(i < n.iter/2){ pb. <- pa }else{ pb. <- pb } B <- sample(0:1,n.,replace=TRUE,prob=c(1-pb.,pb.)) tbl <- matrix(c(table(A),table(B)),ncol=2) ps[i] <- chisq.test(tbl,correct=FALSE)$p.value } hist(ps) ``` p値のうち半分は、上で描いた0寄りの分布であり、残りの半分は0から1の一様分布であって、それを併せた分布となる。 このようなスタディから、報告がなされるとき、p値の分布はどうなるだろうか。 今回は、p値が0.05の場合にのみ報告し、それ以外は厳密に報告しないものとする。 ```{r} ok <- which(ps<0.05) ok.false <- which(ok<(n.iter/2)) ok.true <- which(!(ok<(n.iter/2))) h <- hist(ps[ok],plot=FALSE) hist(ps[ok[ok.true]],breaks=h$breaks,ylim=h$ylim,density=20) hist(ps[ok[ok.false]],breaks=h$breaks,ylim=h$ylim,density=15,col=2,add=TRUE) ``` たくさんの黒(A,Bの2群に違いがあってp値が0.05より小さくなったスタディ)があり、少数の赤(A,Bの2群に違いがないのにp値が0.05より小さくなったスタディ)があることがわかる。 黒のバーはp値が小さい方で長いが赤のバーは高さが一様である。 そのことから、p値が小さいスタディがあれば、それは本当に2群間に違いがあることを反映している割合が高いし、p値が0.05に近ければ、本当は2群に違いがないのに報告されてしまっている可能性も低くないことがわかる。 では、2群に違いがある(有用な治療法である)と判断したもののうち、どれくらいの割合が、「本当に有用」なのだろうか。 ```{r} length(ok.true)/length(ok) ``` そして、本当は有用でないのに、有用であると誤って判断されたのはどれくらいの割合だろうか。 ```{r} length(ok.false)/length(ok) ``` ## 感度・特異度・事前確率・Positive Predictive Value・Negative Predictive Value・False Discovery Rate これが、診断のための検査であったとすると、どう考えればよいだろうか。 まず、2群に違いがあると仮定したスタディ、違いがないと仮定したスタディが50%ずつだったので、検査の前に真に「有効〜病気である」ことの事前確率は50%。 ついで、本当に「差がある〜病気である」場合に、p値が0.05未満になる確率は80%であるから、これは、検査で言えば、感度が80%。 また、本当は「差がない〜病気ではない」場合にp値が0.05未満になる確率は5%であるから、特異度は95%。 このとき、p値が0.05未満になったときに、本当に病気である確率が、上で計算した ```{r} length(ok.true)/length(ok) ``` であり、これはPPVに相当し、 ```{r} length(ok.false)/length(ok) ``` は本当は病気でないのに病気であると判断されてしまった確率で、これはFalse Discovery Rate(FDR)である。 # FDR法 とq値 すべての検定が帰無仮説を満足することを前提にしてp値を補正する例をGWASで見た。 かなりの割合で対立仮説が成り立つ場合に、検定をすると、PPVが問題になり、FDR=1-PPVも気になることを、有用な治療法の研究を例として見た。 たくさんの検定を実施し、それらのうちのいくらかが対立仮説を満足するような状況の典型例は、2群(疾患群vs.健康群、正常組織 vs.異常組織)の間での遺伝子発現量を全遺伝子について比較するトランスクリプトーム解析である。 ２つの群は明らかに異なっており、その2群での遺伝子発現量が違っているのも確実であり、興味があるのは、「どの遺伝子がどれくらい」2群間で違うのか、である。 そのようなとき、すべての遺伝子に関する2群間の違いついて帰無仮説検定をし、そのすべての遺伝子に関して帰無仮説が成り立つことを想定したうえでマルチプルテスティング補正をするのは適当ではない。 ある閾値で違いがありそうな遺伝子を選び出したときに、選ばれた遺伝子にどれくらい本物があり、どれくらい偽物があるかの比率が、納得の行くものであることが大事であり、その比率がわかっていること、わかった上でその結果を解釈することが大事である。 このような場合に、FDRを制御して全遺伝子が「本当に2群間で違う」かどうかを評価する手法を「FDR法」と言う。 FDR法では、たくさんのp値が得られたときに、そのpの値をあるルールで変換して0からの1の値に変える。 その値を補正後p値とかFDR法によるq値と言う。 p値が小さい方が「本物」らしいわけだが、補正後p値も小さい方が「本物らしい。また、2つのp値があったら、その大小関係は、補正後p値でも同じである。 FDR法にはいくつかの方法があるが、まずはBH法と呼ばれる方法でやってみる。 p値と補正後p値との関係は以下のようになる。 値が変わっているので、一直線ではなくなっている。 また、pの値よりも補正後p値の方が大きくなっていることも見て取れる。 ```{r} p.corr <- p.adjust(ps,method="BH") plot(ps,p.corr,cex=0.1) ``` このようにして計算した補正後p値の度数分布を補正前のそれと比べてみる。 全体的に少し大きな値に補正するということが度数分布ではどのように見えるかがわかる。 ```{r} hist(ps,density=20) hist(p.corr,col=2,density=15,add=TRUE) ``` このような補正をした結果、PPVとFDR=1-PPVがいくつになったかを調べてみる。 ```{r} ok.corr <- which(p.corr<0.05) ok.corr.false <- which(ok.corr<(n.iter/2)) ok.corr.true <- which(!(ok.corr<(n.iter/2))) length(ok.corr.true)/length(ok.corr) length(ok.corr.false)/length(ok.corr) # 補正前の値は length(ok.false)/length(ok) ``` FDRが補正により小さくなっていることがわかる。 そうは言うものの、FDRの値がきっちり0.05になるわけではないこともある。 では、きっちりFDRがコントロールできていないのならば、どんな基準で補正がなされているのだろうか。 すべて帰無仮説が成り立っているときにはどうなるかを見る。 ```{r} p.null <- runif(1000) p.corr.null <- p.adjust(p.null,method="BH") plot(p.null,p.corr.null) hist(p.corr.null) ``` ほとんどすべての補正値が1付近になっており、補正後に0.05になるようなものはないことがわかる。 このことから、 >たくさんのp値が帰無仮説に由来していると考えられる場合には、そのすべてが「ありきたり〜補正後p値が1」となるような補正法であることがわかる。 考え方としてはこのような補正法であるが、さまざまな改良がなされている。 その基本とは、補正前のp値の分布が、すべて帰無仮説の場合と比べてどのようにずれているかを評価し、それに応じて、どのくらいの割合で「対立仮説が真である」かを見積もり、補正後のp値(q値)が、「偽の陽性率」とみなせるように補正するものである。 その一つを以下に示す。 ```{r} library(fdrtool) fdr <- fdrtool(ps, statistic="pvalue",verbose=FALSE,plot=FALSE) plot(ps,fdr$qval) # estimated Fdr values hist(fdr$qval) ```
