Tukey-Kramer法

目次

・分散分析ではどの標本間に差があるかまではわからない

・Tukey-Kramer法の検定統計量q

・Tukey-Kramer法と分散分析の整合性ー危険率

・Tukey-Kramer法と分散分析の整合性ー検出力

・Rで行うTukey-Kramer法

・補足1：分散分析後のTukey-Kramer法は多重検定にならないの？

・補足2：検定統計量qの分母の2ってなに？

・補足3：Games-Howell法―不等分散のときの多重比較―

分散分析ではどの標本間に差があるかまではわからない

　分散分析を使うことで、3標本の平均値の間の差を、第1種の過誤を増大させることなく検出できた。しかし、どの標本間に平均値の差があるかまでは分散分析ではわからないのであった。3標本の平均値は異なりました、だけでなく、どの標本に差がありました、と多くの場合は言いたいはずだ。そこで多数の標本の間で、どの標本の間に差があったかを解析する手段が、多重比較multiple comparisonと呼ばれる方法になる。t検定を複数回行う方法も多重比較の一種となるが、以前も言った通り、そのままではfamiliywise type I error (FWER) （= 帰無仮説が正しいときに帰無仮説族の少なくとも1つが誤って棄却される確率） が大きくなるので問題がある。検定を複数回行うときに使える工夫は別の機会に紹介するとして、今回は別の多重比較法を紹介する。その一つがTukey-Kramer（テューキー・クレーマー）法（あるいはサンプルサイズがそろっている場合はTukeyのhonestly significant difference (HSD)、Tukeyの範囲検定）である。Tukey-Kramer法は、すべての標本が等分散であるときに使える。分散分析後に、Tukey-Kramer法の検定をすることでどの標本の間に違いがあるかを検出するのである（本当は分散分析は不要であるが、Rの仕様上、分散分析の結果を使って解析する）。

（分散分析後のTukey-Kramer法の妥当性については補足に移動（2024.7.26））

Tukey-Kramer法の検定統計量q

　この流れもおなじみになってきた。Tukey-Kramer法の検定統計量はqと呼ばれる。このqを、スチューデント化範囲分布Studentized range distributionのもとで、検定を行うのである。けったいな名前な分布名であるが、考え方は今までとそれほど変わらない。この検定統計量qであるが以下のように計算する。

ここで3標本以上あることを前提として、2標本AとBを取り出したとする。XAとXBは各標本の平均、nAとnBは各標本のサンプルサイズである。そして、(sR)^2は分散分析で登場した残差の平均平方和である。この式の形は、分母の分散推定値が2で割られていることを除けば、検定統計量Studentのtとほとんど同じ形だ（なぜこの「2」が登場するのかは、補足2参照）。

　帰無仮説が正しいとき、無作為に2標本を選んで統計量qを計算するとそれはt分布と変わらない。帰無仮説が正しいときの最大の平均を持つ標本と最小の平均を持つ標本で統計量qを計算した時にスチューデント化範囲分布となり、t分布の拡張となっている。スチューデント化範囲分布では、すべての標本を使っているのに、上記の検定統計量qでは2標本だけで計算している点が異なっている。これは、すべての標本が等分散であることを仮定しているので、分母は基本的にどの2標本を使おうと、すべての標本を使おうと大きな乖離が生じないことを仮定している。この意味で、sR^2はt検定における合併不偏分散と同じである（分散が異なる場合は、Games-Howell法を用いる（補足3））。

　帰無仮説の下でのt分布よりもスチューデント化範囲分布は分布の裾が重く、それゆえに帰無仮説が棄却されるためには、2標本のt検定の時よりもより効果量が大きくなくてはならない。直感的には、スチューデント化範囲分布を利用することで、検定で帰無仮説が棄却されにくくなり、多重検定の問題をクリアしているのである。これらのことから、検定の手順が多重比較の手続きに変更されたというよりも、検定統計量tとt分布の多重比較版が検定統計量qとスチューデント化範囲分布であると考えるほうが妥当だ。

　より詳しく説明すると、スチューデント化範囲分布の利用によってFWERを5%に抑える=正しい帰無仮説をすべて棄却しない確率が95%にすることが可能となる。最も大きいq値をqmaxとすれば、スチューデント化範囲分布において「qmaxが棄却域に存在する」=「対比較のq値の少なくとも一つ棄却域に存在する」ことを意味する。qmaxは最も外れ値になりやすいとはいえ、帰無仮説が正しいと仮定としたとき、qmaxが棄却域に存在する確率は5%以下である。それでも棄却域に存在するなら、それは帰無仮説が正しくない可能性が高いからであると判断できる。見方を変えれば、帰無仮説が正しいときにqmax = 対比較のq値のうち少なくとも1つが棄却域に入ってしまい、いすれかの帰無仮説が棄却させる危険率が5%であるともいえる。この危険率はまさに、FWERを5%に抑えることに他ならない。逆に「qmaxが棄却域に存在しない」=「すべての対比較のq値が棄却域に存在しない」=「すべての対比較に関する帰無仮説族を保留する」ことであり、この確率はqmaxが棄却域外に存在する確率95%であり、整合性が取れている。

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

------------------------------------------------------

とりあえず単純に比較した点においては、ともに危険率は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を見ると、それぞれの解析で約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

------------------------------------------------------

危険率と同様に検出力も単純に比較した点においては、同様の傾向がみられる。続いて、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の図示。

------------------------------------------------------

やはり、これも危険率と同様の傾向となる。つまり、分散分析の結果と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

------------------------------------------------------

続いて、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の図示。

------------------------------------------------------

最後に、次に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

------------------------------------------------------

次が最後の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の図示。

------------------------------------------------------

以上で、いずれの対立仮説であっても分散分析と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の描画

------------------------------------------------------

これを見ると標本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標本以上の平均値の比較を行うことができるようになり、解析の幅も広がってきた。次は別の代表的な多重比較法を紹介する。

補足1：分散分析後のTukey-Kramer法は多重検定にならないの？

　今までさんざん、有意水準の補正がないまま複数回検定するのは悪手だ、と言っておきながら、「分散分析後に多重比較します」はおかしいだろ、と思わないだろうか。一回、考え出すと疑念が膨らむ読者もいるだろう。これは複数回検定の制約の問題に反しないのだろうか？

心配ご無用、今まで述べてきた複数回検定の制約には反しないようになっている。まずは上の図の左を見てほしい。例えば、t検定を3回行う場合、互いの検定結果は残りの検定結果に影響がない。つまり、独立な関係の検定となる。この場合、検定の結果を掛け算できて、1回の検定につき有意にならない確率が0.95なら、3回の検定すべてで有意にならない確率を0.95 × 0.95 × 0.95 = 0.857…と計算できる。そして、以前にも言ったように1からこの値を引くことでFWERが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値を記載することで、その後にほかの研究者がメタ解析で活用できる可能性も広がる。そして、正しく統計を理解している査読者に回ったら、上記のように書いておけば事情をある程度汲んでくれる（理解してない査読者に回ることもあるが、そこに合わせる道理はないはずだ）。間違っても、なんとか有意になるように解析方法をあれこれ変えてはいけない。

*************************************************

補足2：検定統計量qの分母の2ってなに？

　検定統計量qは、分母の2を除けば、検定統計量tとほとんど同じ計算である。この点を考えると、「2」は余分な数字に見えるし、どこから生じたのか意味不明に見える。ここで、スチューデント化範囲分布のqの計算を見てみる。標本がS個存在し、それぞれのサンプルサイズをn1～nSとする。Uは全体の不偏分散、Mは各標本のうち最大の平均、mは最小の平均である。

さて、Tukey-Kramer法においては、対比較を考えるので、2（= S）標本の状況を考えていることがわかる。不偏分散は残差平方和である。つまり、以下の通り。

つまり、2は検定統計量qの比較する標本数に由来する数である。それはそうなのだが、サンプルサイズの等しい状況のTukeyのHSDから、サンプルサイズが異なる状況のTukey-Kramer法に拡張するときに、重要な役割を果たした値でもある（永田・𠮷田（1997）, p37）。もう少し、深堀りしよう。

　ここで、√2を除いた部分の式はStudentのt値そのものだから、以下の通り。

Tukey-Kramer法では、検定統計量をqとせずに、tを用いることも多い（実際、glht関数の出力はtだった）。おそらくはt検定との関連性などから、拡張であることを意図的に示すうえで重要な考え方だったろう。tを起点と考えて、qの分布、すなわちスチューデント化範囲分布の1/√2に従うとして検定統計量の比較を行う。

　さて、TukeyのHSDではサンプルサイズが等しいことを仮定しているから、以下のようになる。

このように、サンプルサイズが等しいことを仮定すると√2はキャンセルされて、式の中から消えてしまう。私たちは拡張されたTukey-Kramer法から学んでしまったから、わかりずらくなっているが、もともとの多重比較法はこのq値から始まったのだ。ここを起点に考えると、そもそも√2が出てくる余地がない。

　サンプルサイズが異なる場合に拡張するにあたって、サンプルサイズの違いを式内にどう反映するか、算術平均か、幾何平均か、それとも調和平均か、そのどれが妥当かは必ずしも自明ではない。結果として、調和平均として組み込み、その際に2が登場したのである。

補足3：Games-Howell法―不等分散のときの多重比較―

　Tukey-Kramer法では、各標本が等分散であることを仮定している。では、不等分散ならどうするべきだろうか。その時に用いられるのがGames-Howell（ゴメス・ハウエル）法である。Tukey-Kramer法がStudentのt検定の多重比較版であったように、Games-Howell法はまさにWelchのt検定の多重比較版である。検定統計量qは、以下のように計算する。

UA、UBは、選ばれた各標本の不偏分散である。ここでも、2が登場する以外、Welchのtの値そのものである。自由度も、Weichのt検定同様に、Welch–Satterthwaiteの式を用いて近似する。

青木先生が、プログラムを公開しているので、こちらを利用するのもよいだろう。