Rによる多重比較法
はじめに
本記事では多重比較法のRでの実装を紹介する。具体的には同時分布を利用した手法（Dunnet, Tukey, Steel, Steel-Dwass法）、p値を利用した手法（Bonferroni, Holm, Hochberg, Hommel法）、グラフィカルアプローチを扱う。ただし各手法のアルゴリズムや証明については扱わない。アルゴリズムや証明については、適宜Rのドキュメントもしくは計量生物学36巻の特集等を参照することを推奨する。
【参考文献】
はじめに_多重比較法
多重比較の問題は複数のp値を算出した時に生じることが知られており、医学研究においても多重比較の問題は様々な状況で生じる。例えば、2群のランダム化比較試験での主要評価項目と副次評価項目の検定、サブグループ解析、中間解析あるいは探索的な試験で複数の因子のアウトカムへの影響の検討などが挙げられる。
本稿では基本的な多重比較法のRのプログラムの紹介を行う。しかし改めて調べてみると、多重比較法について体系的に書かれた文献はネット上では意外と少ない（ように感じた）。また多重比較法は全体図が掴みにくいため、本記事内で扱う手法を明らかにする意も込めて筆者なりに簡単にではあるが整理を試みた。ただし表現は必ずしも厳密ではないので注意して頂きたい。
まず多重比較法は大きく「検定統計量の同時分布を利用する手法」と「p値のみを利用する手法（検定統計量の周辺分布の利用）」に分けられ、それぞれの手法について概要を示す。
※（読み飛ばしてもらっても構わないが、）本記事では、少なくとも1つの帰無仮説が間違っている場合にGlobal帰無仮説を棄却する状況、つまりat-least one手順でのFWERの強制御を想定している。可能な限り文脈で識別可能なように努めたが、FWERと仮説群を構成する各帰無仮説に対する有意水準を明示的に区別せず記載している。
検定統計量の同時分布を利用する手法
例えば、Placebo群、低用量群、高用量群の3群のランダム化比較試験における比較を考える。Placebo群と低用量群の比較と、Placebo群と高用量群の比較をそれぞれt検定で行うと2回の比較がなされる。その際にしばしば、帰無仮説が正しい下でどちらかの仮説が誤って棄却される確率は、
検定統計量の同時分布を利用する手法はさらにパラメトリックな手法とノンパラメトリックな手法、リサンプリング等による同時分布の推定の手法に分けられる。
パラメトリックな手法
パラメトリックな手法というのはアウトカムの分布に仮定を置く手法であり、アウトカムに正規分布の仮定を置き先ほど示したような対照群とのt検定を複数回繰り返す場合の手法としてはDunnet法が知られている。他にも同様に正規分布の仮定の下で対比較をする場合の手法はTukey法がある。
ノンパラメトリックな手法
パラメトリックな手法に対して、アウトカムに分布の仮定を置かず順序に対する同時分布により多重比較を調整する方法をノンパラメトリックな手法と呼ぶ。Dunnet法に対応するノンパラメトリックな手法としてSteel法、Tukey法に対応する手法としてSteel-Dwass法が知られている。
同時分布の推定
また検定統計量の同時分布を何らかの方法で推定する手法も考えられる。リサンプリングによるシミュレーションや並べ替えを考えるPermutation法等がある。
p値のみを利用する手法
検定統計量の同時分布を利用する手法に対して算出されたp値だけをもとに調整する方法も存在する。検定統計量の同時分布を利用する手法に比べ相関を考慮しないため検出力の低下は避けられないが、同時分布の誤特定によらずαを制御することができ、また同時分布を考えにくい状況でも適用できるためよく利用される。例えば、Bonferroni法はp値のみを利用する手法の一つである。
p値のみを利用する手法でも特定の分布の仮定はしないが、p値間に「非負の相関がある」という仮定を置くような手法は存在する。
本記事内では検定統計量の同時分布を利用する手法としてDunnet, Tukey, Steel, Steel-Dwass法を、p値のみを利用する手法としてBonferroni, Holm, Hochberg, Hommel法を紹介する。これらの手法は帰無仮説間の構造を考慮しない無構造仮説群に対する多重比較である。その後、構造仮説群に対する手法としてグラフィカルアプローチを紹介する。
筆者自身、多重比較法の勉強量が乏しく現在提案されている手法のうち紹介しているのはごく一部である。しかし無構造仮説群に対する基本的な多重比較を理解した上で、構造化仮説群にはグラフィカルアプローチを適応することである程度は広範な状況に対応できるのではないかと考えている。
本記事を読む際の注意事項だが、多重比較法においてはp値を固定してαを調整する方法と逆にαを固定してp値を調整する方法があり、これらは同等の検定結果を示す。例えば、Bonferroni法で
目次
- Package
- 検定統計量の同時分布を利用する手法
- データ生成
- パラメトリックな手法
- Dunnet法
- Tukey法
- ノンパラメトリックな手法
- Steel法
- Steel-Dwass法
- p値のみを利用する手法
- Bonferroni, Holm, Hochberg, Hommel法
- 重み付き Bonferroni, Holm, Hochberg, Hommel法
- 構造化仮説群（グラフィカルアプローチ）
- 固定順序法/ Fallback法/ SerialなChain法
- CycleなChain法
- Holm法の再現
- まとめ
- （メモ）今後、追加する可能性ある解析
Package
利用するPackageは以下の通り。gMCPはJavaをインストールしておく必要がある。
- tidyverse
- multcomp
- nparcomp
- Mediana
- gMCP
#package
library(tidyverse)
library(multcomp)
library(nparcomp)
library(Mediana)
library(gMCP)
検定統計量の同時分布を利用する手法
データ生成
仮想的に、4群で各群50名、アウトカムが正規分布に従うデータを発生させた。
#データ生成
set.seed(1234)
n <- 200
th1 <-
data.frame(ID = c(1:n)) %>%
mutate(Group = case_when(
ID <= n/4 ~ "A",
ID <= 2*n/4 ~ "B",
ID <= 3*n/4 ~ "C",
ID <= n ~ "D"))
th1$Group <- as.factor(th1$Group)
y <- c(
round(rnorm(n=50, mean=0, sd=1),2),
round(rnorm(n=50, mean=0.2, sd=1),2),
round(rnorm(n=50, mean=0.4, sd=1),2),
round(rnorm(n=50, mean=-0.2, sd=1),2)
)
th2 <- cbind(th1, y)
パラメトリックな手法
パラメトリックな手法として、multcompを利用してDunnet法とTukey法を実装した。
Dunnet法
アウトカムに正規分布の仮定を置き、特定の対照群とその他の群の比較を行う場合の多重調整法であるDunnet法は以下のように行える。
#Dunnet
amod <- aov(y ~ Group, data = th2)
Dunnet <-(glht(amod, linfct = mcp(Group = "Dunnet")))
summary(Dunnet)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: aov(formula = y ~ Group, data = th2)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 0.7942 0.2003 3.966 <0.001 ***
## C - A == 0 0.8736 0.2003 4.362 <0.001 ***
## D - A == 0 0.3158 0.2003 1.577 0.27
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Tukey法
アウトカムに正規分布の仮定を置き、すべての群間での対比較を行う場合の多重調整法であるTukey法は以下のように行える。
#Tukey
Tukey <- (glht(amod, linfct = mcp(Group = "Tukey")))
summary(Tukey)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = y ~ Group, data = th2)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 0.7942 0.2003 3.966 <0.001 ***
## C - A == 0 0.8736 0.2003 4.362 <0.001 ***
## D - A == 0 0.3158 0.2003 1.577 0.3942
## C - B == 0 0.0794 0.2003 0.396 0.9788
## D - B == 0 -0.4784 0.2003 -2.389 0.0828 .
## D - C == 0 -0.5578 0.2003 -2.785 0.0296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
ノンパラメトリックな手法
一応、nparcompを利用したSteel法とSteel-Dwass法を示すが、これらの手法は一義な解析ではなく、分布の仮定や分散の仮定などの違いが存在する。そのため実際に多重比較法としてこれらの手法を利用する場合には、きちんと理解したうえで適切なプログラムを用いることを推奨する。
Steel法
Dunnet法に対応する手法であり、正規分布の仮定のかわりに順序統計量を利用するSteel法は以下のように実装できる。
#Steel
Steel <-nparcomp(y ~ Group, data=th2, asy.method = "mult.t",
type = "Dunnett", control = "A",alternative = "two.sided",info = FALSE)
summary(Steel)
##
## #------------Nonparametric Multiple Comparisons for relative contrast effects----------#
##
## - Alternative Hypothesis: True relative contrast effect p is not equal to 1/2
## - Estimation Method: Global Pseudo ranks
## - Type of Contrast : Dunnett
## - Confidence Level: 95 %
## - Method = Multi - T with 92 DF
##
## - Estimation Method: Pairwise rankings
##
## #---------------------------Interpretation--------------------------------------------#
## p(a,b) > 1/2 : b tends to be larger than a
## #-------------------------------------------------------------------------------------#
##
## #----Data Info-------------------------------------------------------------------------#
## Sample Size
## 1 A 50
## 2 B 50
## 3 C 50
## 4 D 50
##
## #----Contrast--------------------------------------------------------------------------#
## A B C D
## B - A -1 1 0 0
## C - A -1 0 1 0
## D - A -1 0 0 1
##
## #----Analysis--------------------------------------------------------------------------#
## Comparison Estimator Lower Upper Statistic p.Value
## 1 p( A , B ) 0.727 0.604 0.850 4.423957 7.952613e-05
## 2 p( A , C ) 0.780 0.664 0.895 5.808727 1.865093e-07
## 3 p( A , D ) 0.606 0.469 0.744 1.849751 1.687196e-01
##
## #----Overall---------------------------------------------------------------------------#
## Quantile p.Value
## 1 2.395503 1.865093e-07
##
## #--------------------------------------------------------------------------------------#
Steel-Dwass法
Tukey法に対応する手法であり、正規分布の仮定のかわりに順序統計量を利用するSteel-Dwass法は以下のように実装できる。
#Steel_Dwass
Steel_Dwass <-nparcomp(y ~ Group, data=th2, asy.method = "mult.t",
type = "Tukey",alternative = "two.sided",info = FALSE)
summary(Steel_Dwass)
##
## #------------Nonparametric Multiple Comparisons for relative contrast effects----------#
##
## - Alternative Hypothesis: True relative contrast effect p is not equal to 1/2
## - Estimation Method: Global Pseudo ranks
## - Type of Contrast : Tukey
## - Confidence Level: 95 %
## - Method = Multi - T with 88 DF
##
## - Estimation Method: Pairwise rankings
##
## #---------------------------Interpretation--------------------------------------------#
## p(a,b) > 1/2 : b tends to be larger than a
## #-------------------------------------------------------------------------------------#
##
## #----Data Info-------------------------------------------------------------------------#
## Sample Size
## 1 A 50
## 2 B 50
## 3 C 50
## 4 D 50
##
## #----Contrast--------------------------------------------------------------------------#
## A B C D
## B - A -1 1 0 0
## C - A -1 0 1 0
## D - A -1 0 0 1
## C - B 0 -1 1 0
## D - B 0 -1 0 1
## D - C 0 0 -1 1
##
## #----Analysis--------------------------------------------------------------------------#
## Comparison Estimator Lower Upper Statistic p.Value
## 1 p( A , B ) 0.727 0.592 0.862 4.4239575 1.289571e-04
## 2 p( A , C ) 0.780 0.653 0.906 5.8087268 2.837200e-07
## 3 p( A , D ) 0.606 0.455 0.758 1.8497508 2.676542e-01
## 4 p( B , C ) 0.544 0.389 0.698 0.7413576 9.075590e-01
## 5 p( B , D ) 0.381 0.233 0.530 -2.1063419 1.638790e-01
## 6 p( C , D ) 0.343 0.196 0.491 -2.7947842 3.253740e-02
##
## #----Overall---------------------------------------------------------------------------#
## Quantile p.Value
## 1 2.630279 2.8372e-07
##
## #--------------------------------------------------------------------------------------#
p値のみを利用する手法
Bonferroni, Holm, Hochberg, Hommel法
ここでは、Bonferroni法、Holm法、Hochberg法、Hommel法の4つの多重比較法を実装した。一般に、検出力はBonferroni法から順にHolm法、Hochberg法、Hommel法と大きくなる。BonferroniとHolmはp値に対してなんら仮定を置かないが、Hochberg法、Hommel法はp値に非負の相関が見られるような仮定を置いた上でαを制御する。
他にも基本的な方法としてSidak法やSimes法があるが、Sidak法はBonferroni法とほとんど検出力に差がないため紹介しない。またSimes法はHochberg, Hommel法を理解する上では不可欠だが、Simes法ではGlobal帰無仮説に対する検定でありどの帰無仮説が誤っているかまでは検定できないため紹介しない。
Bonferroni, Holm, Hochberg, Hommel法は得られたp値から調整を行う手法であり、今回は4つの仮説
rawpが調整する前のp値で、各種法での調整p値を算出した。算出された調整p値を有意水準5%と比較して棄却するか判断すればよい。
rawp <- c(0.0082,0.0174, 0.0042,0.0180)
names(rawp) <- c("H1", "H2", "H3", "H4")
proc <- c("BonferroniAdj", "HolmAdj", "HochbergAdj", "HommelAdj")
adjp <- sapply(proc,function(x){
AdjustPvalues(rawp, proc = x)
})
row.names(adjp) <- c("H1", "H2", "H3", "H4")
adjp <- t(adjp)
round(adjp, 4)
## H1 H2 H3 H4
## BonferroniAdj 0.0328 0.0696 0.0168 0.0720
## HolmAdj 0.0246 0.0348 0.0168 0.0348
## HochbergAdj 0.0180 0.0180 0.0168 0.0180
## HommelAdj 0.0180 0.0180 0.0164 0.0180
重み付き Bonferroni, Holm, Hochberg, Hommel法
上記の例では、
全体を通じて第一種の過誤をαに抑える場合に、どの仮説にαをどの程度配分するかを決めることができる。上記の例では、すべての仮説に均等にαを配分しているため4つの仮説に対して1/4ずつ割り当てられている。各仮説に割り当てた
rawp <- c(0.0082,0.0174, 0.0042,0.0180)
names(rawp) <- c("H1", "H2", "H3", "H4")
proc <- c("BonferroniAdj", "HolmAdj", "HochbergAdj", "HommelAdj")
weight <- c(0.5 ,0.2, 0.2, 0.1)
non_weight <- c(1/4 ,1/4, 1/4, 1/4)
weighted_adjp <- sapply(proc,function(x){
AdjustPvalues(rawp, proc = x, par = parameters(weight=weight))
})
row.names(weighted_adjp) <- c("H1", "H2", "H3", "H4")
weighted_adjp <- t(weighted_adjp)
round(adjp, 4)
round(weighted_adjp, 4)
## H1 H2 H3 H4
## BonferroniAdj 0.0164 0.0870 0.0210 0.1800
## HolmAdj 0.0164 0.0261 0.0164 0.0261
## HochbergAdj 0.0131 0.0180 0.0131 0.0180
## HommelAdj 0.0131 0.0180 0.0117 0.0180
構造化仮説群（グラフィカルアプローチ）
重み付き Bonferroni, Holm, Hochberg, Hommel法では各仮説に割り当てるαで優先度を表現したが、より帰無仮説同士の構造がはっきりした場合も存在する。例えば、第三相試験のRCTであれば、主要評価項目と副次評価項目では主要評価項目で有意な結果がでた場合に限り副次評価項目を検討したいというようなことも考えられる。
構造化仮説群の多重調整の方法も複数提案されているが、柔軟でわかりやすいグラフィカルアプローチを紹介する[1]。グラフィカルアプローチではn個の帰無仮説に対する各有意水準とn*nのtransitionmatrixを事前に規定する必要がある。transitionmatrixの
固定順序法/ Fallback法/ SerialなChain法
固定順序法とFallback法では一方向で枝分かれのない順序を指定する。
固定順序法では、H1から順にαで検定を行い、棄却された場合のみ次の検定に進む。一方で、Fallback法では事前にαをそれぞれの帰無仮説に割り付けておく。H1から順に検定を行い、棄却された場合のみ、その検定のαを次の検定に持ち越すことができる。transition1が固定順序法とFallback法でのtransitionmatrixである。
SerialなChain法では、枝分かれするような順序を指定する。ただし一方向のみでCycleは発生しない。
Fallback法と同様に事前にαをそれぞれの帰無仮説に割り付けておく。w13が0の場合は、単にFallback法と同様になる。上のグラフを表現したtransitionmatrixが、transition2である。
rawp <- c(0.0082,0.0174, 0.0042)
names(rawp) <- c("H1", "H2", "H3")
weight1 <- c(1,0,0)
weight2 <- c(0.8, 0.1, 0.1)
transition1 = matrix(c(0, 1, 0, 0, 0, 1, 0, 0, 0),3,3,byrow=TRUE)
transition2 = matrix(c(0, 0.5, 0.5, 0, 0, 1, 0, 0, 0),3,3,byrow=TRUE)
transition1
## [,1] [,2] [,3]
## [1,] 0 1 0
## [2,] 0 0 1
## [3,] 0 0 0
transition2
## [,1] [,2] [,3]
## [1,] 0 0.5 0.5
## [2,] 0 0.0 1.0
## [3,] 0 0.0 0.0
adjp_fix <- AdjustPvalues(rawp, proc = "ChainAdj",
par = parameters
(weight = weight1,
transition = transition1))
adjp_fallback <- AdjustPvalues(rawp, proc = "ChainAdj",
par = parameters
(weight = weight2,
transition = transition1))
adjp_Serialchain <- AdjustPvalues(rawp, proc = "ChainAdj",
par = parameters
(weight = weight2,
transition = transition2))
adjp <- rbind(adjp_fix, adjp_fallback, adjp_Serialchain)
colnames(adjp) <- c("H1", "H2", "H3")
round(adjp, 4)
## H1 H2 H3
## adjp_fix 0.0082 0.0174 0.0174
## adjp_fallback 0.0102 0.0193 0.0193
## adjp_Serialchain 0.0102 0.0348 0.0102
CycleなChain法
SerialなChain法とは異なり、両方向の矢線が存在している場合の多重比較法で、かなり柔軟なアプローチになる。様々な状況に対応することができるが、一例としてHolm法の再現を行った。
Holm法の再現
rawp <- c(0.0082,0.0174, 0.0042,0.0180)
transition3 <- rbind(
H1=c(0,1/3,1/3,1/3),
H2=c(1/3,0,1/3,1/3),
H3=c(1/3,1/3,0,1/3),
H4=c(1/3,1/3,1/3,0))
round(transition3, 4)
## [,1] [,2] [,3] [,4]
## H1 0.0000 0.3333 0.3333 0.3333
## H2 0.3333 0.0000 0.3333 0.3333
## H3 0.3333 0.3333 0.0000 0.3333
## H4 0.3333 0.3333 0.3333 0.0000
weight3 <- c(1/4, 1/4, 1/4, 1/4)
graph <- new("graphMCP", m=transition3, weights=weight3)
adjp <- gMCP(graph, pvalues = rawp)@adjPValues
names(adjp) <- c("H1", "H2", "H3", "H4")
round(adjp, 4)
## H1 H2 H3 H4
## 0.0246 0.0348 0.0168 0.0348
まとめ
本記事内では、無構造仮説群に対する多重比較法として検定統計量の同時分布を利用する手法とp値のみを利用する手法を紹介した。構造化仮説群に対する手法としてはグラフィカルアプローチを紹介した。
臨床試験の計画においてもグラフィカルアプローチを適用すればかなり広範な状況に対応ができると思われる。ただし本記事では閉検定手順等は触れておらず、グラフィカルアプローチでαが制御できることの証明は別途文献を参照して頂きたい。
多重比較法は状況により複雑になりやすく、冒頭でも述べたが筆者自身より勉強する必要を強く感じた。もし誤りを見つけられた方はご連絡頂ければ幸いである。
（メモ）今後、追加する可能性ある解析
- リサンプリング、並べ替えによる同時分布の推定
- ゲートキーピング法
- グラフィカルアプローチの例の追加
- 中間解析
サポートして頂けるとモチベーションに繋がりますのでぜひ宜しくお願いします。データ解析や臨床研究でのご相談があれば、お気軽にTwiiterもしくはメールにてご連絡下さい。
作成者：Masahiro Kondo
作成日：2022/3/11
連絡先：m.kondo1042(at)gmail.com
-
Bretz,F.,Maurer,W.,Brannath,W.,andPosch,M.(2009).A Graphicalapproach tosequentiallyrejective multipletestprocedures.Statisticsinmedicine,28, 586–604. ↩︎
Discussion