Kruskal-Wallis検定


標本が正規分布に従わないときの分散分析の代替法

　分散分析は標本の正規性が満たされないと利用できないという制約がある（ただし、ある程度、頑健である）。これまでのt検定とMann-WhitneyのU検定の関係同様に、分散分析にもノンパラメトリックな代替法が存在する。Kruskal-Wallis検定は、標本が正規分布に従わないときに利用が推奨される方法である。これまでのノンパラメトリック検定と同様、標本を順位に変換することで、特に外れ値に対して頑健な結果を返すことができる。Kruskal-Wallis検定の帰無仮説は、分散分析と類似しており、すべての標本の中央値が等しいというものだ。また、分散分析と同じく、Kruskal-Wallis検定によって帰無仮説が棄却されたときの解釈は、すべての標本の中央値が等しいわけではない、というもので、具体的にどの標本間で中央値が異なるかを明らかにしてくれるわけではない。もし、どの標本間で異なるかを検出したいのであれば、複数回検定をしてから有意水準をHolmの補正するか、Steel-Dwass検定を行おう。


Kruskal-Wallis検定の危険率

　Kruskal-Wallis検定は下記のような手順に従って検定統計量Hを算出する。

サンプルサイズが大きいとき、検定統計量Hは自由度を標本数 - 1とするカイ二乗分布に従うことを利用して検定（漸近検定）を行う。サンプルサイズが小さいときは、正確なp値計算のために、精度の高いH分布を使い、正確検定を行うことが推奨されている。

　Mann-WhitneyのU検定の時と同じように、標本の分布が二峰分布、コーシー分布の時の、危険率や検出力を確認してみよう。まずは、標本が二峰分布に従うときを考える。Kruskal-Wallis検定とともに分散分析も行い、比較を行う。


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

library(plotn)

library(kSamples)

#install.packages("kSamples", repos="http://cran.ism.ac.jp/")


library(VennDiagram)

#install.packages("VennDiagram", repos="http://cran.ism.ac.jp/")


rbimodal <- function(n, total_mean = 60, mean1 = 45, rate = 1/3, sd1 = 14.5, sd2 = 11){#二峰分布の生成関数

samples <- NULL

for(i in 1:n){

b <- rbinom(1,1,rate)

if(b == 0){

samples <- c(samples, rnorm(1, mean1, sd1))

} else {

samples <- c(samples, rnorm(1, (total_mean - mean1 * (1 - rate))/rate, sd2))

}

}

return(samples)

}


#二峰分布

nA <- 30

nB <- 30

nC <- 30

meA <- 50

meB <- 50

meC <- 50


pA <- NULL

pK <- NULL

Fvalue <- NULL

H <- NULL

for (i in 1:10000){

mA <- rbimodal(nA, total_mean = meA, mean1 = meA - 15)#二峰分布の生成。どんな分布か各自で描画してみよう

mB <- rbimodal(nB, total_mean = meB, mean1 = meB - 15)

mC <- rbimodal(nC, total_mean = meC, mean1 = meC - 15)

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)

#ここからF値の算出

mT <- mean(c(mA,mB,mC))

m_mA <- mean(mA)

m_mB <- mean(mB)

m_mC <- mean(mC)

dfT <- nA + nB + nC - 1

dfF <- 3 - 1

dfR <- dfT - dfF

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)/dfR

Fvalue <- c(Fvalue, sF/sR)

#ここまでF値の算出

pA <- c(pA, fit2[[1]]$`Pr(>F)`[1])

fit3 <- kruskal.test(d ~ g, data = d)

pK <- c(pK, fit3$p.value)

#ここからH値の算出

r <- rank(c(mA, mB, mC))

N <- nA + nB + nC

H <- c(H, 12/(N*(N + 1)) * ((sum(r[1:nA])^2)/nA + (sum(r[nA + 1:nB])^2)/nB + (sum(r[nA + nB + 1:nC])^2)/nC) - 3 * (N + 1))

#ここまでH値の算出

}


p <- data.frame(pA, pK)


p2 <- data.frame(p = c(pA, pK), g = c(rep("ANOVA",10000),rep("KW",10000)))


histn(Fvalue, xlab = "統計量F", ylab = "頻度", freq = F)#図1の描画

xx <- 0:40

yy <- df(xx, df1 = dfF, df2 = dfR)

overdraw(points(xx, yy, type = "l", col = "red"))


histn(H, xlab = "統計量H", ylab = "頻度", freq = F)#図2の描画

xx <- 0:40

yy <- dchisq(xx, df = dfF)

overdraw(points(xx, yy, type = "l", col = "red"))


histn(p2$p ~ p2$g, xlab = "P value", ylab = "頻度")#図3の描画


sigA <- rownames(p)[p$pA < 0.05]

sigK <- rownames(p)[p$pK < 0.05]


sigList = list(ANOVA = sigA, KW = sigK)


venn.diagram(sigList, filename = "Fig_error_bimodal_ANOVA_KW.png", fill = c(2,4),

alpha = 0.4, disable.logging = T)#図4の描画。ディレクトリにファイルが生成


sum(pA < 0.05)#分散分析の危険率

## [1] 514


sum(pK < 0.05)#Kruskal-Wallis検定の危険率

## [1] 507

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

図1 二峰分布に対するF分布。赤線（df(..., df1 = dfF, df2 = dfR)）は帰無仮説が正しいときのF分布

図2 二峰分布に対するH分布。赤線（qchisq(..., df = dfF)）は帰無仮説が正しく、サンプルサイズが大きいときに近似されるカイ二乗分布

図3 二峰分布に対する分散分析（青）とKruskal-Wallis検定（橙）のp値の分布

図4 二峰分布に対する分散分析とKruskal-Wallis検定で棄却された標本の包含関係

標本が形が同じで平均も等しい二峰分布に従うとき、上記のように分散分析とKruskal-Wallis検定のいずれにおいても危険率は5%程度となる。F分布もH分布も理論予測通りの挙動であり、大きな問題は見受けられない。したがって、分散分析は、コーシー分布のような外れ値を含みやすい分布でなければ有効に使える可能性が高い。Mann-WhitneyのU検定の時の議論と同様に、平均値は中心極限定理の影響を受けるため、平均値を使う検定は、サンプルサイズが大きければ、母集団の形に頑健になるのだろうと推測される。また、基本的に分散分析とKruskal-Wallis検定で、帰無仮説が棄却される標本は、危険率5%のなかの約80%が共通である（430/514、430/507）。いずれの検定でも問題ないだろう。

　次に標本がコーシー分布に従うときを考えよう。


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

#コーシー分布

nA <- 30

nB <- 30

nC <- 30

meA <- 50

meB <- 50

meC <- 50


pA <- NULL

pK <- NULL

Fvalue <- NULL

H <- NULL

for (i in 1:10000){

mA <- rcauchy(nA, location = meA) #locationで中央値を指定する

mB <- rcauchy(nB, location = meB)

mC <- rcauchy(nC, location = meC)

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)

#ここからF値の算出

mT <- mean(c(mA,mB,mC))

m_mA <- mean(mA)

m_mB <- mean(mB)

m_mC <- mean(mC)

dfT <- nA + nB + nC - 1

dfF <- 3 - 1

dfR <- dfT - dfF

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)/dfR

Fvalue <- c(Fvalue, sF/sR)

#ここまでF値の算出

pA <- c(pA, fit2[[1]]$`Pr(>F)`[1])

fit3 <- kruskal.test(d ~ g, data = d)

pK <- c(pK, fit3$p.value)

#ここからH値の算出

r <- rank(c(mA, mB, mC))

N <- nA + nB + nC

H <- c(H, 12/(N*(N + 1)) * ((sum(r[1:nA])^2)/nA + (sum(r[nA + 1:nB])^2)/nB + (sum(r[nA + nB + 1:nC])^2)/nC) - 3 * (N + 1))

#ここまでH値の算出

}


p <- data.frame(pA, pK)


p2 <- data.frame(p = c(pA, pK), g = c(rep("ANOVA",10000),rep("KW",10000)))


histn(Fvalue, xlab = "統計量F", ylab = "頻度", freq = F)#図5の描画

xx <- seq(1,40,length= 200)

yy <- df(xx, df1 = dfF, df2 = dfR)

overdraw(points(xx, yy, type = "l", col = "red"))


histn(H, xlab = "統計量H", ylab = "頻度", freq = F)#図6の描画

xx <- seq(1,40,length= 200)

yy <- dchisq(xx, df = dfF)

overdraw(points(xx, yy, type = "l", col = "red"))


histn(p2$p ~ p2$g, xlab = "P value", ylab = "頻度")#図7の描画


sigA <- rownames(p)[p$pA < 0.05]

sigK <- rownames(p)[p$pK < 0.05]


sigList = list(ANOVA = sigA, KW = sigK)


venn.diagram(sigList, filename = "Fig_error_cauchy_ANOVA_KW.png", fill = c(2,4),

alpha = 0.4, disable.logging = T)#図8の描画


sum(pA < 0.05)#分散分析の危険率

## [1] 178


sum(pK < 0.05)#Kruskal-Wallis検定の危険率

## [1] 479

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

図5 コーシー分布に対するF分布。赤線（df(..., df1 = dfF, df2 = dfR)）は帰無仮説が正しいときのF分布

図6 コーシー分布に対するH分布。赤線（qchisq(..., df = dfF)）は帰無仮説が正しく、サンプルサイズが大きいときに近似されるカイ二乗分布

図7 コーシー分布に対する分散分析（青）とKruskal-Wallis検定（橙）のp値の分布

図8 コーシー分布に対する分散分析とKruskal-Wallis検定で棄却された標本の包含関係

一見すれば分散分析の危険率が低くなっており、良い検定になっているようにも見えるが、図5を確認するとF分布から逸脱が見られる。また、図7を見てもp値の分布が不安定になっており、信用にかける。一方で、Kruskal-Wallis検定は危険率5%を維持しつつ、H分布も安定だ。これまでの様々な検定の議論のなかから、気づいてきた人もいるかもしれないが、危険率5%と設定しているにもかかわらず、それを大きく上回る場合、下回る場合、いずれにおいても検定の精度が安定していないことが多い。どんな標本の状態でも危険率5%を維持できることこそが、頑健な検定といえるのである。


Kruskal-Wallis検定の検出力

　続いて検出力について考えよう。Kruskal-Wallis検定の帰無仮説が、すべての標本の中央値が等しい、というものだったので、それに対する3標本での具体的な対立仮説は次の3通りとなる：（1）A = B > C、（2）A > B = C、（3）A > B > C。これは分散分析の時と同様だ。今回は、（2）のA > B = Cの時だけ検討しよう。では、二峰分布の時の検出力は以下のようにコードできる。


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

#二峰分布

nA <- 30

nB <- 30

nC <- 30

meA <- 70

meB <- 50

meC <- 50


pA <- NULL

pK <- NULL

Fvalue <- NULL

H <- NULL

for (i in 1:10000){

mA <- rbimodal(nA, total_mean = meA, mean1 = meA - 15)

mB <- rbimodal(nB, total_mean = meB, mean1 = meB - 15)

mC <- rbimodal(nC, total_mean = meC, mean1 = meC - 15)

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)

#ここからF値の算出

mT <- mean(c(mA,mB,mC))

m_mA <- mean(mA)

m_mB <- mean(mB)

m_mC <- mean(mC)

dfT <- nA + nB + nC - 1

dfF <- 3 - 1

dfR <- dfT - dfF

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)/dfR

Fvalue <- c(Fvalue, sF/sR)

#ここまでF値の算出

pA <- c(pA, fit2[[1]]$`Pr(>F)`[1])

fit3 <- kruskal.test(d ~ g, data = d)

pK <- c(pK, fit3$p.value)

#ここからH値の算出

r <- rank(c(mA, mB, mC))

N <- nA + nB + nC

H <- c(H, 12/(N*(N + 1)) * ((sum(r[1:nA])^2)/nA + (sum(r[nA + 1:nB])^2)/nB + (sum(r[nA + nB + 1:nC])^2)/nC) - 3 * (N + 1))

#ここまでH値の算出

}


p <- data.frame(pA, pK)


p2 <- data.frame(p = c(pA, pK), g = c(rep("ANOVA",10000),rep("KW",10000)))


histn(Fvalue, xlab = "統計量F", ylab = "頻度", freq = F)#図9の描画

xx <- seq(0,40,length= 200)

yy <- df(xx, df1 = dfF, df2 = dfR)

overdraw(points(xx, yy, type = "l", col = "red"))


histn(H, xlab = "統計量H", ylab = "頻度", freq = F)#図10の描画

xx <- seq(0,40,length= 200)

yy <- dchisq(xx, df = dfF)

overdraw(points(xx, yy, type = "l", col = "red"))


histn(p2$p ~ p2$g, xlab = "P value", ylab = "頻度")#図11の描画


sigA <- rownames(p)[p$pA < 0.05]

sigK <- rownames(p)[p$pK < 0.05]


sigList = list(ANOVA = sigA, KW = sigK)


venn.diagram(sigList, filename = "Fig_power_bimodal_ANOVA_KW.png", fill = c(2,4),

alpha = 0.4, disable.logging = T)#図12の描画

sum(pA < 0.05)

## [1] 8912


sum(pK < 0.05)

## [1] 8768

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

図9 平均の異なる二峰分布に対するF分布。赤線（df(..., df1 = dfF, df2 = dfR)）は帰無仮説が正しいときのF分布

図10 平均の異なる二峰分布に対するH分布。赤線（qchisq(..., df = dfF)）は帰無仮説が正しく、サンプルサイズが大きいときに近似されるカイ二乗分布

図11 平均の異なる二峰分布に対する分散分析（青）とKruskal-Wallis検定（橙）のp値の分布

図12 平均の異なる二峰分布に対する分散分析とKruskal-Wallis検定で棄却された標本の包含関係

危険率の時と同様に、二峰分布のとき、分散分析もKruskal-Wallis検定も類似した結果を返す。平均が20程度異なるときの検出力は、ともに87～89%程度であり、また帰無仮説が棄却される標本も重複が大きい。いずれの手法を用いても解釈は一致する可能性が高いだろう。

　続いてコーシー分布の時を考える。


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

#コーシー分布

nA <- 30

nB <- 30

nC <- 30

meA <- 70

meB <- 50

meC <- 50


pA <- NULL

pK <- NULL

Fvalue <- NULL

H <- NULL

for (i in 1:10000){

mA <- rcauchy(nA, location = meA) #locationで中央値を指定する

mB <- rcauchy(nB, location = meB)

mC <- rcauchy(nC, location = meC)

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)

#ここからF値の算出

mT <- mean(c(mA,mB,mC))

m_mA <- mean(mA)

m_mB <- mean(mB)

m_mC <- mean(mC)

dfT <- nA + nB + nC - 1

dfF <- 3 - 1

dfR <- dfT - dfF

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)/dfR

Fvalue <- c(Fvalue, sF/sR)

#ここまでF値の算出

pA <- c(pA, fit2[[1]]$`Pr(>F)`[1])

fit3 <- kruskal.test(d ~ g, data = d)

pK <- c(pK, fit3$p.value)

#ここからH値の算出

r <- rank(c(mA, mB, mC))

N <- nA + nB + nC

H <- c(H, 12/(N*(N + 1)) * ((sum(r[1:nA])^2)/nA + (sum(r[nA + 1:nB])^2)/nB + (sum(r[nA + nB + 1:nC])^2)/nC) - 3 * (N + 1))

#ここまでH値の算出

}


p <- data.frame(pA, pK)


p2 <- data.frame(p = c(pA, pK), g = c(rep("ANOVA",10000),rep("KW",10000)))


histn(Fvalue, xlab = "統計量F", ylab = "頻度", freq = F)#図13の描画

xx <- seq(0,40,length= 200)

yy <- df(xx, df1 = dfF, df2 = dfR)

overdraw(points(xx, yy, type = "l", col = "red"))


histn(H, xlab = "統計量H", ylab = "頻度", freq = F)#図14の描画

xx <- seq(0,40,length= 200)

yy <- dchisq(xx, df = dfF)

overdraw(points(xx, yy, type = "l", col = "red"))


histn(p2$p ~ p2$g, xlab = "P value", ylab = "頻度")#図15の描画


sigA <- rownames(p)[p$pA < 0.05]

sigK <- rownames(p)[p$pK < 0.05]


sigList = list(ANOVA = sigA, KW = sigK)


venn.diagram(sigList, filename = "Fig_power_cauchy_ANOVA_KW.png", fill = c(2,4),

alpha = 0.4, disable.logging = T)#図16の描画


sum(pA < 0.05)

## [1] 8472


sum(pK < 0.05)

## [1] 10000

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

図13 中央値の異なるコーシー分布に対するF分布。赤線（df(..., df1 = dfF, df2 = dfR)）は帰無仮説が正しいときのF分布

図14 中央値の異なるコーシー分布に対するH分布。赤線（qchisq(..., df = dfF)）は帰無仮説が正しく、サンプルサイズが大きいときに近似されるカイ二乗分布

図15 中央値の異なるコーシー分布に対する分散分析（青）とKruskal-Wallis検定（橙）のp値の分布

図16 中央値の異なるコーシー分布に対する分散分析とKruskal-Wallis検定で棄却された標本の包含関係

上記の結果を見て、私から正直なコメントをするとすれば、どちらがよりよい分布の形となっているかは判断がつかない。傾向としてF分布は極めてばらつきが大きく、安定していないようすが見受けられる。H分布は連続分布となっているかは怪しいところだ。今後私の勉強で、平均値の差から予測される分布を描けるようになったら、評価も変わるだろうか。とりあえず、検出力だけで見ると、Kruskal-Wallis検定のほうが高いことが変わる。少なくとも、危険率に関しては、分散分析は安定した結果を出せないので、外れ値を含む標本ではKruskal-Wallis検定を行うのが賢明であろう。


Rで行うKruskal-Wallis検定

　上記のような背景のもと、Kruskal-Wallis検定を行ってみよう。下記のようなサンプルサイズ30の標本が3つあるとする。分布を描いてみるとわかるが、外れ値が存在するデータだ。これを有意水準5%でどれかの標本に有意差が存在するかを確かめよう。


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

mA <- c(69.157, 70.189, 69.962, 70.194, 71.362,

70.836, 67.288, -76.316, 69.832, 70.18,

69.573, 69.173, 71.125, 73.144, 69.888,

69.698, 69.928, 66.933, 68.538, 70.91,

56.037, 71.279, 75.24, 70.903, 66.965,

68.856, 69.071, 68.595, 68.79, 70.52)

mB <- c(50.35, 50.421, 50.805, 28.743, 48.35,

48.471, 87.512, 50.211, 58.163, 50.148,

42.519, 33.793, 47.016, 49.605, 50.961,

49.524, 48.681, 47.495, 50.952, 51.26,

43.279, 98.529, 50.167, 52.348, 50.54,

49.673, 49.782, 50.882, 49.113, 48.29)

mC <- c(50.85, 50.099, 49.944, 50.081, 51.522,

50.293, 50.826, 50.256, 47.586, 51.457,

46.923, 48.748, 49.954, 50.587, 48.736,

46.086, 49.405, 50.473, 50.885, 51.313,

52.406, 51.677, 50.157, 47.842, 45.6,

53.415, 51.406, 49.881, 50.346, 48.097)


nA <- length(mA)

nB <- length(mB)

nC <- length(mC)


d <- data.frame(d = c(mA, mB, mC),

g = c(rep("mA", nA),

rep("mB", nB),

rep("mC", nC))

)


boxplotn(d$d ~ d$g, ylab = "data", jitter.method = "swarm")#図17の描画

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

図17 標本A、B、Cの描画。外れ値の存在に注意しよう。

標本Aは残り2つよりも中央値が大きい傾向があるが、一方で残り2標本よりも小さな外れ値が存在する。おそらく、平均値を検定する分散分析では標本Aの分布の中心を過小評価してしまうだろう。こんなときはKruskal-Wallis検定が有効だ。下記のコマンドで実行する。


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

fit <- kruskal.test(d ~ g, data = d)


fit

##

## Kruskal-Wallis rank sum test

##

## data: d by g

## Kruskal-Wallis chi-squared = 44.719, df = 2, p-value = 1.947e-10

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


Kruskal-Wallis chi-squaredは上記で議論してきた検定統計量Hのことで、以下のように計算できる。（なお、同順位があるときは計算が一致しない、どう計算するんだ……？）


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

r <- rank(c(mA, mB, mC))

N <- nA + nB + nC

H0 <- 12/(N*(N + 1)) * ((sum(r[1:nA])^2)/nA + (sum(r[nA + 1:nB])^2)/nB + (sum(r[nA + nB + 1:nC])^2)/nC) - 3 * (N + 1)


H0

## [1] 44.71922

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


dfは自由度のことで、これは標本数 - 1で計算する。今回は3標本あるので、2となる。最後、p-valueは、統計量Hがdfに従うカイ二乗分布に従うと仮定して算出されており、下記のように計算できる。今回は、有意水準5%を下回っているので、標本のどれかは中央値が異なりそうだ。


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

pchisq(H0, df = 2, lower.tail = F)

## [1] 1.946906e-10

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


サンプルサイズが十分に大きいときは、上記のカイ二乗分布への近似が成り立つため、p値の計算にカイ二乗分布が使える。とりあえず、下記のコードを実行することで帰無仮説が正しいときに、近似が成り立つことを示す。


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

nA <- 30

nB <- 30

nC <- 30

H <- NULL

for(i in 1:300000){

sA <- rnorm(nA, mean = 0, sd = 20)

sB <- rnorm(nB, mean = 0, sd = 20)

sC <- rnorm(nC, mean = 0, sd = 20)

r <- rank(c(sA, sB, sC))

N <- nA + nB + nC

H <- c(H, 12/(N*(N + 1)) * ((sum(r[1:nA])^2)/nA + (sum(r[nA + 1:nB])^2)/nB + (sum(r[nA + nB + 1:nC])^2)/nC) - 3 * (N + 1))

}


histn(H, freq = F, xlab = "統計量H")#図18の描画

xx <- 0:40

yy <- dchisq(0:40, df = 2)

overdraw(points(xx, yy, type = "l", col = "red"))

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

図18 帰無仮説が正しいときの統計量Hの分布。赤線は自由度2のカイ二乗分布。

たしかに近似の精度が非常に良いことがわかる。この近似を利用しているため、kruskal.test()関数のデフォルトは漸近検定である。一方で、サンプルサイズが小さいときは近似が成立しないため、統計量Hの正確な分布が必要になる。Kruskal-Wallis検定の正確検定を行うときは、kSamplesパッケージのqn.test()関数を使おう。この関数はシミュレーションによって正確なH分布を作成する。


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

mA <- c(52, 58, 54, 53, 43)

mB <- c(49, 49, 50, 51, 49)

mC <- c(51, 50, 49, 49, 48)


nA <- length(mA)

nB <- length(mB)

nC <- length(mC)


d <- data.frame(d = c(mA, mB, mC),

g = c(rep("mA", nA),

rep("mB", nB),

rep("mC", nC))

)


boxplotn(d$d ~ d$g, ylab = "data", jitter.method = "swarm")#図19の描画

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

図19 標本A、B、Cの描画

今回はサンプルサイズ5である。中央値は異なりそうだが、果たしてどうだろうか？


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

fit1 <- kruskal.test(d ~ g, data = d)

fit1

##

## Kruskal-Wallis rank sum test

##

## data: d by g

## Kruskal-Wallis chi-squared = 3.5599, df = 2, p-value = 0.1687


fit2 <- qn.test(mA, mB, mC, test="KW", Nsim = 1000000, method="exact")

fit2

## Kruskal-Wallis k-sample test.

##

## Number of samples: 3

## Sample sizes: 5, 5, 5

## Number of ties: 6

##

## Null Hypothesis: All samples come from a common population.

##

## test statistic asympt. P-value exact P-Value

## 3.5600 0.1687 0.1740

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


特にqn.test()の出力を確認すると、asympt. P-valueとあるのが、カイ二乗近似を用いたときのp値である。これはkruskal.test()のp値と一致している。exact P-Valueが正確検定によるp値でわずかに違うことがわかる。小標本の時は可能であれば、qn.test()を使うとよいだろう。ただし、サンプルサイズが大きいときは、シミュレーションが重たくなり、パソコンがスタックする可能性があるので、qn.test()でmethod = "exact"を用いるのは避けたほうが良い。今回は、5%水準で有意差があると結論できないことになる。