Fisher's Least Significant Difference (LSD) testというpost hoc analysisを聞いた。
というのも、N群のデータでどれとどれに有意差があるか、ということを検定したいときに、単純に通りのt検定をすると多重検定補正が必要で云々、となるが、生物学や医学の論文でよくやるのが、ANOVAをして、ひとまずN群のなかで何かは有意差がありそう、というぼやっとした検定をして、ANOVAで有意差があればペアワイズのt検定をしましょう、という感じになっている。
という解析をしたら、査読者から
「それってFisher's LSD test をしたほうがいいんじゃね?統計解析担当者と確認してよ」
というコメントをもらったのでこれってどういうこと?という相談を受けた。
この投稿先はIFで10点以上はあるかなりの有名どこなので、査読者もけっこう知識がありそうだが、かといってこういう有名雑誌でも検定法、特に多重検定の補正についてはたぶんよくわかってない。
Rではagricolaeというパッケージでできるので、これを使った。データサンプルはこちらからぱくった。
LSD testの概念としては、
N群をすべてプールしてLSD、つまり最小の有意差を生み出せる差を計算する。
2対比較をしておいて、LSDを超える差があった場合に、有意水準で有意差がある、とする。
ANOVAからのペアワイズ検定より、検出力は高い。
ただし、この場合は多重検定補正は行われていない(を調整しない)。
ということらしい。ただし、LSD.test 関数では多重検定補正法を入れられる。
データで試してみる。
library(agricolae) # データ dat <- data.frame(Contact = c(21, 20, 26, 46, 35, 13, 41, 30, 42, 26), Hit = c(23, 30, 34, 51, 20, 38, 34, 44, 41, 35), Bump = c(35, 35, 52, 29, 54, 32, 30, 42, 50, 21), Collide = c(44, 40, 33, 45, 45, 30, 46, 34, 49, 44), Smash = c(39, 44, 51, 47, 50, 45, 39, 51, 39, 55) ) M <- colMeans(dat) # 平均 question <- rep(colnames(dat), each=nrow(dat)) # ラベル dat <- data.frame(y=unlist(dat), question=question) # データフレーム化 v <- aov(y ~ question, data=dat) # anova # 補正法を試す pm <- c("none", "holm", "hochberg", "bonferroni", "BH", "BY", "fdr") lsd <- mapply(function(x) LSD.test(v, "question", p.adj=x), pm, SIMPLIFY=FALSE) lsdm <- mapply(function(x) x$statistics$LSD, lsd)
none holm hochberg bonferroni BH BY fdr 8.056414 11.808317 8.056414 11.808317 8.056414 9.908163 8.056414
単にLSDをしただけならば、8の差があれば有意となる。Bonferroni などの補正法を行うと、8よりももっと差がないと有意にならない。
abs(outer(M, M, "-"))
Contact Hit Bump Collide Smash Contact 0 5 8 11 16 Hit 5 0 3 6 11 Bump 8 3 0 3 8 Collide 11 6 3 0 5 Smash 16 11 8 5 0
Smash と Contact, Smash と Hit は平均にかなり差があるので、これは見た感じ有意となりそう。何も補正をしないLSDは8なので、8より大きいこれらの対比較は有意、ということにする。
ここで、普通にペアワイズに検定を繰り返す。何も補正をしなければわりと有意になる比較が多いが、
pairwise.t.test(dat$y, dat$question, p.adj="none")
Pairwise comparisons using t tests with pooled SD data: dat$y and dat$question Bump Collide Contact Hit Collide 0.45716 - - - Contact 0.05156 0.00855 - - Hit 0.45716 0.14060 0.21776 - Smash 0.05156 0.21776 0.00023 0.00855 P value adjustment method: none
Bonferroni にするとSmash vs Hit が有意でなくなる。これはLSD=11.8だが、実際の差は11だから絶妙に有意にならない。
pairwise.t.test(dat$y, dat$question, p.adj="bonferroni")
Pairwise comparisons using t tests with pooled SD data: dat$y and dat$question Bump Collide Contact Hit Collide 1.0000 - - - Contact 0.5156 0.0855 - - Hit 1.0000 1.0000 1.0000 - Smash 0.5156 1.0000 0.0023 0.0855 P value adjustment method: bonferroni
検出力が高いらしいのでシミュレーションをして確かめようと思ったけど保留(たぶんしない。
