Tukey-Kramer法
分散分析ではどの標本間に差があるかまではわからない
分散分析を使うことで、3標本の平均値の間の差を、第1種の過誤を増大させることなく検出できた。しかし、どの標本間に平均値の差があるかまでは分散分析ではわからないのであった。3標本の平均値は異なりました、だけでなく、どの標本に差がありました、と多くの場合は言いたいはずだ。そこで多数の標本の間で、どの標本の間に差があったかを解析する手段が、多重比較multiple comparisonと呼ばれる方法になる。t検定を複数回行う方法も多重比較の一種となるが、以前も言った通り、そのままでは過誤が大きくなるので問題がある。検定を複数回行うときに使える工夫は別の機会に紹介するとして、今回は別の多重比較法を紹介する。その一つがTukey-Kramer法（あるいはTukeyのhonestly significant difference (HSD)）である。分散分析後に、Tukey-Kramer法の検定をすることでどの標本の間に違いがあるかを検出するのである（本当は分散分析は不要であるが、Rの仕様上、分散分析の結果を使って解析する）。
そのまま、本題に移ってもよいが、ちょっと足を止めて考えてほしいことがある。今までさんざん、有意水準の補正がないまま複数回検定するのは悪手だ、と言っておきながら、「分散分析後に多重比較します」はおかしいだろ、と思わないだろうか。一回、考え出すと疑念が膨らむ読者もいるだろう。これは複数回検定の制約の問題に反しないのだろうか？
心配ご無用、今まで述べてきた複数回検定の制約には反しないようになっている。まずは上の図の左を見てほしい。例えば、t検定を3回行う場合、互いの検定結果は残りの検定結果に影響がない。つまり、独立な関係の検定となる。この場合、検定の結果を掛け算できて、1回の検定につき有意にならない確率が0.95なら、3回の検定すべてで有意にならない確率を0.95 × 0.95 × 0.95 = 0.857…と計算できる。そして、以前にも言ったように1からこの値を引くことで危険率が14%に上ることになる。一方で、上図の右のように、分散分析とTukey-Kramer法をはじめとする多重比較法はほとんど包含関係にある。つまり、分散分析でp < 0.05となって初めて、多重比較でもp < 0.05が検出される。逆に分散分析で有意差が出なければ、多重比較を行っても有意差は出ないのだ。ただし、これはあくまでも解析データが仮定に厳密に従っている、つまり分散分析なら等分散の正規分布から生成された標本群であるときで、そこから外れたデータの時は、分散分析で有意だが、多重比較では有意な標本ペアが見つからないことがやその逆がある。これは仕方のないことなので、このような結果になったからと言ってがっかりしないでよい。一応、Sheffe法と呼ばれる多重比較法は、分散分析と結果が一致するようになっているようだが、今やあまり見かけない方法に思える。
*************************************************
ここで、実践的な話をしよう（ただし、個人的な見解であることも念頭においてほしい）。もし、始めから多重比較をするつもりで解析をするのであれば、多重比較の結果だけを報告してもよい。また、分散分析で有意だが、多重比較では有意な標本ペアが見つからないときに、論文に書く場合は、「分散分析では有意だったが（p < ○○）、多重比較では有意な標本ペアがなかった（標本A-B: p = ○○; 標本B-C: p = ○○; ……）」と潔く書いてしまおう（大体はp = 0.05を少し上回った程度のペアがあるはずである）。ただし、結果を見てから報告の仕方を変えるのはやめたほうが良い。近年は、「p < 0.05かどうかが重要」「p < 0.05でないなら筆者の主張に妥当性はない」という、いわゆる「p値主義」的な雑誌の応対も減りつつある（ないとは言わないが……）。あくまでもp < 0.05は勝手に決めた基準であるし、有意か有意ではないかと、注目している平均の差に意味があるかは全く別問題だ。有意水準を上回ろうともその平均値の差に意味があると信じているなら（もちろん妥当な議論も併せたうえで）、堂々と解析結果を開示するほうが、研究者として正しいスタンスである。出版バイアスpublication biasと呼ばれる、有意な結果ばかりが出版に至る、というバイアスが問題視されるようにもなっており、その状況を打破するためにも上記のようなスタンスの論文が増えていく必要がある。さらに、正確なp値を記載することで、その後にほかの研究者がメタ解析で活用できる可能性も広がる。そして、正しく統計を理解している査読者に回ったら、上記のように書いておけば事情をある程度汲んでくれる（理解してない査読者に回ることもあるが、そこに合わせる道理はないはずだ）。間違っても、なんとか有意になるように解析方法をあれこれ変えてはいけない。
*************************************************
Tukey-Kramer法の検定統計量q
この流れもおなじみになってきた。Tukey-Krammer法の検定統計量はqと呼ばれる。このqを、スチューデント化範囲分布Studentized range distributionのもとで、検定を行うのである。けったいな名前な分布名であるが、考え方は今までとそれほど変わらない。この検定統計量qであるが以下のように計算する。
ここで3標本以上あることを前提として、2標本AとBを取り出したとする。XAとXBは各標本の平均、nAとnBは各標本のサンプルサイズである。そして、(sR)^2は分散分析で登場した残差の平均平方和である。この式の形は、分母に来る分散が違えども、検定統計量tと全く同じ形だ。Tukey-Kramer法はt検定を3標本以上に拡張したものなのである。これまでの2標本の時と決定的に違うことがあるとすれば、帰無仮説の下でのq分布は必ずしもスチューデント化範囲分布にはならない。ただ無作為に2標本を選んで統計量qを計算すると、t分布と変わらないからである。帰無仮説の下でのt分布よりもスチューデント化範囲分布は分布の裾が重く、それゆえに帰無仮説が棄却されるためには、2標本のt検定の時よりもより効果量が大きくなくてはならない。スチューデント化範囲分布を利用することで、検定で帰無仮説が棄却されにくくなり、多重検定の問題をクリアしているのである。スチューデント化範囲分布となるのは、帰無仮説が正しく、3標本以上のとき、最大の平均を持つ標本と最小の平均を持つ標本で統計量qを計算した時である。これらのことから、検定の手順が多重比較の手続きに変更されたというよりも、検定統計量tとt分布の多重比較版が検定統計量qとスチューデント化範囲分布に改良されたというイメージが近いだろう。
Tukey-Kramer法と分散分析の整合性ー危険率
さて、ここからは、Tukey-Kramer法を実際に行ってみよう。すでに紹介したが、分散分析とTukey-Kramer法の結果はおおむね一致する。今回は、q分布を作るのではなく、分散分析とTukey-Kramer法をおこなって、それぞれがどれくらいの危険率となり、そして一致するかを検討してみよう。では、平均値の等しい３つの標本に対して分散分析を行う状況を考える。すべての標本が平均50、標準偏差20、サンプルサイズ30にそろえてある。
------------------------------------------------------
library(plotn)
library(multcomp)
#install.packages("multcomp", repos="http://cran.ism.ac.jp/")
#後述のTukey-Kramer法、cld()関数を使うのに必要
library(VennDiagram)
#install.packages("VennDiagram", repos="http://cran.ism.ac.jp/")
#後述のVenn図を作るのに必要
nA <- 30
nB <- 30
nC <- 30
meA <- 50
meB <- 50
meC <- 50
s <- 20
pA <- NULL
pM <- NULL
for (i in 1:10000){
mA <- rnorm(nA, mean = meA, sd = s)
mB <- rnorm(nB, mean = meB, sd = s)
mC <- rnorm(nC, mean = meC, sd = s)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
d$g <- as.factor(d$g)
fit1 <- aov(d ~ g, data = d) #分散分析
fit2 <- summary(fit1)
fit3 <- glht(fit1, linfct = mcp(g = "Tukey")) #Tukey-Kramer法。分散分析の結果fit1を利用している。
fit4 <- summary(fit3)
pA <- c(pA, fit2[[1]]$`Pr(>F)`[1]) #分散分析のp値の取得
pM <- rbind(pM, fit4$test$pvalues[1:3]) #Tukey-Kramer法のp値の取得。A-B、B-C、C-Aのペアの3つ分がある
}
pM <- as.data.frame(pM)
colnames(pM) <- names(fit4$test$tstat)
pM$pmin <- apply(pM, MARGIN = 1, FUN = min) #Tukey-Kramer法のA-B、B-C、C-Aのペアのp値のうち、最小のものを選ぶことで、あとで最低でも1つのペアで有意になるものを選ぶことができる
p <- cbind(pA, pM)
head(p) #同じ列のp値は同じデータに行った解析の結果を示す
## pA B - A C - A C - B pmin
## 1 0.5197138 0.9469076 0.7045524 0.5081796 0.5081796
## 2 0.1818622 0.1630603 0.4645053 0.7913402 0.1630603
## 3 0.5285201 0.6565696 0.9814569 0.5409426 0.5409426
## 4 0.5604227 0.9329313 0.7608621 0.5409838 0.5409838
## 5 0.4936425 0.4930665 0.6556916 0.9632016 0.4930665
## 6 0.5232869 0.8040279 0.4921099 0.8665841 0.4921099
p2 <- data.frame(p = c(p$pA, p$pmin), g = c(rep("ANOVA",10000),rep("Tukey",10000)))
histn(p2$p ~ p2$g, xlab = "P value", ylab = "頻度") #図1の描画
sum(p$pA < 0.05) #分散分析の危険率
## [1] 494
sum(p$pmin < 0.05) #Tukey-Kramer法の危険率
## [1] 484
------------------------------------------------------
図1 分散分析およびTukey-Kramer法のp値の分布。橙は分散分析、青はTukey-Kramer法を示す。
とりあえず単純に比較した点においては、ともに危険率は5%程度であり、同様の傾向がみられる。だが、ここで問題なのは互いの包含関係である。つまり、p値が分散分析で5%以下なら、Tukey-Kramer法でも5%以下というような関係が見られないと、検定の多重性に引っかかってしまう。そこで、Venn図を作成することで互いの包含関係を確認しよう。
------------------------------------------------------
sigA <- rownames(p)[p$pA < 0.05] #分散分析でp < 0.05となる列名
sigT <- rownames(p)[p$pmin < 0.05] #Tukey-Kramer法でp < 0.05となる列名
sigList = list(ANOVA = sigA, Tukey = sigT)
venn.diagram(sigList, filename = "Fig_error.png", fill = c(2,4),
alpha = 0.4, disable.logging = T) #図2の図示。ファイルが出力される
#上記の関数でsigAとsigTの包含関係を図示
#例えばsigA = {1,2,4}、sigT = {3,4,5}なら、sigA固有の要素が2つ、sigT固有の要素が2つ、sigAとsigTの共有要素が1つあるようなVenn図が描かれる。
------------------------------------------------------
図2 分散分析およびTukey-Kramer法でp < 0.05となった標本群の包含関係
図2を見ると、それぞれの解析で約500回、p < 0.05となっていたのだが、それらは独立ではなく、分散分析で有意となる場合はTukey-Kramer法でも有意となっていることがわかる。ここが共通しているということは、逆に分散分析で有意ではないなら、Tukey-Kramer法でも有意にならないということだ。一部は片方の解析だけで有意になる場合があるようだが、それは10000回のシミュレートのうち(45 + 35)/10000 < 0.01であることもわかり、ほとんど起こらないことがわかる。以上から、分散分析とそれを足掛かりにしたTukey-Kramer法は検定の多重性を生じないことがわかる。
Tukey-Kramer法と分散分析の整合性ー検出力
続いて、Tukey-Kramer法の検出力を確認しよう。分散分析の時と同様、A = B = Cに対する対立仮説は、どれかの標本の平均値が異なっている、となる。3標本の場合、平均値を大きいものからA、B、Cとしたとき、（1）A = B > C、（2）A > B = C、（3）A > B > Cである。まずは、A = B > Cの場合だ。
------------------------------------------------------
nA <- 30
nB <- 30
nC <- 30
meA <- 70
meB <- 70
meC <- 50
s <- 20
pA <- NULL
pM <- NULL
for (i in 1:10000){
mA <- rnorm(nA, mean = meA, sd = s)
mB <- rnorm(nB, mean = meB, sd = s)
mC <- rnorm(nC, mean = meC, sd = s)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
d$g <- as.factor(d$g)
fit1 <- aov(d ~ g, data = d)
fit2 <- summary(fit1)
fit3 <- glht(fit1, linfct = mcp(g = "Tukey"))
fit4 <- summary(fit3)
pA <- c(pA, fit2[[1]]$`Pr(>F)`[1])
pM <- rbind(pM, fit4$test$pvalues[1:3])
}
pM <- as.data.frame(pM)
colnames(pM) <- names(fit4$test$tstat)
pM$pmin <- apply(pM, MARGIN = 1, FUN = min)
p <- cbind(pA, pM)
head(p)
## pA B - A C - A C - B pmin
## 1 3.744175e-04 0.9759637 1.067775e-03 2.129108e-03 1.067775e-03
## 2 5.391683e-05 0.1173968 3.298466e-05 2.390391e-02 3.298466e-05
## 3 2.062436e-02 0.9937609 3.590441e-02 4.661819e-02 3.590441e-02
## 4 7.472923e-06 0.7624063 2.801156e-04 1.821244e-05 1.821244e-05
## 5 1.456008e-03 0.9119229 9.105122e-03 2.577581e-03 2.577581e-03
## 6 1.147745e-03 0.7126148 1.300451e-03 1.480138e-02 1.300451e-03
p2 <- data.frame(p = c(p$pA, p$pmin), g = c(rep("ANOVA",10000),rep("Tukey",10000)))
histn(p2$p ~ p2$g, xlab = "P value", ylab = "頻度") #図3の描画
sum(p$pA < 0.05)
## [1] 9829
sum(p$pmin < 0.05)
## [1] 9797
------------------------------------------------------
図3 分散分析およびTukey-Kramer法のp値の分布。橙は分散分析、青はTukey-Kramer法を示す。A = B > Cの場合。
危険率と同様に検出力も単純に比較した点においては、同様の傾向がみられる。続いて、Venn図を作成して、互いの包含関係を確認しよう。
------------------------------------------------------
sigA <- rownames(p)[p$pA < 0.05]
sigT <- rownames(p)[p$pmin < 0.05]
sigList = list(ANOVA = sigA, Tukey = sigT)
venn.diagram(sigList filename = "Fig_power1.png", fill = c(2,4),
alpha = 0.4, disable.logging = T) #図4の図示。
------------------------------------------------------
図4 分散分析およびTukey-Kramer法でp < 0.05となった標本群の包含関係。A = B > Cの場合。
やはり、これも危険率と同様の傾向となる。つまり、分散分析の結果とTukey-Kramer法の結果は一致する。では、次にA > B = Cを検討しよう。
------------------------------------------------------
nA <- 30
nB <- 30
nC <- 30
meA <- 70
meB <- 50
meC <- 50
s <- 20
pA <- NULL
pM <- NULL
for (i in 1:10000){
mA <- rnorm(nA, mean = meA, sd = s)
mB <- rnorm(nB, mean = meB, sd = s)
mC <- rnorm(nC, mean = meC, sd = s)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
d$g <- as.factor(d$g)
fit1 <- aov(d ~ g, data = d)
fit2 <- summary(fit1)
fit3 <- glht(fit1, linfct = mcp(g = "Tukey"))
fit4 <- summary(fit3)
pA <- c(pA, fit2[[1]]$`Pr(>F)`[1])
pM <- rbind(pM, fit4$test$pvalues[1:3])
}
pM <- as.data.frame(pM)
colnames(pM) <- names(fit4$test$tstat)
pM$pmin <- apply(pM, MARGIN = 1, FUN = min)
p <- cbind(pA, pM)
head(p)
## pA B - A C - A C - B pmin
## 1 3.892853e-06 6.478990e-06 2.233392e-04 0.66267628 6.478990e-06
## 2 1.975085e-04 5.858637e-04 1.406601e-03 0.96538681 5.858637e-04
## 3 2.691033e-07 1.960839e-06 5.111391e-06 0.97636028 1.960839e-06
## 4 1.610996e-04 9.451706e-02 9.715830e-05 0.06233767 9.715830e-05
## 5 1.752000e-05 1.997538e-05 1.898574e-03 0.40633554 1.997538e-05
## 6 1.909180e-03 8.632149e-02 1.303981e-03 0.29433647 1.303981e-03
p2 <- data.frame(p = c(p$pA, p$pmin), g = c(rep("ANOVA",10000),rep("Tukey",10000)))
histn(p2$p ~ p2$g, xlab = "P value", ylab = "頻度") #図5の描画
sum(p$pA < 0.05)
## [1] 9814
sum(p$pmin < 0.05)
## [1] 9801
------------------------------------------------------
図5 分散分析およびTukey-Kramer法のp値の分布。橙は分散分析、青はTukey-Kramer法を示す。A > B = Cの場合。
続いて、Venn図を作成する。
------------------------------------------------------
sigA <- rownames(p)[p$pA < 0.05]
sigT <- rownames(p)[p$pmin < 0.05]
sigList = list(ANOVA = sigA, Tukey = sigT)
venn.diagram(sigList filename = "Fig_power2.png", fill = c(2,4),
alpha = 0.4, disable.logging = T) #図4の図示。
------------------------------------------------------
図6 分散分析およびTukey-Kramer法でp < 0.05となった標本群の包含関係。A > B = Cの場合。
最後に、次にA > B > Cを検討しよう。
------------------------------------------------------
nA <- 30
nB <- 30
nC <- 30
meA <- 70
meB <- 60
meC <- 50
s <- 20
pA <- NULL
pM <- NULL
for (i in 1:10000){
mA <- rnorm(nA, mean = meA, sd = s)
mB <- rnorm(nB, mean = meB, sd = s)
mC <- rnorm(nC, mean = meC, sd = s)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
d$g <- as.factor(d$g)
fit1 <- aov(d ~ g, data = d)
fit2 <- summary(fit1)
fit3 <- glht(fit1, linfct = mcp(g = "Tukey"))
fit4 <- summary(fit3)
pA <- c(pA, fit2[[1]]$`Pr(>F)`[1])
pM <- rbind(pM, fit4$test$pvalues[1:3])
}
pM <- as.data.frame(pM)
colnames(pM) <- names(fit4$test$tstat)
pM$pmin <- apply(pM, MARGIN = 1, FUN = min)
p <- cbind(pA, pM)
head(p)
## pA B - A C - A C - B pmin
## 1 3.744175e-04 0.9759637 1.067775e-03 2.129108e-03 1.067775e-03
## 2 5.391683e-05 0.1173968 3.298466e-05 2.390391e-02 3.298466e-05
## 3 2.062436e-02 0.9937609 3.590441e-02 4.661819e-02 3.590441e-02
## 4 7.472923e-06 0.7624063 2.801156e-04 1.821244e-05 1.821244e-05
## 5 1.456008e-03 0.9119229 9.105122e-03 2.577581e-03 2.577581e-03
## 6 1.147745e-03 0.7126148 1.300451e-03 1.480138e-02 1.300451e-03
p2 <- data.frame(p = c(p$pA, p$pmin), g = c(rep("ANOVA",10000),rep("Tukey",10000)))
histn(p2$p ~ p2$g, xlab = "P value", ylab = "頻度") #図7の描画
sum(p$pA < 0.05)
## [1] 9829
sum(p$pmin < 0.05)
## [1] 9797
------------------------------------------------------
図7 分散分析およびTukey-Kramer法のp値の分布。橙は分散分析、青はTukey-Kramer法を示す。A > B > Cの場合。
次が最後のVenn図である。
------------------------------------------------------
sigA <- rownames(p)[p$pA < 0.05]
sigT <- rownames(p)[p$pmin < 0.05]
sigList = list(ANOVA = sigA, Tukey = sigT)
venn.diagram(sigList filename = "Fig_power3.png", fill = c(2,4),
alpha = 0.4, disable.logging = T) #図8の図示。
------------------------------------------------------
図8 分散分析およびTukey-Kramer法でp < 0.05となった標本群の包含関係。A > B > Cの場合。
以上で、いずれの対立仮説であっても分散分析とTukey-Kramer法の解析結果は一致することが確認できた。これで検定の多重性問題、分散分析と長い解説を経て、3標本以上における妥当な比較方法を手に入れることができた。
Rで行うTukey-Kramer法
上記のような背景のもと、Tukey-Kramer法を行ってみよう。今回は分散分析の時に解析をした標本について、どの標本ペアに有意差があったのかをTukey-Kramer法によって明らかにする。条件設定を改めて述べよう。以下のようなサンプルサイズ30の標本が3つあるとする。ともに同じ標準偏差の正規分布からデータを生成していると考えてもらってよい。有意水準は5%としよう。まずは、解析の前にデータの様子を確認だ。
------------------------------------------------------
mA <- c(79, 66, 47, 17, 78, 41, 83, 91, 9, 43, 64, 41, 58, 53, 62,
55, 30, 46, 64, 34, 26, 47, 55, 23, 51, 51, 26, 67, 40, 70)
mB <- c(37, 64, 70, 41, 43, 38, 46, 62, 21, 64, 31, 73, 99, 46, 29,
53, 19, 68, 70, 51, 103, 30, 61, 7, 46, 69, 40, 39, 53, 49)
mC <- c(60, 50, 80, 57, 59, 84, 78, 50, 51, 68, 68, 58, 69, 104, 61,
74, 69, 99, 46, 91, 60, 75, 82, 65, 47, 81, 91, 70, 71, 63)
nA <- length(mA)
nB <- length(mB)
nC <- length(mC)
dfT <- nA + nB + nC - 1
dfF <- 3 - 1
dfR <- dfT - dfF
dfT #全体の自由度
## [1] 89
dfF #要因の自由度
## [1] 2
dfR #残差の自由度
## [1] 87
ST <- sum((c(mA,mB,mC) - mT)^2)
SF <- sum(((m_mA - mT)^2)*nA, ((m_mB - mT)^2)*nB, ((m_mC - mT)^2)*nC)
sT <- sum((c(mA,mB,mC) - mT)^2)/dfT
sF <- sum(((m_mA - mT)^2)*nA, ((m_mB - mT)^2)*nB, ((m_mC - mT)^2)*nC)/dfF
SR <- (sT*dfT - sF*dfF)
sR <- (sT*dfT - sF*dfF)/dfR
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
boxplotn(d$d ~ d$g, ylab = "data", jitter.method = "swarm") #図9の描画
------------------------------------------------------
図9 標本A、B、Cの描画
これを見ると標本Cでは平均値が大きそうな様子がわかる。では、実際に違いそうかを解析してみよう。
------------------------------------------------------
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
d$g <- as.factor(d$g)
fit1 <- aov(d ~ g, data = d) #分散分析
fit2 <- summary(fit1)
fit3 <- TukeyHSD(fit1) #上記の方法とは異なるTukey-Kramer法。上記の方法は下記に改めて行う。
fit2
## Df Sum Sq Mean Sq F value Pr(>F)
## g 2 7007 3503 9.617 0.000168 ***
## Residuals 87 31692 364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit3
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = d ~ g, data = d)
##
## $g
## diff lwr upr p adj
## B-A 0.1666667 -11.584055 11.91739 0.9993696
## C-A 18.8000000 7.049278 30.55072 0.0007371
## C-B 18.6333333 6.882611 30.38406 0.0008275
------------------------------------------------------
まず、aov()関数で分散分析を行い、fit1に格納している。このfit1の結果を用いてTukeyHSD()関数を使ってTukey-Kramer法を行っている。fit3を入力して中身を確認してみよう。一番左列に比較のペアが書かれており、例えば、B-AはBとAのペアを示し、次の列の diffはBからAを引いた値を示している。以下のように計算する。
------------------------------------------------------
mT <- mean(c(mA,mB,mC))
m_mA <- mean(mA)
m_mB <- mean(mB)
m_mC <- mean(mC)
#diff
m_mB - m_mA #B - A
## [1] 0.1666667
m_mC - m_mA #C - A
## [1] 18.8
m_mC - m_mB #C - B
## [1] 18.63333
------------------------------------------------------
次のlwrとuprはスチューデント化範囲分布における95%信頼区間を示しており、以下のように計算する。
------------------------------------------------------
#lwr-upr
Qq <- qtukey(p = 0.95, nmeans = 3, df = dfR)
m_mB - m_mA + Qq * sqrt((sR/2)*(1/nA + 1/nB)) #B - Aのupr
## [1] 11.91739
m_mB - m_mA - Qq * sqrt((sR/2)*(1/nA + 1/nB)) #B - Aのlwr
## [1] -11.58406
m_mC - m_mA + Qq * sqrt((sR/2)*(1/nC + 1/nA)) #C - Aのupr
## [1] 30.55072
m_mC - m_mA - Qq * sqrt((sR/2)*(1/nC + 1/nA)) #C - Aのlwr
## [1] 7.049278
m_mC - m_mB + Qq * sqrt((sR/2)*(1/nB + 1/nC)) #C - Bのupr
## [1] 30.38406
m_mC - m_mB - Qq * sqrt((sR/2)*(1/nB + 1/nC)) #C - Bのlwr
## [1] 6.882611
------------------------------------------------------
最後にP adjはいわゆるp値のことで、スチューデント化範囲分布を利用して以下のように計算できる。
------------------------------------------------------
qAB <- (m_mB - m_mA)/sqrt((sR/2)*(1/nB + 1/nA)) #B - Aのq値
qCA <- (m_mC - m_mA)/sqrt((sR/2)*(1/nA + 1/nC)) #C - Aのq値
qBC <- (m_mC - m_mB)/sqrt((sR/2)*(1/nC + 1/nB)) #C - Bのq値
ptukey(q = qAB, df = dfR, nmeans = 3, lower.tail = F) #B - Aのp値
## [1] 0.9993696
ptukey(q = qCA, df = dfR, nmeans = 3, lower.tail = F) #C - Aのp値
## [1] 0.0007371075
ptukey(q = qBC, df = dfR, nmeans = 3, lower.tail = F) #C - Bのp値
## [1] 0.0008274895
------------------------------------------------------
BCペアとCAペアでp値が有意水準を下回ったので、標本BとC、CとAの間には平均値に違いがありそうだ、と結論できる。一方で、標本AとBに関しては帰無仮説を棄却できなさそうだ、と言える。これは図9の状態と一致しており、解析は成功だといえよう。
さて、別の方法でもTukey-Kramer法をおこなうことができる。multcompパッケージのglht()（general linear hypothesis test）を使う方法で、こちらの方が解析の汎用性が高いので、おすすめをする。10000回シミュレーションした時の方法である。
------------------------------------------------------
fit4 <- glht(fit1, linfct = mcp(g = "Tukey")) #Tukey-Kramer法
fit5 <- summary(fit4)
fit5
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = d ~ g, data = d)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 0.1667 4.9280 0.034 0.999370
## C - A == 0 18.8000 4.9280 3.815 0.000736 ***
## C - B == 0 18.6333 4.9280 3.781 0.000815 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
------------------------------------------------------
Estimateは、平均値の違いを示しており、上記で計算したものと一致する。次列のStd.Errorは次のように計算する。上記のq値を計算した時とちょっと式の形が違うので注意しよう（なぜ違うのかは謎）。今回はすべてサンプルサイズが同じため、同じ値を示す。
------------------------------------------------------
sqrt((sR)*(1/nB + 1/nA)) #B - AのSE
## 4.928003
sqrt((sR)*(1/nC + 1/nB)) #C - BのSE
## 4.928003
sqrt((sR)*(1/nA + 1/nC)) #C - AのSE
## 4.928003
------------------------------------------------------
t valueは実はt分布のt値ことのであり、上記のStd. Errorの計算がtの計算と同じであることに由来している。以下のように計算する。
------------------------------------------------------
tAB <- (m_mB - m_mA)/sqrt((sR)*(1/nB + 1/nA))
tBC <- (m_mC - m_mB)/sqrt((sR)*(1/nC + 1/nB))
tCA <- (m_mC - m_mA)/sqrt((sR)*(1/nA + 1/nC))
tAB
## [1] 0.03382033
tBC
## [1] 3.781113
tCA
## [1] 3.814933
------------------------------------------------------
最後のPr(>|t|)はいわゆるp値のことで、上記のTukeyHSD()と近い値であるものの、微妙に違う値となっている。本当は計算方法を紹介したかったが、突き止められなかった。TukeyHSD()とかなり近い値なので実用上は問題ないが、少し気持ち悪い気もする。しかし、glht()はさまざまな解析に対応しており、今後も活用するのでここで紹介した。TukeyHSD()にできないことが次のようなことだ。
------------------------------------------------------
cld(fit5)
## A B C
## "a" "a" "b"
------------------------------------------------------
cldはcompact letter displayの略称で、glht()の解析結果をcld()に通すと、標本A、B、Cに対応するアルファベットが生成される。このアルファベットは平均値の違いを示しており、同じアルファベット間では有意差がなく、異なるアルファベット間では有意差があることを示している。つまり、今回は標本Cだけが平均値が異なることを示している。このような関数に対応しているので、非常にglht()コマンドは便利である。このように、3標本以上の平均値の比較を行うことができるようになり、解析の幅も広がってきた。次は別の代表的な多重比較法を紹介する。