手当り次第に効果を探すあなたのための多重比較補正
調査・研究は泥縄式になりがちですが、安易に統計的仮説検定で偽陽性が入り込みやすくなります。多重比較補正で、まやかしの統計的有意性を抑制しましょう。
手当り次第に効果を探す研究はぼちぼちあります。例えば薬害には吐き気や蕁麻疹など多岐に渡る症状がありますが、どれか一つだけ現れても薬害です。
有意水準が\(\alpha\)の統計的仮説検定を\(n\)個の正しい帰無仮説に対して行ったとき、それぞれのp値が独立として、ひとつ以上の仮説が統計的に有意になる確率ファミリーワイズエラー率（FWER）は
\[\begin{align} 1-(1-\alpha)^n &= 1 - (1 - {}_n\!C_1 \alpha + {}_n\!C_2 \alpha^2 - \cdots + (-1)^n {}_n\!C_n \alpha^n )\\ &\approx {}_n\!C_1 \alpha = n\alpha \end{align}\]
となります。慣習的に\(\alpha\)は\(0.1\)以下の小さい値になり、二乗以上の項は微々たるものなので無視しました。有意水準\(\alpha\)で検定しているのに、実際はその\(n\)倍ぐらい緩く評価することになります。これを補正するのが多重比較補正です。
もっとも基本的な多重比較補正のBonferroni補正は\(\alpha\)を\(n\)で割って、個々の検定の棄却域にします。単純ですが理に適っているように思えますね。個々の検定の有意水準を\(\alpha/n\)にすれば、どれかが棄却されれば有意という検定全体の有意水準はおおよそ\(\alpha\)にできそうです。しかし、そうは簡単に行きません。検出力が低く、偽陰性が出やすいことが知られています。
Bonferroni補正は、ある仮説が棄却されるとき他の仮説は棄却されない無い場合のみ全体の有意水準が\(\alpha\)になり、そうでない場合はもっと低くなります1。極端な場合で、n個のp値がそれぞれ\(\alpha/n\)より大きく、\(\alpha\)より小さいとしましょう。発熱と倦怠感のように相関がある現象をそれぞれ分析していれば、そういう事もありえます。補正なしでもすべてが有意であれば、全体でも有意と見做すのが直観的です。しかし、多重比較補正をかけると帰無仮説をどれも棄却できないことになります。Bonferroni補正はまた、一番低いp値の分布にあわせた補正です。二番目、三番目に低いp値は異なる分布に従うので、Bonferroni補正では偽陰性の確率が高くなります。
検出力の低さ（偽陰性の多さ）をどう減らすかが問題になり、多重比較補正は多様な手法が提案されています。また、p値と比較する有意水準を補正する他に、検定統計量の最大値の分布から全体のp値を計算する方法もあります。
1 p値を補正する手法
Rのp.adjustでサポートしている手法を俯瞰します。大別すると、p値をFWERを保証する値に変換する手法と、偽陽性率（FDR2）を保証する値に変換する手法に分かれます。検証型の研究ではFWERが問題になり、発見的研究の場合はFDRが問題にされることが多いようです。
1.1 Bonferroni補正
既に説明した有意水準\(\alpha\)を仮説の数\(n\)で割る補正です。ブールの不等式（確率測度の劣加法性）から、p値が独立であればFWERは\(\alpha\)になり、相関があれば\(\alpha\)未満になることが言えます。\(\mbox{FWER} \le \alpha\)が数理的に保証されています。
Bonferroni補正をかけたあとのFWERは以下のようになります。
\[\begin{align} \mbox{FWER} &= P \bigg ( \bigcup^{n_0}_{j=1} \bigg \{ p_j \leq \frac{\alpha}{n} \bigg\} \bigg ) \end{align}\]
ここで\(n_0\)は、正しい帰無仮説（対立仮説が間違っている場合）の数で、\(n_0 \leq n\)です。\(j\)はどの帰無仮説かをあらわす添字、\(p_j\)は仮説\(j\)のp値です。\(1\)から\(n_0\)までの正しい帰無仮説のp値が\(\alpha/n\)以下となる事象の和集合となる事象が起きる確率となります。
\[\begin{align} P \bigg ( \bigcup^{n_0}_{j=1} \bigg \{ p_j \leq \frac{\alpha}{n} \bigg\} \bigg ) &\leq \sum^{n_0}_{j=1} P\bigg( \bigg \{ p_j \leq \frac{\alpha}{n} \bigg \} \bigg) \end{align}\]
確率測度の劣加法性から、和集合の確率は、和集合の要素の事象の確率の和よりも小さくなります。複数の事象が同時に発生する確率があれば、どれか一つ以上の事象が生じる本当の確率は、同時発生を考慮しない個々の事象の発生確率を足した値よりも小さくなるということです。
\[\begin{align} \sum^{n_0}_{j=1} P\bigg( \bigg \{ p_j \leq \frac{\alpha}{n} \bigg \} \bigg) &= \sum^{n_0}_{j=1} \frac{\alpha}{n} = n_0 \frac{\alpha}{n} \leq n \frac{\alpha}{n} = \alpha \end{align}\]
事象\(p_j \leq \alpha/n\)が生じる確率は、p値が有意水準\(\alpha/n\)以下になる確率なので、有意水準の定義から\(\alpha/n\)です。あとは\(n_0 \leq n\)に注意すると、\(\mbox{FWER} \leq \alpha\)が導出できます。
RでBonferroni補正を使うのは簡単です。
[1] 0.048 0.060 0.104 0.152p値を\(5\)倍にしてくれていますね。このp値の集合では、5%有意1つ、10%有意1つです。
しかし、現在では積極的に使うべきではありません。
1.2 Holm-Bonferroni補正
Holm補正とも呼びます。昇順で並べたp値の順番ごとに棄却域を変える手法です。p値の集合の要素を値の小さい方から検定していきます。\(j\)番目のp値の棄却域は\(\alpha/j\)とします。だんだんと棄却域が大きくなり、有意になりやすくなります。一度でも帰無仮説が棄却されなければ、終了です。残りのp値は検定しません。棄却され、かつ未検定のp値が残っていれば、次に小さいp値の検定を行います。
Bonferroni補正よりも一様に（2番目以降の）検出力が高いです。また、二番目以降の棄却域は広くなりますが、どのp値の検定においても\(\mbox{FWER} \le \alpha\)が数理的に保証されます3。
[1] 0.048 0.048 0.052 0.052このp値の集合では、5%有意2つ、10%有意2つになっています。
1.3 Benjamini-Hochberg補正
Benjamini-Hochberg補正は、BH法と略されることもあります。降順で並べたp値の順番ごとに棄却域を変える手法です。p値の集合の要素を値が大きい方から検定していきます。\(j\)番目のp値の棄却域は\(\alpha/j\)とします。だんだんと棄却域が小さくなり、有意になりづらくなります。一度でも帰無仮説が棄却されれば、終了です。残りのp値に対応する帰無仮説も棄却されたものとします。棄却されず、かつ未検定のp値が残っていれば、次に小さいp値の検定を行います。
BH法は補正前の数字がFWERになることを保証しません。ただし、FDR以下の値になることが数理的に保証されています。\(n\)を帰無仮説の数、\(n_0\)はを真である帰無仮説の数とします。\(R\)を棄却された帰無仮説の数、\(V\)を棄却された帰無仮説で真であった数とします。FDPを\(V/\max(1, R)\)と、FDRを\(E[\mbox{FDP}]\)と定義します。FDPが観測値、FDRはその期待値になります。\(R\)は\(0\)になりうるので、分母が変則的に定義しています。\(R \geq V\)なので、\(R=0\)ならば\(V=0\)となり、棄却された帰無仮説がなければFDPは\(0\)です。FDPは帰無仮説ごとに分解できます。 \[ \mbox{FDP} = \sum_{j=1}^{n_0} \frac{V_j}{\max(1, R)} \] \(j\)は帰無仮説を表す添字で、\(V_j\)は帰無仮説ごとの棄却数です。\(0\)か\(1\)になります。級数の内側の項を書き直します。 \[\begin{align} &= \sum_{j=1}^{n_0} \sum_{k=0}^{n} \frac{V_j}{\max(1, R)}I_{R=k} \end{align}\] \(I_{R=k}\)は\(R\)と\(k\)が等しいときは\(1\)、異なるときは\(0\)になる指示関数です。級数の内側に級数を考えてややこしいわけですが、棄却される帰無仮説の数ごとに\(1\)か\(0\)かを計算し、合計します。合算しても\(1\)か\(0\)にしかなりませんが。 \[\begin{align} &= \sum_{j=1}^{n_0} \sum_{k=1}^{n} \frac{V_j}{\max(1, R)}I_{R=k} \\ \end{align}\] \(k=0\)のときは\(V_j=0\)なので、\(k=1\)から級数の計算をしても同じ値になります。 \[\begin{align} &= \sum_{j=1}^{n_0} \sum_{k=1}^{n} \frac{V_j}{k}I_{R=k} \end{align}\] \(R \ne k\)のときは\(I_{R=k}=0\)なので、\(\max(1, R)\)は\(k\)に置き換えても同じです。 \[\begin{align} &= \sum_{j=1}^{n_0} \sum_{k=1}^{n} \frac{I_{p_j < \frac{k}{n}\alpha}}{k}I_{R=k} \end{align}\] \(V_j\)は\(1\)か\(0\)かの棄却された帰無仮説の中で真である数ですから、棄却される条件が満たされているか否かの指示関数と等しくなります。\(p_j\)は帰無仮説\(j\)のp値です。BH法では、ひとつの帰無仮説が棄却されると、残り全部の帰無仮説も棄却して終了となります。\(k\)は、帰無仮説\(j\)を含めた残っている帰無仮説の数になります。BH法の有意水準は、検定せずに残っている帰無仮説が多ければ多いほど有意水準が高くなるため、\(p_j < k\alpha/n\)と言う条件のついた指示関数になります。 FDRを計算しましょう。 \[\begin{align} \mbox{FDR} &= E[\mbox{FDP}] = E\bigg[ \sum_{j=1}^{n} \sum_{k=1}^{n} \frac{I_{p_j < \frac{k}{n}\alpha}}{k}I_{R=k} \bigg] \\ &=\sum_{j=1}^{n_0} \sum_{k=1}^{n} \frac{ E[ I_{p_j < \frac{k}{n}\alpha} ] }{k}I_{R=k} \end{align}\] 確率変数は\(p_j\)だけです。\(R\)が確率変数ではないことに注意しましょう。 \[\begin{align} &=\sum_{j=1}^{n_0} \sum_{k=1}^{n} \frac{ P({p_j < \frac{k}{n}\alpha}) }{k}I_{R=k} \\ &=\sum_{j=1}^{n_0} \sum_{k=1}^{n} \frac{ \frac{k}{n}\alpha }{k}I_{R=k} \end{align}\] p値が有意水準未満の確率は、有意水準です。 \[\begin{align} &=\sum_{j=1}^{n_0} \sum_{k=1}^{n} \frac{\alpha}{n} I_{R=k} \\ &=\sum_{j=1}^{n_0} \frac{\alpha}{n} = \frac{n_0}{n} \alpha \leq \alpha \end{align}\] \(\sum_{k=1}^{n}\)の中で\(I_{R=k}\)が\(1\)の項は一つ、かつ\(I_{R=k}\)以外に\(k\)が使われていないので、\(\sum_{k=1}^{n}\)と\(I_{R=k}\)は消せます。
これで、\(\mbox{FDR} \leq \alpha\)が示せました。なお、偽陽性率が100%だとFDRはLWERに一致し、そうでない場合はFDRはLWERよりも小さくなります
多数の症状を調べて薬害があるか否か断定するような、研究仮説の検証には向かない手法になります。一方、検出力が高く、すべてのp値が\(\alpha/2\)より大きく\(\alpha\)よりも小さいような場合でも、有意になります。薬害の可能性がないか探していく発見的な研究においては、有用な方法です。
[1] 0.03000000 0.03000000 0.03466667 0.03800000すべて5%有意になりました。
1.4 Benjamini-Yekutieli
Benjamini-YekutieliはBY法とも記述されます。BH補正を改良したものと説明されています（松田 (2008)）が、保守的になるようです。
[1] 0.06250000 0.06250000 0.07222222 0.07916667すべて10%有意になりました。
1.5 Hochberg補正とHommel補正
p.adjustではHochberg (1988)とHommel
(1988)で提案された補正方法も、methodにhochbergもしくはhommelを指定することで使う事ができます。これらは「検定統計量間に独立ないし正の相関が仮定できる場合にしか
FWER を制御することができない」（杉谷・森川
(2017)）そうですが、条件を満たしているときは検出力が高いようです。
1.6
p.adjustではサポートされない手法
多重比較補正の手法は多様で、p.adjustでサポートされていないŠidák補正なども、分野によってはよく使われるようです。サポートしているパッケージが見つからない場合もあるかも知れませんが、計算手順はそう複雑ではないことが多いです。
またゲートキーピング戦略といって、重要性に応じて帰無仮説に順番をつけ、多重比較補正なしに検定を行っていき、棄却できない帰無仮説が出たらそこで終了して、残りの帰無仮説は棄却できないと見做す方法もあります。この戦略では、FWERが個々の有意水準以下に保証され、かつ検出力も高くなりますが、帰無仮説の重要性に順番をつける学術的な理屈が必要になります。
2 検定統計量から補正p値を計算する手法
p値からFWERもしくはFDRを保証する値を計算する手法を見てきましたが、FWERとなる検定統計量の最大値の分布のp値を計算する手法もあります。サブサンプルとサブサンプルの違いを探す多群比較でよく使われています。
2.1 Tukey Honestly Significant Difference
Tukeyの範囲検定、Tukey法とも呼ばれます。サブサンプルとサブサンプルの全組み合わせ（many-to-many）の二群比較のt統計量の最大値の分布から、p値を計算する方法です。サブサンプルごとに要素数が異なる場合は、自由度で補正したTukey-Kramer法を用います。
慣習的にANOVAと併用されることが多いですが、RではANOVAの後処理のようなコマンド体系になっていますが、多群比較が目的であれば前段でANOVAを参照する必要はありません。
Speciesはirisデータフレーム中の群をあらわすfactor型の列名です。
ヘルプに明記は無いですが、Tukey-Kramer法になっているそうで、サンプルサイズの不均一性は気にしなくても大丈夫です。
2.2 ガブリエル比較区間
信頼区間は多群比較に使えないので、プロットで群間比較を一望したい場合はt検定に対応したガブリエル比較区間が便利です。また、ガブリエル比較区間はTurkey HSDにあわせて区間調整をしてくれます。
2.3 Dunnett補正
サブサンプルの数をnとして、Tukey HSDはn対n比較になります。しかし、薬剤の効果の検証などでは統制群と介入群A、統制群と介入群B、統制群と介入群C…と言うようにmany-to-one比較になることが通例で、Tukey HSDは過度の補正になります。many-to-oneに対応するのがDunnett補正です。
Rではmultcompパッケージでサポートされています。irisデータセットで試してみましょう。
library(multcomp)
summary(glht(aov(Sepal.Width ~ Species, iris),
linfct = mcp(Species = "Dunnett")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Dunnett Contrasts
Fit: aov(formula = Sepal.Width ~ Species, data = iris)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
versicolor - setosa == 0 -0.65800 0.06794 -9.685 < 1e-10 ***
virginica - setosa == 0 -0.45400 0.06794 -6.683 9.07e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)Speciesはirisデータフレーム中の群をあらわすfactor型の列名です。
2.3.1 コントロールの指定
統制群の指定をしたい事があると思いますが、なぜかmultcompではその指定をサポートしておらず、最初の要素が使われます。そこで、群をあらわすfactor型の変数をつくりなおします。
library(multcomp)
# 元の1:setosa, 2:versicolor, 3:virginicaという並びを、
# 1:versicolor, 2:virginica, 3:setosaに入れ替える
sorted_gn <- factor(
as.character(iris$Species),
levels = levels(iris$Species)[c(3, 1, 2)])
summary(r_glht_sorted <- glht(aov(Sepal.Width ~ sorted_gn, iris),
linfct = mcp(sorted_gn = "Dunnett")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Dunnett Contrasts
Fit: aov(formula = Sepal.Width ~ sorted_gn, data = iris)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
setosa - virginica == 0 0.45400 0.06794 6.683 9.07e-10 ***
versicolor - virginica == 0 -0.20400 0.06794 -3.003 0.00607 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)3 利用上の注意
多重比較補正は、複数の統計的仮説検定のどれかが有意になれば主張（研究仮説）が示されるようなとき、発見的に手当り次第に効果を探したときには偽要請を抑制する有用な手法です。ただし、盲目的に利用すると弊害が出ます。
頑強性のテストなどで、説明変数が異なる計量モデルを試した場合を考えます。建前として検証に用いる計量モデルは一つです。統制変数の有無があっても、異なるモデルの同じ説明変数の係数のp値には往々にして強い相関があるので、Bonferroni型の補正は強すぎます。実際、複数のモデルで推定して、すべてのモデルで統計的有意になった、主要なモデルで有意になり頑強性テストでも大きくは変化しなかったと言うような記述で済ます事が多く、暗に緩くFDRで検証するようなことが多いです。