六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

多重比較補正のはなし

最近になって、データ分析界隈で多重比較補正が話題に上ることが増えていると聞きまして。一方で、僕自身も何を隠そう研究者時代の専門分野が長年多重比較補正の問題に悩まされてきた分野だったこともあって、かなり若い頃から多重比較補正については色々勉強したり実践したり検討したりしてきたものでした。


ちなみに下記のリンクはその分野で広く使われている多重比較補正の方法論に重大な瑕疵があるのではないかと指摘した2年前の論文で、曰く「理論上は偽陽性(false positives)5%で済むはずのものが実際には偽陽性が最大70%に達する可能性がある」とのこと*1

事によっては15年間に渡る40000件の研究が実は偽陽性まみれだったという可能性もあるという話で、多重比較補正って怖いんだなぁとつくづく思う次第です。


閑話休題。このように多重比較補正というのは「偽陽性」という統計データ分析においては厄介な問題と密接に関わるものでもあります。そんなわけで、この記事では昔勉強した内容を思い出しながら多重比較補正にまつわるあれこれを書き綴ってみようとかと思います。特に学術的にきっちりとした議論をするつもりはないので、内容としてはほんのさわり程度に留めます。


そうそう、今回の話題も僕自身頑張って勉強していた頃から既に6年以上が経って曖昧な記憶を思い出しながら書いている有様なので、間違っているところが多々あるかもしれません。いつも通りご意見ご指摘などあればどしどしお寄せいただけると有難いですm(_ _)m


多重比較における偽陽性とは何か


偽陽性(false positive)とは、平たく言えば「本来陽性ではないのに誤って陽性とみなされてしまうこと(みなされてしまったもの)」のことです。詳しくは例えば英語版WikipediaMultiple comparisons problem - Wikipediaの項あたりをお読みいただきたいのですが、以下の表のVに当たるものです。なお以下の表はWikipediaから引用したもので、検定の数がm個ある状況で検定仮説もH_1, H_2, \cdots, H_mと並列して存在すると仮定した場合のものだそうです。

帰無仮説が正しい(H_0) 対立仮説が正しい(H_A) 合計
検定が有意である V S R
検定が有意ではない U T m-R
合計 m_0 m - m_0  m
  •  m: 仮説検定の総数(多重比較の個数)
  •  m_{0}: 真に帰無仮説に従う検定数(未知パラメータ)
  •  m-m_{0}: 真に対立仮説に従う検定数
  •  V: 偽陽性の数(第一種の過誤、"false discoveries"とも称する)
  •  S: 真陽性の数("true discoveries"とも称する)
  •  T: 偽陰性の数(第二種の過誤)
  •  U: 真陰性の数
  •  R = V + S: 帰無仮説が棄却された検定数(true / negative問わず"discoveries"とも称する)

この表を見ればすぐ分かるように、V(第一種の過誤)とT(第二種の過誤)は「正しくない*2検定結果の数」であり、正しい多重検定結果を得たければ出来るだけ減らしたいものであるというわけです。なお、この辺の詳しい話は、例えば東大出版会『自然科学の統計学』(通称青本)の第6章を読むと出てきます。


このような考え方が出てくる背景としては、言うまでもなく統計的手法が原則として「無作為抽出されたサンプルに対する分析」であるということであり、そこを起点としてネイマンとピアソンの「第一種の過誤と第二種の過誤」という概念が生まれたという歴史的経緯にも拠っています。個人的な理解を書いてしまうと「サンプルの抽出自体に既に原理的に解決不能なばらつきが紛れ込んでしまうのがこの自然の森羅万象というものなので、そのばらつきに振り回され得ることを加味した上で出来るだけ間違った結論を排除すべきであり、間違った結論に達してしまう確率の容認可能なレベルはあらかじめ決めておく必要がある」ということなのかなと。


しかしながら、これは裏を返すと例えば有意水準5%と定めた場合「完全なる偶然からAとBという2つの全く同一の分布から生成されたサンプル同士を比較してもこれを100回繰り返せば5% = 5回は差があるという間違った結論に至り得る」ということを意味します。100回のうち5回なら、容認するという人もいるでしょう。けれども、これが100000回のうち5000回となれば。。。とても受け容れられないという人も多いのではないでしょうか。


そこで、「様々な事情から同じような比較を何千回何万回と繰り返した(多重比較)場合であっても適切に偽陽性を減らす」ことを目的として、様々な方法が考えられてきたわけです。それらを総称して、「多重比較補正」と一般に呼びます。


偽陽性の簡単な例


というのは理論的な解説で、これだけ書いてもピンと来ない人が多いかもしれないので簡単な例を挙げます。まず、以下のように[-5, 5]の一様乱数から生成した20サンプルを1組として、適当に10000組生成します。次に、これを適当にランダムに5000組ずつに二手に分けます。そして、二手に分けたサンプル5000組から1組ずつ取ってきてはt検定するというのを5000ペア分繰り返して、p値を保存します。最後に、p値5000個のヒストグラムを描画し、p < 0.05を満たす個数をカウントします。


以上の操作を実行したのがこちらのRコードです。

> set.seed(101)
> d <- matrix(runif(200000, -5, 5), ncol = 20, nrow = 10000, byrow = T)
> set.seed(102)
> idx <- sample(10000, 5000, replace=F)
> d1 <- d[idx, ]
> d2 <- d[-idx, ]
> p_d <- c()
> for (i in 1:5000){
+      tmp <- t.test(d1[i, ], d2[i, ])
+      p_d <- c(p_d, tmp$p.value)
+ }
> hist(p_d)
> length(which(p_d < 0.05))
[1] 227

f:id:TJO:20180611000445p:plain

全く同じ分布から生成したサンプル同士でt検定しているので、本来ならばt検定の真の結果は「有意差なし」になるはずです。しかしながら、結果はご覧の通り。。。5000組中227組が「有意差あり」となってしまっています。これはおかしな話であり、同時にこれこそが「偽陽性」そのものであると言えます。実際、5000組の4.54%が偽陽性になってしまっており、有意水準5%という当初の閾値をある程度反映していることも分かります。


こんなことが本当に起こり得るのか?というと、特に実験科学においては割と頻繁に起こり得る話です。例えば、先述した僕の前の専門分野であるヒト脳機能画像研究の世界では、ヒトの脳を機能的MRIによっておよそ2万くらいの画素(voxel)に分け*3、それぞれのvoxelごとに線形モデルを利用して脳の活動量が「課題遂行時>安静時」となる部位を探し出すということをしていました。この場合2万組のvoxelにおいて特定条件間での活動量比較を同時に行うため、上記の例のように完全な偶然であっても例えば1000個ぐらいのvoxelは有意差を示してしまう可能性があります。バイオインフォマティクス分野でも同様の問題があるらしく、例えば遺伝子発現解析などでも多重比較補正を使うべきだという議論があったりするようです。


もちろん基礎研究分野に限らず、例えば大規模なwebサービスでクリエイティブを多数のページで一括変更した場合にPV / CVがどれくらい伸びたか?みたいな比較を行うのも、個々のページにつき1つ以上A/Bテストで比較するということを行うならば立派な多重比較補正の対象になり得ます。時々web開発絡みで多重比較補正の話題が出てくるのは、恐らくその辺の事情によるのでしょう。


ところで。ここで例えばp値の閾値を0.05ではなく5000で割ったp < 0.05/5000 = 1e-05に変えてみると、こうなります。

> length(which(p_d < 0.05/5000))
[1] 0

偽陽性ゼロとなりました。これは後述するBonferroni補正そのものの例ですが、このようにして検定の方式自体に対して様々な工夫をすることによって偽陽性を出来るだけ少なくするというのが、多重比較補正の基本的なコンセプトです。


偽陽性を補正する方法


では、どうやって偽陽性を補正するかというその具体的な方法ですが。Wikipedia記事の続きを見れば分かるように、考え方は色々あります。続きの記述に従えば、大まかに

のどちらをコントロールするか、によって分けられます。以下詳細はWikipedia記事のページ内アンカーにリンクを張ってありますので、詳細はリンク先をご確認ください。


1つ目のFWERですが、上記の表に従えば FWER = Pr(V \geq 1) = 1 - Pr(V = 0)と表せます。これを抑えるということは、直接偽陽性の割合自体を減らしてしまえという考え方だとも言えます。FWERに関するWikipedia記事によれば以下の方法論があります。

これ以外にもResampling procedureすなわちブートストラップ法によって閾値を決定する方法も紹介されていますが、個人的にはこちらの方法は特に多重比較の件数が多い場合には計算量が大き過ぎてキツい気がするので、ちょっとどうかなぁという気がしなくもないです(経験談*4)。


2つ目のFDRですが、これは FDR = E[ Q ], ~ Q = V / R = V / (V + S) と表せます。これを抑えるということは、得られた陽性数における偽陽性の割合を低くするという考え方だとも言えます。FDRに関するWikipedia記事によれば以下の方法論があります。


なお、Wikipedia記事ではFalse coverage rate (FCR)を抑えるという考え方についても触れられていますが、こちらは僕自身があまり詳しくないことと、その実装をあまり見かけないということもあるので、ここでは割愛させていただきます。


余談ながら、研究者時代に僕が多用していたのはBonferroni補正とFDR control(つまりBH補正)でした。理由は簡単で、大体どのデータ解析ソフトウェアにもその2つは実装されていたからですorz ただし当時FDR controlを使っていた経験から言うと、ステップワイズ法の挙動はp値の分布にかなり依存するのではないかと思っています。と言うのも、例えば上記のリンク先にもあるBH補正の定義を見ると

  • m組ある検定のp値を小さい順に並べて P_{(1)}, P_{(2)}, \cdots, P_{(m)}と表す。対応する帰無仮説 H_{(1)}, H_{(2)}, \cdots, H_{(m)}とする
  • FDR controlレベルとして \alphaを置き(普通は0.05)、 P_{(k)} \leq \frac{k}{m} \alphaを満たす最大の kを探す
  •  H_{(1)}, H_{(2)}, \cdots, H_{(k)}を棄却し、 P_{(k)}をFDR  \alpha閾値に対応するp値とみなす

となっていて、 P_{(k)}の値の選ばれ方は P_{(1)}, P_{(2)}, \cdots, P_{(m)}の分布に強く依存するからです。HolmとHochbergの方法も同様のステップワイズ法を使っているので、恐らくp値の分布次第という結果になるのだと思われます。


Rで多重比較補正を実践するには


Rではデフォルトで同梱されている{stats}パッケージ*5にp.adjustというメソッドとして実装されていて、もちろん中身を見ることもできます。こんな感じです。上記のWikipediaによる説明よりはこちらを読んだ方が分かりやすいという人も少なくないのではないかと思います。

function (p, method = p.adjust.methods, n = length(p)) 
{
    method <- match.arg(method)
    if (method == "fdr") 
        method <- "BH"
    nm <- names(p)
    p <- as.numeric(p)
    p0 <- setNames(p, nm)
    if (all(nna <- !is.na(p))) 
        nna <- TRUE
    p <- p[nna]
    lp <- length(p)
    stopifnot(n >= lp)
    if (n <= 1) 
        return(p0)
    if (n == 2 && method == "hommel") 
        method <- "hochberg"
    p0[nna] <- switch(method, bonferroni = pmin(1, n * p), holm = {
        i <- seq_len(lp)
        o <- order(p)
        ro <- order(o)
        pmin(1, cummax((n - i + 1L) * p[o]))[ro]
    }, hommel = {
        if (n > lp) p <- c(p, rep.int(1, n - lp))
        i <- seq_len(n)
        o <- order(p)
        p <- p[o]
        ro <- order(o)
        q <- pa <- rep.int(min(n * p/i), n)
        for (j in (n - 1):2) {
            ij <- seq_len(n - j + 1)
            i2 <- (n - j + 2):n
            q1 <- min(j * p[i2]/(2:j))
            q[ij] <- pmin(j * p[ij], q1)
            q[i2] <- q[n - j + 1]
            pa <- pmax(pa, q)
        }
        pmax(pa, p)[if (lp < n) ro[1:lp] else ro]
    }, hochberg = {
        i <- lp:1L
        o <- order(p, decreasing = TRUE)
        ro <- order(o)
        pmin(1, cummin((n - i + 1L) * p[o]))[ro]
    }, BH = {
        i <- lp:1L
        o <- order(p, decreasing = TRUE)
        ro <- order(o)
        pmin(1, cummin(n/i * p[o]))[ro]
    }, BY = {
        i <- lp:1L
        o <- order(p, decreasing = TRUE)
        ro <- order(o)
        q <- sum(1L/(1L:n))
        pmin(1, cummin(q * n/i * p[o]))[ro]
    }, none = p)
    p0
}
<bytecode: 0x1030b42d0>
<environment: namespace:stats>


このp.adjustメソッドを適当にデモ用データセットに使って見ましょう。まず、最初のデモ用データセットは以下のコードで作ることにします。[10, 15]と[10.5, 15.5]という僅かにずれた2つの一様分布から10000組ずつ20サンプルのデータセットを生成しています。クソコードなのはお見逃し下さい(泣)。

> set.seed(201)
> x1 <- matrix(runif(200000, 10, 15), ncol = 20, nrow = 10000, byrow = T)
> set.seed(202)
> x2 <- matrix(runif(200000, 10.5, 15.5), ncol = 20, nrow = 10000, byrow = T)
> p <- c()
> for (i in 1:10000){
+     tmp <- t.test(x1[i, ], x2[i, ])
+     p <- c(p, tmp$p.value)
+ }
> hist(p)
> length(which(p < 0.05))
[1] 1846

f:id:TJO:20180611223055p:plain

わずかに重なっている分、p値が低い方に分布が偏っていることが分かります。実際、補正なしだとp < 0.05となる組は1846個と同じ分布同士で比べた場合よりもかなり多いです。


p.adjustメソッドはmethod引数に"holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"を選ぶことが出来ます。そこで、それぞれの手法を選んだ上で多重比較補正後にp < 0.05となる組が幾つあるかを算出してみます。ただし実装的にBHとfdrは同じであり、noneは文字通り無補正なので、fdrとnoneは省略します。

> p_holm <- p.adjust(p, method = 'holm')
> length(which(p_holm < 0.05))
[1] 0
> p_hochberg <- p.adjust(p, method = 'hochberg')
> length(which(p_hochberg < 0.05))
[1] 0
> p_hommel <- p.adjust(p, method = 'hommel')
> length(which(p_hommel < 0.05))
[1] 0
> p_bonferroni <- p.adjust(p, method = 'bonferroni')
> length(which(p_bonferroni < 0.05))
[1] 0
> p_BH <- p.adjust(p, method = 'BH')
> length(which(p_BH < 0.05))
[1] 42
> p_BY <- p.adjust(p, method = 'BY')
> length(which(p_BY < 0.05))
[1] 0

BH法のみ42組が残り、他の補正法では全て有意な検定は0となりました。では、他のシミュレーションデータではどうでしょうか。例えば元の一様分布のサンプリング区間のずれを1.0に増やすとこうなります。

> set.seed(201)
> x1 <- matrix(runif(200000, 10, 15), ncol = 20, nrow = 10000, byrow = T)
> set.seed(202)
> x2 <- matrix(runif(200000, 11, 16), ncol = 20, nrow = 10000, byrow = T)
> p <- c()
> for (i in 1:10000){
+ tmp <- t.test(x1[i, ], x2[i, ])
+ p <- c(p, tmp$p.value)
+ }
> hist(p)
> length(which(p < 0.05))
[1] 5623
> p_holm <- p.adjust(p, method = 'holm')
> length(which(p_holm < 0.05))
[1] 41
> p_hochberg <- p.adjust(p, method = 'hochberg')
> length(which(p_hochberg < 0.05))
[1] 41
> p_hommel <- p.adjust(p, method = 'hommel')
> length(which(p_hommel < 0.05))
[1] 43
> p_bonferroni <- p.adjust(p, method = 'bonferroni')
> length(which(p_bonferroni < 0.05))
[1] 41
> p_BH <- p.adjust(p, method = 'BH')
> length(which(p_BH < 0.05))
[1] 4080
> p_BY <- p.adjust(p, method = 'BY')
> length(which(p_BY < 0.05))
[1] 300

f:id:TJO:20180612091045p:plain

先ほどとは異なり、補正法ごとの差が出てきています。やはりBH法がかなり甘く出ている一方で、FDR系の手法であるBY法がFWER系の手法群よりも甘く出ているという結果になっています。


そこで、今度はデモ用データセットを一様分布ではなく正規分布に変え、ついでに乱数シードも変えてみます。

> set.seed(301)
> x1 <- matrix(rnorm(200000, 12, 3), ncol = 20, nrow = 10000, byrow = T)
> set.seed(302)
> x2 <- matrix(rnorm(200000, 11, 6), ncol = 20, nrow = 10000, byrow = T)
> p <- c()
> for (i in 1:10000){
+ tmp <- t.test(x1[i, ], x2[i, ])
+ p <- c(p, tmp$p.value)
+ }
> hist(p)
> length(which(p < 0.05))
[1] 1027
> p_holm <- p.adjust(p, method = 'holm')
> length(which(p_holm < 0.05))
[1] 1
> p_hochberg <- p.adjust(p, method = 'hochberg')
> length(which(p_hochberg < 0.05))
[1] 1
> p_hommel <- p.adjust(p, method = 'hommel')
> length(which(p_hommel < 0.05))
[1] 1
> p_bonferroni <- p.adjust(p, method = 'bonferroni')
> length(which(p_bonferroni < 0.05))
[1] 1
> p_BH <- p.adjust(p, method = 'BH')
> length(which(p_BH < 0.05))
[1] 3
> p_BY <- p.adjust(p, method = 'BY')
> length(which(p_BY < 0.05))
[1] 0

f:id:TJO:20180612091602p:plain

一つ前の例とは異なり、BH法が相変わらず多め(と言っても3組だけですが)に出る一方で、FDR系のBY法がむしろFWER系の手法群よりも厳しくなり0組という結果を返しています。


最後に、同じように正規分布からデモ用データセットを生成した例ですが、これはちょっと様相が変わって見えます。

> set.seed(501)
> x1 <- matrix(rnorm(200000, 20, 8), ncol = 20, nrow = 10000, byrow = T)
> set.seed(502)
> x2 <- matrix(rnorm(200000, 18, 4), ncol = 20, nrow = 10000, byrow = T)
> p <- c()
> for (i in 1:10000){
+ tmp <- t.test(x1[i, ], x2[i, ])
+ p <- c(p, tmp$p.value)
+ }
> hist(p)
> length(which(p < 0.05))
[1] 1677
> p_holm <- p.adjust(p, method = 'holm')
> length(which(p_holm < 0.05))
[1] 3
> p_hochberg <- p.adjust(p, method = 'hochberg')
> length(which(p_hochberg < 0.05))
[1] 3
> p_hommel <- p.adjust(p, method = 'hommel')
> length(which(p_hommel < 0.05))
[1] 3
> p_bonferroni <- p.adjust(p, method = 'bonferroni')
> length(which(p_bonferroni < 0.05))
[1] 3
> p_BH <- p.adjust(p, method = 'BH')
> length(which(p_BH < 0.05))
[1] 4
> p_BY <- p.adjust(p, method = 'BY')
> length(which(p_BY < 0.05))
[1] 0

f:id:TJO:20180612092311p:plain

BY法以外ではほとんど結果が変わりません。しかも、比較的甘いはずのBY法のみが0組という結果を返しています。これらの実験結果から、p値の分布を変えることで多重比較補正法はかなり振る舞いが変わるということが見て取れます。


余談ながら。わざと最初に得られたp値のうち例えば任意の区間(0.2以下とか0.8以上とか)をupsamplingして水増しするなどしてp値の分布を大きく歪ませてみたらどうなるのかなと思って試してみたら、それぞれの多重比較補正手法の定義*6を考えれば当たり前の話ですが「p < 0.05の個数が変わらない限りはそれほど大きくは結果は変わらない」ということが分かりました。大した手間はかからないので、興味のある方はお手元でお試しください。


最後に


これは例えばAndrew Gelmanあたりも言っている話のようですが、そもそも多重比較という考え方が頻度主義的な点推定に拠った考え方であり、点推定である以上は何かしらの偏りを持ってしまう可能性が否定できない。故に、ベイジアン的枠組みによって「分布」の情報を推定することでパラメータの信頼性に関わる情報をもっと増やすべきだ、という話もあるようです(以下の記事ではeffect sizeの推定を例として挙げています)。

普段からRStanとか好んで使っている身としてはそりゃその通りだと思うのですが、世の中一般に行われている統計分析の一つとして多重比較補正について知っておいても損はないかなと、思う次第です。。。

*1:興味のある人は論文に加えてここを読むとよろしいかと Statistical parametric mapping (SPM) - Scholarpedia

*2:真ではないという意味

*3:実際には10万画素ぐらいあるがヒト脳実質をカバーする範囲だけに解剖学的データから絞り込むということをしていた

*4:7年前に書いた論文でブートストラップ法を使った多重比較補正をやったことがある

*5:起動時に自動的に読み込まれる

*6:特にステップワイズ法は「順位」と「生p値」との関係だけが重要と言っているようなものなので