初めに
ANOVA + 多重比較はBiologyでは頻繁に使用される解析で、論文でも高頻度で見かけます。
一方で適切な使用ができているのかというとよくわからない場合も多いです。例えばよく取りざたされる話では多重比較をANOVAの事後検定として用いるのは良いだのどうだの という話があります。意外と解析者もよくわからずに使うことも多いのではないでしょうか。
自分としては多重比較があまり良い方法とは思えず、他の方法はないだろうかといつも考えております。 総当たりの検定は要するに仮説がない解析であるわけで、基本的にはクリアな仮説を立てて、わかりやすい手法で解析するのが最も重要かなと思っております。
この考えを1元配置の例に当てはめると、 1. 複数群のデータに対して、仮説なしの多重比較 (全群比較) 2. 複数群のデータに対して、何らかの手法で仮説を生成してから2郡比較 という2つを考えられます。この2つをサンプルサイズの観点で比べてみました。
言い訳
筆者の統計学の理解レベルは学部生程度 (かそれ以下)です。さらに個人的には頻度論よりもベイズのほうが良いと思っていますが、現状の生物学の論文はほとんど頻度論 (というか仮説検定) で書かれているので以下は頻度論の話だけにします。
本題
こんなグラフがあるとします。 (一元配置)
多分普通はOne way ANOVAをして4群比較
という形になるかと思いますが、こういう場合にTukey-HSDを使用すると薬剤B vs 薬剤Cというぜんぜん興味のない検定をやるがために検定の回数が増えてしまい、差がありそうなコントロール vs 薬剤Aが有意ではなくなることがあると思います。
こういうデータに対して、ANOVA後の解析としては、みたい箇所だけ2郡比較を行って検定回数でBoferroni補正を行うのはどうなのか というのが頭をよぎりました。 ただこれはANOVAの結果 (というかグラフ) から仮説を立てて、その後に同じデータで検定をしているという意味でHARKingに近いと思います。
実験前時点で仮説が存在しないので、ANOVAをやっているわけであり、元から仮説が存在するならそこだけを検証する実験を組むのが良いですよね。
そういう意味で多重比較は仮説なしで行うことしかできない という弱さがある気がします。実際に海外の教科書ではTukey HSDなどの多重比較検定を “unplanned test” というふうに呼んでいました。多数の群になると金科玉条のように多重比較が出てくるのにモヤッとしていましたが、まさにこの表現が正しい気がします。
そうなると上記のケースで言えば、考えられる代替案として、ANOVA (もしくは図示) はあくまで仮説生成に使用して、その後の仮説検証を行うのは別のデータセットで行う という方が検出力を上げる意味でも良いのではないか と思いました。
具体的に一例目の一元配置の例で考えるのであれば、4群でのデータを取り、 図示とANOVAで薬剤の効果がありそうだというところまでを主張して、薬剤Aに効果がありそうだという仮説を作成します。その後にctrlと薬剤Aに差があることを、別の実験を組んで検証する という流れです。
この方が最終的な結論を出すまでの必要サンプルサイズが小さくなる気がするのですが、どうなのでしょうか。検証してみたいと思います。
検証1 Tukey HSD
# example1 ANOVA + Tukey HSD
# in for loop, sample size n is increase until p < 0.05 tukey HSD between ctrl vs drugA
# every time total n is stored in vector
# repeat 10000 times
# report mean and sd of n
# starting sample size is 5 for each group
N <- numeric(3000)
for (i in 1:3000){
n <- 3
p <- 1
set.seed(i)
while (p > 0.05){
n <- n + 1
ctrl <- rnorm(n, 10, 2)
drugA <- rnorm(n, 12, 2)
drugB <- rnorm(n, 10, 2)
drugC <- rnorm(n, 10, 2)
df <- data.frame(ctrl, drugA, drugB, drugC) %>% gather(key = "Condition", value = "value")
aov_res <- aov(value ~ Condition, data = df)
tukey_res <- TukeyHSD(aov_res)
p <- tukey_res$`Condition`[,4][1]
}
N[i] <- n * 4
}
mean(N)
sd(N)
これで結果は
mean(N)
[1] 31.576
となりました。
だいたい各群8例ずつで有意になるようです。
検証2 ANOVAで有意になったら2群間のt検定を行う
こちらではまず4群でANOVAを行い、有意になったらctrl vs drug Aの2群比較を行ってみます。
# example2 ANOVA + T-test
# In for loop, sample size n22 is increase until p < 0.05 in F test
# If p < 0.05, then do T-test between ctrl vs. drugA and increase sample size n21 until p < 0.05
# repeat 3000 times
# every time total sample size (N21 + n22) is stored in vector
# report mean and sd of n
N21 <- numeric(3000)
N22 <- numeric(3000)
for (i in 1:3000){
n <- 3
set.seed(i)
p <- 1
while (p > 0.05){
n <- n + 1
ctrl <- rnorm(n, 10, 2)
drugA <- rnorm(n, 12, 2)
drugB <- rnorm(n, 10, 2)
drugC <- rnorm(n, 10, 2)
df <- data.frame(ctrl, drugA, drugB, drugC) %>% gather(key = "Condition", value = "value")
aov_res <- aov(value ~ Condition, data = df)
p <- summary(aov_res)[[1]][[5]][1]
}
n22 <- n * 4
N22[i] <- n22
n <- 3
p = 1
while (p > 0.05){
n <- n + 1
ctrl <- rnorm(n, 10, 2)
drugA <- rnorm(n, 12, 2)
t_res <- t.test(ctrl, drugA)
p <- t_res$p.value
}
N21[i] <- n22 + n *2
}
# total sample size
mean(N21)
# sample size for ANOVA
mean(N22)
これで結果を得ると
mean(N21)
[1] 37.30467
となりました。予想に反してこちらではサンプルサイズが各群9以上は必要なようです。
ANOVAが有意になるのに必要なサンプルサイズを見てみると
mean(N22)
[1] 24.71733
と各群6例以上を必要としていることがわかります。
すなわち2つ目の方法ではANOVAの部分でかなりのサンプルサイズを費やしていることがわかります。
結論
結果として最初から全群比較をした方が必要サンプルサイズは小さいということがわかりました。 今回は全群比較がTukeyでしたが、これがほかの手法であるとまた違うと思いますし、アンバランスな (各群のサンプルサイズが異なる) データだとまた変わってくるかもしれません。さらに言えば、仮説生成のためにANOVAが必須だとも思いません。極端な話適当なところでグラフを描いてみて、それで仮説生成をすれば2番目に必要なサンプルサイズを小さくすることができるかもしれません。
ただ個人的には意外と最初から多重比較してもサンプルサイズはそんなに大きくならないのだなーというので少し勉強になりました。
