Steel-Dwass検定

標本が正規分布に従わないときのTukey-Kramer法の代替法

　分散分析と同じく、Tukey-Kramer法は標本の正規性が満たされないと利用できないという制約がある。ただし、後述するように、これまた分散分析と同様に、正規性がなくてもある程度頑健な結果を返し、正規分布に従うか否かよりも外れ値に対して安定した結果が出ない。外れ値に対して頑健な結果を出してくれる多重比較法が、Steel-Dwass検定と呼ばれる方法である。これまでのノンパラメトリック検定と同様、標本を順位に変換して検定を行う。Steel-Dwass検定の帰無仮説は、Tukey-Kramer法と類似しており、取り出した2標本の中央値が等しいというものだ。Kruskal-Wallis検定と異なり、具体的にどの標本間で中央値が異なるかを明らかにすることができる。

Steel-Dwass検定の危険率

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

サンプルサイズが大きいとき、検定統計量SDを標本数 - 1、自由度∞のスチューデント化範囲分布の下で検定（漸近検定）を行う。Tukey-Kramer法の時と同様、帰無仮説のもとでSDは厳密にスチューデント化範囲分布に従うわけではなく、SD分布よりもスチューデント化範囲分布のほうが分布の裾が広い。それゆえ、有意水準5%であっても、2標本の検定時よりもより大きな効果量がないと、帰無仮説は棄却されない（後述する）。サンプルサイズが小さいときは、正確なp値計算のために、精度の高いSD分布を使い、正確検定を行うことが推奨されている。

　では、標本の分布が二峰分布、コーシー分布の時の、危険率や検出力を確認してみよう。まずは、標本が二峰分布に従うときを考える。Steel-Dwass検定とともにTukey-Kramer法も行い、比較を行う。Steel-Dwass検定を行うNSM3パッケージが提供されているが、計算時間がかかるので、より簡便に使える青木先生のプログラムを利用するのが良いだろう。

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

library(plotn)

library(multcomp)

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

library(NSM3)

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

library(VennDiagram)

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

source("http://aoki2.si.gunma-u.ac.jp/R/src/Steel-Dwass.R", encoding="euc-jp")#青木先生のSteel Dwass検定用プログラム

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)

}

dtukey <- function(q, nmeans, df, nranges = 1){#スチューデント化範囲分布の描画のために定義

  d <- ptukey(q = q, nmeans = nmeans, 

              df = df, nranges = nranges, 

              lower.tail = TRUE, log.p = FALSE)

  d <- d - c(0, d[1:(length(q)-1)])

  return(d)

}

#A = B = C

#二峰分布

nA <- 30

nB <- 30

nC <- 30

meA <- 50

meB <- 50

meC <- 50

pT <- NULL

pS <- 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 <- glht(aov(d ~ g, data = d), linfct = mcp(g = "Tukey")) #Tukey-Kramer法

  fit2 <- summary(fit1)

  pT <- rbind(pT, fit2$test$pvalues[1:3])

  fit3 <- Steel.Dwass(d$d, as.numeric(d$g)) #Steel Dwass検定

  pS <- rbind(pS, fit3[,2])

}

pT <- as.data.frame(pT)

colnames(pT) <- paste0(names(fit2$test$tstat), "_Tukey")

pT_min <- apply(pT, MARGIN = 1, FUN = min)

pT_test <- apply(pT, MARGIN = 1, FUN = function(x){sum(x < 0.05)})

pT <- data.frame(pT, pT_min, pT_test)

pS <- as.data.frame(pS)

colnames(pS) <- paste0(names(fit2$test$tstat), "_Steel")

pS_min <- apply(pS, MARGIN = 1, FUN = min)

pS_test <- apply(pS, MARGIN = 1, FUN = function(x){sum(x < 0.05)})

pS <- data.frame(pS, pS_min, pS_test)

p <- cbind(pT, pS)

head(p)

##   B...A_Tukey C...A_Tukey C...B_Tukey     pT_min pT_test B...A_Steel

## 1  0.52891407   0.9747142  0.40145914 0.40145914       0   0.4546140

## 2  0.98285525   0.9993711  0.97578533 0.97578533       0   0.9980737

## 3  0.05407208   0.8275402  0.01160633 0.01160633       1   0.1496002

## 4  0.66643168   0.3387561  0.84560406 0.33875610       0   0.6578201

## 5  0.93775826   0.9828373  0.98522109 0.93775826       0   0.9902866

## 6  0.99967410   0.7905309  0.77664324 0.77664324       0   0.9827975

##   C...A_Steel C...B_Steel     pS_min pS_test

## 1   0.9956711  0.37001755 0.37001755       0

## 2   0.9989160  0.99807373 0.99807373       0

## 3   0.9159243  0.01791296 0.01791296       1

## 4   0.4113538  0.83262992 0.41135379       0

## 5   0.9880220  0.98802197 0.98802197       0

## 6   0.7401318  0.74013175 0.74013175       0

histn(p$pT_test, xlab = "棄却された検定数") #図1の描画

histn(p$pS_test, xlab = "棄却された検定数") #図2の描画

sigT <- rownames(p)[p$pT_min < 0.05]

sigS <- rownames(p)[p$pS_min < 0.05]

sigList = list(Tukey = sigT, Steel = sigS)

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

                    alpha = 0.4, disable.logging = T, imagetype = "png") #図3の描画

sum(p$pT_min < 0.05) #Tukey-Kramer法の危険率

## [1] 485

sum(p$pS_min < 0.05) #Steel-Dwass検定の危険率

## [1] 489

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

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

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

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

#コーシー分布

nA <- 30

nB <- 30

nC <- 30

meA <- 50

meB <- 50

meC <- 50

pT <- NULL

pS <- 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 <- glht(aov(d ~ g, data = d), linfct = mcp(g = "Tukey")) #Tukey-Kramer法

  fit2 <- summary(fit1)

  pT <- rbind(pT, fit2$test$pvalues[1:3])

  fit3 <- Steel.Dwass(d$d, as.numeric(d$g))

  pS <- rbind(pS, fit3[,2])

}

pT <- as.data.frame(pT)

colnames(pT) <- paste0(names(fit2$test$tstat), "_Tukey")

pT_min <- apply(pT, MARGIN = 1, FUN = min)

pT_test <- apply(pT, MARGIN = 1, FUN = function(x){sum(x < 0.05)})

pT <- data.frame(pT, pT_min, pT_test)

pS <- as.data.frame(pS)

colnames(pS) <- paste0(names(fit2$test$tstat), "_Steel")

pS_min <- apply(pS, MARGIN = 1, FUN = min)

pS_test <- apply(pS, MARGIN = 1, FUN = function(x){sum(x < 0.05)})

pS <- data.frame(pS, pS_min, pS_test)

p <- cbind(pT, pS)

head(p)

##   B...A_Tukey C...A_Tukey C...B_Tukey    pT_min pT_test B...A_Steel C...A_Steel

## 1   0.5054522   0.9944348   0.4452702 0.4452702       0   0.6109858   0.5454790

## 2   0.9639944   0.5386819   0.7001126 0.5386819       0   0.5547844   0.3863022

## 3   0.5462627   0.2254435   0.8147196 0.2254435       0   0.5547844   0.1168564

## 4   0.6354320   0.8939736   0.8906071 0.6354320       0   0.8326299   0.9382548

## 5   0.6917835   0.8559412   0.3708133 0.3708133       0   0.5177498   0.9956711

## 6   0.8398594   0.9079078   0.5899644 0.5899644       0   0.8247520   0.8839604

##   C...B_Steel    pS_min pS_test

## 1   0.2470195 0.2470195       0

## 2   0.9529482 0.3863022       0

## 3   0.3781168 0.1168564       0

## 4   0.9827975 0.8326299       0

## 5   0.8628276 0.5177498       0

## 6   0.9989160 0.8247520       0

histn(p$pT_test, xlab = "棄却された検定数") #図4の描画

histn(p$pS_test, xlab = "棄却された検定数") #図5の描画

sigT <- rownames(p)[p$pT_min < 0.05]

sigS <- rownames(p)[p$pS_min < 0.05]

sigList = list(Tukey = sigT, Steel = sigS)

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

                    alpha = 0.4, disable.logging = T, imagetype = "png") #図6の描画

sum(p$pT_min < 0.05) #Tukey-Kramer法の危険率

## [1] 148

sum(p$pS_min < 0.05) #Steel-Dwass検定の危険率

## [1] 471

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

標本がコーシー分布に従うとき、Tukey-Kramer法では危険率が過少となる。これは分散分析とKruskal-Wallis検定の比較時と同様に、外れ値の影響を受けた結果であると考えられる。一方で、Steel-Dwass検定は危険率5%を維持できており、外れ値に頑健な検定となっている。

Steel-Dwass検定の検出力

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

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

#A > B = C

#二峰分布

nA <- 30

nB <- 30

nC <- 30

meA <- 70

meB <- 50

meC <- 50

pT <- NULL

pS <- 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 <- glht(aov(d ~ g, data = d), linfct = mcp(g = "Tukey")) #Tukey-Kramer法

  fit2 <- summary(fit1)

  pT <- rbind(pT, fit2$test$pvalues[1:3])

  fit3 <- Steel.Dwass(d$d, as.numeric(d$g))

  pS <- rbind(pS, fit3[,2])

}

pT <- as.data.frame(pT)

colnames(pT) <- paste0(names(fit2$test$tstat), "_Tukey")

pT_min <- apply(pT, MARGIN = 1, FUN = min)

pT_test <- apply(pT, MARGIN = 1, FUN = function(x){sum(x < 0.05)})

pT <- data.frame(pT, pT_min, pT_test)

pS <- as.data.frame(pS)

colnames(pS) <- paste0(names(fit2$test$tstat), "_Steel")

pS_min <- apply(pS, MARGIN = 1, FUN = min)

pS_test <- apply(pS, MARGIN = 1, FUN = function(x){sum(x < 0.05)})

pS <- data.frame(pS, pS_min, pS_test)

p <- cbind(pT, pS)

head(p)

##    B...A_Tukey C...A_Tukey C...B_Tukey       pT_min pT_test  B...A_Steel

## 1 0.6795544149 0.284508375   0.7728473 0.2845083755       0 0.7311986551

## 2 0.3918431764 0.017506943   0.3059978 0.0175069430       1 0.4284445665

## 3 0.0001094564 0.002839598   0.6236144 0.0001094564       2 0.0003354254

## 4 0.0009172489 0.003976322   0.8922044 0.0009172489       2 0.0006409715

## 5 0.0197839577 0.002713130   0.7804224 0.0027131297       2 0.0569683444

## 6 0.0027450446 0.002650088   0.9999546 0.0026500879       2 0.0091433148

##   C...A_Steel C...B_Steel       pS_min pS_test

## 1 0.545478992   0.7489947 0.5454789918       0

## 2 0.024098506   0.3620064 0.0240985064       1

## 3 0.010030036   0.6856599 0.0003354254       2

## 4 0.005409158   0.8839604 0.0006409715       2

## 5 0.016424232   0.7836414 0.0164242321       1

## 6 0.001844001   0.9657769 0.0018440006       2

histn(p$pT_test, xlab = "棄却された検定数") #図7の描画

histn(p$pS_test, xlab = "棄却された検定数") #図8の描画

sigT <- rownames(p)[p$pT_min < 0.05]

sigS <- rownames(p)[p$pS_min < 0.05]

sigList = list(Tukey = sigT, Steel = sigS)

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

                    alpha = 0.4, disable.logging = T, imagetype = "png") #図9の描画

sum(p$pT_min < 0.05) #Tukey-Kramer法の検出力

## [1] 8861

sum(p$pS_min < 0.05) #Steel-Dwass検定の検出力

## [1] 8729

table(p$pT_test)

## 

##    0    1    2    3 

## 1139 2439 6370   52

table(p$pS_test)

## 

##    0    1    2    3 

## 1271 2975 5709   45

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

期待される棄却検定数は2であるので、正確な検出力は、Tukey-Kramer法では63.7%、Steel-Dwass検定では57.1%である。また、図9では棄却された標本は共通部分が大きく、判断に大きな差は生まれないと考えられる。以上から、標本が二峰分布であっても、Tukey-Kramer法でも十分に有効に活用できる。

　次に標本がコーシー分布に従うときの検出力を確認する。

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

#コーシー分布

nA <- 30

nB <- 30

nC <- 30

meA <- 70

meB <- 50

meC <- 50

pT <- NULL

pS <- 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 <- glht(aov(d ~ g, data = d), linfct = mcp(g = "Tukey")) #Tukey-Kramer法

  fit2 <- summary(fit1)

  pT <- rbind(pT, fit2$test$pvalues[1:3])

  fit3 <- Steel.Dwass(d$d, as.numeric(d$g))

  pS <- rbind(pS, fit3[,2])

}

pT <- as.data.frame(pT)

colnames(pT) <- paste0(names(fit2$test$tstat), "_Tukey")

pT_min <- apply(pT, MARGIN = 1, FUN = min)

pT_test <- apply(pT, MARGIN = 1, FUN = function(x){sum(x < 0.05)})

pT <- data.frame(pT, pT_min, pT_test)

pS <- as.data.frame(pS)

colnames(pS) <- paste0(names(fit2$test$tstat), "_Steel")

pS_min <- apply(pS, MARGIN = 1, FUN = min)

pS_test <- apply(pS, MARGIN = 1, FUN = function(x){sum(x < 0.05)})

pS <- data.frame(pS, pS_min, pS_test)

p <- cbind(pT, pS)

head(p)

##    B...A_Tukey  C...A_Tukey C...B_Tukey       pT_min pT_test  B...A_Steel

## 1 4.404599e-03 9.242106e-02   0.4847279 4.404599e-03       1 8.617518e-11

## 2 7.205347e-14 2.730038e-13   0.9999890 7.205347e-14       2 2.319150e-09

## 3 0.000000e+00 0.000000e+00   0.5416670 0.000000e+00       2 1.911923e-10

## 4 6.017409e-14 6.283862e-14   0.9995297 6.017409e-14       2 2.900047e-08

## 5 8.379897e-03 9.194812e-04   0.7628874 9.194812e-04       2 1.594646e-09

## 6 1.273358e-05 1.784614e-05   0.9923525 1.273358e-05       2 8.617518e-11

##    C...A_Steel C...B_Steel       pS_min pS_test

## 1 8.617518e-11   0.9574359 8.617518e-11       2

## 2 1.594646e-09   0.8480014 1.594646e-09       2

## 3 2.328466e-10   0.4029229 1.911923e-10       2

## 4 1.163740e-10   0.9969918 1.163740e-10       2

## 5 2.230179e-08   0.9766591 1.594646e-09       2

## 6 8.617518e-11   0.9980737 8.617518e-11       2

histn(p$pT_test, xlab = "棄却された検定数") #図10の描画

histn(p$pS_test, xlab = "棄却された検定数") #図11の描画

sigT <- rownames(p)[p$pT_min < 0.05]

sigS <- rownames(p)[p$pS_min < 0.05]

sigList = list(Tukey = sigT, Steel = sigS)

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

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

sum(p$pT_min < 0.05)

## [1] 8471

sum(p$pS_min < 0.05)

## [1] 10000

table(p$pT_test)

## 

##    0    1    2    3 

## 1529  692 7715   64

table(p$pS_test)

## 

##    2    3 

## 9809  191

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

標本がコーシー分布に従うとき、正確な検出力はTukey-Kramer法では77.2%、Steel-Dwass検定では98.1%である。このときは、Steel-Dwass検定の検出力のほうが、Tukey-Kramer法より上回る。危険率のときの議論も含めると、外れ値が含まれる標本を検定するときは、Steel-Dwass検定を行う方が賢明だろう。

　Rで行うSteel-Dwass検定

　上記のような背景のもと、Steel-Dwass検定を行ってみよう。下記のようなサンプルサイズ30の標本が3つあるとする。これはKruskal-Wallis検定のときに使った標本と同じものである。これを有意水準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))

                )

d$g <- as.factor(d$g)

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

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

標本Aは残り2つよりも中央値が大きい傾向があるが、一方で残り2標本よりも小さな外れ値が存在する。おそらく、平均値を検定するTukey-Kramer法では標本Aの分布の中心を過小評価してしまうだろう。こんなときはSteel-Dwass検定が有効だ。青木先生のプログラムを利用して下記のコマンドで実行する。

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

fit <- Steel.Dwass(d$d, as.numeric(d$g))

fit

##             t            p

## 1:2 5.3371776 2.827819e-07

## 1:3 6.2094587 1.594646e-09

## 2:3 0.7096524 7.577823e-01

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

as.numeric(d$g)と打って、mA、mB、mCのどれに対応しているかは確認してほしいが、mAは1、mBは2、mCは3に対応している。すると、mAとmB、mAとmCでは有意水準5%を下回っており、mBとmCは有意水準を上回っている。ゆえに、mAだけが、mBやmCと中央値が異なるだろうと判断できる。tとpは下記のように計算できる。スチューデント化範囲分布は自由度∞、標本数3として考えている。

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

NAB <- nA + nB

NBC <- nB + nC

NCA <- nC + nA

rAB <- rank(c(mA,mB))

RAB <- sum(rAB[1:nA])

rBC <- rank(c(mB,mC))

RBC <- sum(rBC[1:nB])

rCA <- rank(c(mC,mA))

RCA <- sum(rCA[1:nC])

EAB <- nA * (NAB + 1)/2

VAB <- nA * nB * (sum(rAB^2) - NAB * (NAB + 1)^2/4)/(NAB * (NAB - 1))

EBC <- nB * (NBC + 1)/2

VBC <- nB * nC * (sum(rBC^2) - NBC * (NBC + 1)^2/4)/(NBC * (NBC - 1))

ECA <- nC * (NCA + 1)/2

VCA <- nC * nA * (sum(rCA^2) - NCA * (NCA + 1)^2/4)/(NCA * (NCA - 1))

tAB <- abs(RAB - EAB)/sqrt(VAB)

tBC <- abs(RBC - EBC)/sqrt(VBC)

tCA <- abs(RCA - ECA)/sqrt(VCA)

tAB

## [1] 5.337178

tBC

## [1] 0.7096524

tCA

## [1] 6.209459

ptukey(tAB*sqrt(2), df = Inf, nmeans = 3, lower.tail = FALSE)

## [1] 2.827819e-07

ptukey(tBC*sqrt(2), df = Inf, nmeans = 3, lower.tail = FALSE)

## [1] 0.7577823

ptukey(tCA*sqrt(2), df = Inf, nmeans = 3, lower.tail = FALSE)

## [1] 1.594646e-09

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

サンプルサイズが十分に大きいときは、上記のt値に√2を掛けた値をスチューデント化範囲分布のもとで検定を行う。上記で述べたように、t値はスチューデント化範囲分布にそのまま従うわけではないことを確認してみよう。

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

nA <- 30

nB <- 30

meA <- 50

meB <- 50

s <- 20

RAB <- NULL

for (i in 1:300000){

  mA <- rnorm(nA, mean = meA, sd = s)

  mB <- rnorm(nB, mean = meB, sd = s)

  nA <- length(mA)

  nB <- length(mB)

  NAB <- nA + nB

  rAB <- rank(c(mA,mB))

  RAB <- c(RAB, sum(rAB[1:nA]))

}

EAB <- nA * (NAB + 1)/2

VAB <- nA * nB * (sum(rAB^2) - NAB * (NAB + 1)^2/4)/(NAB * (NAB - 1))

tAB <- abs(RAB - EAB)/sqrt(VAB)

histn(tAB * sqrt(2), freq = F, xlab = "統計量R(A-B)")#図14の描画

dens <- hist(tAB * sqrt(2), plot = F)$density

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

yy <- dtukey(xx, df = Inf, nmeans = 3) * max(dens)/max(dtukey(xx, df = Inf, nmeans = 3))#スチューデント化範囲分布の大きさの補正

zz <- dt(xx, df = Inf) *  max(dens)/max(dt(xx, df = Inf))#t分布の大きさの補正

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

         points(xx, zz, type = "l", col = "blue"))

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

図14のように、t√2の分布よりも赤線で示したスチューデント化範囲分布のほうが分布の裾が重くなっている。ゆえに、スチューデント化範囲分布のもとで検定行うと、帰無仮説が棄却されにくくなる。なお、参考に自由度∞のt分布をx軸方向に√2倍した分布も示した。SDはt分布によく近似されることがわかる。このt分布とスチューデント化範囲分布の関係は、Tukey-Kramer法の検定統計量qの時も同じことが言える。

　一方で、サンプルサイズが小さいときは精度が落ちるため、統計量SDの正確な分布が必要になる。Steel-Dwass検定の正確検定を行うときは、NSM3パッケージのpSDCFlig()関数を使おう。この関数はシミュレーションによって正確なSD分布を作成する。ただし、サンプルサイズ5、標本数3でもmethod = "Exact"とすると極めて時間がかかる。method = "Monte Carlo"としよう。

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

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

                )

d$g <- as.factor(d$g)

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

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

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

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

fit1 <- Steel.Dwass(d$d, as.numeric(d$g))

fit1

##             t         p

## 1:2 1.5860416 0.2516429

## 1:3 1.5714682 0.2580430

## 2:3 0.3365809 0.9394663

fit2 <- pSDCFlig(d$d, as.numeric(d$g), method = "Monte Carlo")

fit2

## Ties are present, so p-values are based on conditional null distribution. 

## Group sizes: 5 5 5 

## Using the Monte Carlo (with 10000 Iterations) method: 

##  

## For treatments 1 - 2, the Dwass, Steel, Critchlow-Fligner W Statistic is -2.243. 

## The smallest experimentwise error rate leading to rejection is 0.2793 .

##   

## For treatments 1 - 3, the Dwass, Steel, Critchlow-Fligner W Statistic is -2.2224. 

## The smallest experimentwise error rate leading to rejection is 0.2855 .

##   

## For treatments 2 - 3, the Dwass, Steel, Critchlow-Fligner W Statistic is -0.476. 

## The smallest experimentwise error rate leading to rejection is 0.9441 .

##

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

Steel.Dwass()関数の出力は、pSDCFlig()関数においてmethod = "Asymptotic"としたときの結果と一致する。なお、Steel.Dwass()関数のほうが、計算は格段に速い。pSDCFlig()関数の出力のp値はSteel.Dwass()関数とは異なっており、こちらの方が小標本時正確な検定ができる。今回は、5%水準で有意差があると結論できない。