スチューデント化範囲分布

目次

・スチューデント化範囲分布Studentized range distribution

・多重比較での役割

スチューデント化範囲分布Studentized range distribution

　この大変けったいな名前の分布は、これまでのTukey-Krammer法やSteel-Dwass検定で、帰無分布として活用してきた。多重比較するうえで欠かせない分布であるにもかかわらず、この分布に関する日本語情報が少なく、私自身も分布の性質をよく理解できていない。本ページはこの分布の理解を深めるために作られたものである。

　スチューデント範囲分布の検定統計量はqと呼ばれ、正規分布に従う母集団から得られた標本がS個（S ≧ 2）あるときに、次のように計算をする。

このように求めた検定統計量qは、標本数S、自由度n1 + n2 + …… + nS - Sのパラメータを持つスチューデント化範囲分布に従う。計算方法を見るとわかるが、これは統計量tの計算とほとんど同じである。特に、標本数が2の時、この2標本のうち、どちらかが最大、あるいは最小の平均値を持つわけだから、実質、この検定統計量はtである（分母における不偏分散の計算と、標本数でさらに割られている点が違うが）。標本数が3以上の時、このqの分布は平均値の差の最大値の分布、と考えることができるだろう。複数回サンプリングを行ったときのことを踏まえて直感的に考えると、同じ正規分布からサンプリングしても、標本数が増えれば増えるほど極端に大きなor小さな平均が生じる可能性が高くなる。すなわち、分子は大きくなりやすくなるだろう。分母は実質、サンプリング元の正規分布の分散の計算なので、標本数に関わらずほとんど一定の値である。以上の議論から、スチューデント化範囲分布は、標本数が多くなるほど、より平均が大きな分布となるはずである。これを次はシミュレーションで確かめよう。まずは標本数2の場合である。

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

library(plotn)

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)

} #なぜかRにはスチューデント化範囲分布のp値を出す関数はあれども、確率密度を出す関数がないので疑似的に作成した。

n <- 30 #基本のサンプルサイズ

nr <- -5:5

n_rdm <- n + nr #サンプルサイズをばらつかせる

m <- 50 #母集団の平均

s <- 20 #母集団の標準偏差

S <- 2 #標本数

q <- NULL

for(i in 1:100000){

  tmp <- NULL

  ns <- NULL

  for(i in 1:S){

    n_tmp <- sample(n_rdm, 1) #サンプルサイズの決定

    ns <- c(ns, n_tmp)

    tmp <- c(tmp, rnorm(n = n_tmp, mean = m, sd = s)) #正規分布からサンプリング

  }

  g <- rep(1:S, times = ns)

  tmp_mean <- tapply(tmp, g, mean) #標本ごとに平均を求める

  Mmin <- min(tmp_mean) #最小標本平均

  Mmax <- max(tmp_mean) #最大標本平均

  H <- S/sum(1/ns) #サンプルサイズの調和平均

  q <- c(q, (Mmax - Mmin)/(sd(tmp)/sqrt(H)) #検定統計量qの計算

         )

}

df <- sum(ns) - S

histn(q, freq = F, xlab = "検定統計量q", ylab = "確率密度") #図1の作図

dens <- hist(q, plot = F)$density

xx <- seq(0, 10, length = 1000) #理論上のスチューデント化範囲分布を表示

yy <- dtukey(xx, nmeans = S, df = df)*max(dens)/max(dtukey(xx, nmeans = S, df = df))

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

xx <- seq(0, 10, length = 1000) #近いt分布を表示

yy <- dt(xx, df = df) * max(dens)/max(dt(xx, df = df))

overdraw(points(xx * sqrt(2), yy, type = "l", col = "blue"))

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

このとき、この分布はほとんどt分布に近いものとなる。実際、縮尺などを調整したt分布を重ね書きすると、重なることが確認できる。

　次に標本数3の場合を考える。

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

S <- 3

q <- NULL

for(i in 1:100000){

  tmp <- NULL

  ns <- NULL

  for(i in 1:S){

    n_tmp <- sample(n_rdm, 1)

    ns <- c(ns, n_tmp)

    tmp <- c(tmp, rnorm(n = n_tmp, mean = m, sd = s))

  }

  g <- rep(1:S, times = ns)

  tmp_mean <- tapply(tmp, g, mean)

  Mmin <- min(tmp_mean)

  Mmax <- max(tmp_mean)

  H <- S/sum(1/ns)

  q <- c(q, (Mmax - Mmin)/(sd(tmp)/sqrt(H))

         )

}

df <- sum(ns) - S

histn(q, freq = F, xlab = "検定統計量q", ylab = "確率密度") #図2の作図

dens <- hist(q, plot = F)$density

xx <- seq(0, 10, length = 1000) #理論上のスチューデント化範囲分布を表示

yy <- dtukey(xx, nmeans = S, df = df)*max(dens)/max(dtukey(xx, nmeans = S, df = df))

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

xx <- seq(0, 10, length = 1000) #近いt分布を表示

yy <- dt(xx, df = df) * max(dens)/max(dt(xx, df = df))

overdraw(points(xx * sqrt(2), yy, type = "l", col = "blue"))

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

すると、0付近のqが大きく減少し、全体として検定統計量qの平均値が増加した。これは、3標本あるときに、「すべての標本の平均値が近い ⇒ q = 0」となる状況が少ないことを示している。この分布を利用すれば、3標本のうちの最大平均値と最小平均値の差が求めたときに、その差が生じる確率を明らかにすることができる。

　次に標本数10の場合を考える。

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

S <- 10

q <- NULL

for(i in 1:100000){

  tmp <- NULL

  ns <- NULL

  for(i in 1:S){

    n_tmp <- sample(n_rdm, 1)

    ns <- c(ns, n_tmp)

    tmp <- c(tmp, rnorm(n = n_tmp, mean = m, sd = s))

  }

  g <- rep(1:S, times = ns)

  tmp_mean <- tapply(tmp, g, mean)

  Mmin <- min(tmp_mean)

  Mmax <- max(tmp_mean)

  H <- S/sum(1/ns)

  q <- c(q, (Mmax - Mmin)/(sd(tmp)/sqrt(H))

         )

}

df <- sum(ns) - S

histn(q, freq = F, xlab = "検定統計量q", ylab = "確率密度") #図3の作図

dens <- hist(q, plot = F)$density

xx <- seq(0, 10, length = 1000) #理論上のスチューデント化範囲分布を表示

yy <- dtukey(xx, nmeans = S, df = df)*max(dens)/max(dtukey(xx, nmeans = S, df = df))

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

xx <- seq(0, 10, length = 1000) #近いt分布を表示

yy <- dt(xx, df = df) * max(dens)/max(dt(xx, df = df))

overdraw(points(xx * sqrt(2), yy, type = "l", col = "blue"))

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

すると、さらに0付近のqが大きく減少し、全体として検定統計量qの平均値が増加した。

　最後に母集団の影響を受けないのか、確認しておこう。標本数2で母集団の正規分布の平均を変えてみた。

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

n <- 30

nr <- -5:5

n_rdm <- n + nr

m <- 20

s <- 20

S <- 2

q <- NULL

for(i in 1:100000){

  tmp <- NULL

  ns <- NULL

  for(i in 1:S){

    n_tmp <- sample(n_rdm, 1)

    ns <- c(ns, n_tmp)

    tmp <- c(tmp, rnorm(n = n_tmp, mean = m, sd = s))

  }

  g <- rep(1:S, times = ns)

  tmp_mean <- tapply(tmp, g, mean)

  Mmin <- min(tmp_mean)

  Mmax <- max(tmp_mean)

  H <- S/sum(1/ns)

  q <- c(q, (Mmax - Mmin)/(sd(tmp)/sqrt(H))

         )

}

df <- sum(ns) - S

histn(q, freq = F, xlab = "検定統計量q", ylab = "確率密度") #図4の作図

dens <- hist(q, plot = F)$density

xx <- seq(0, 10, length = 1000) #理論上のスチューデント化範囲分布を表示

yy <- dtukey(xx, nmeans = S, df = df)*max(dens)/max(dtukey(xx, nmeans = S, df = df))

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

xx <- seq(0, 10, length = 1000) #近いt分布を表示

yy <- dt(xx, df = df) * max(dens)/max(dt(xx, df = df))

overdraw(points(xx * sqrt(2), yy, type = "l", col = "blue"))

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

影響はなさそうである。

　次に標本数2で母集団の正規分布の標準偏差を変えてみた。

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

n <- 30

nr <- -5:5

n_rdm <- n + nr

m <- 50

s <- 50

S <- 2

q <- NULL

for(i in 1:100000){

  tmp <- NULL

  ns <- NULL

  for(i in 1:S){

    n_tmp <- sample(n_rdm, 1)

    ns <- c(ns, n_tmp)

    tmp <- c(tmp, rnorm(n = n_tmp, mean = m, sd = s))

  }

  g <- rep(1:S, times = ns)

  tmp_mean <- tapply(tmp, g, mean)

  Mmin <- min(tmp_mean)

  Mmax <- max(tmp_mean)

  H <- S/sum(1/ns)

  q <- c(q, (Mmax - Mmin)/(sd(tmp)/sqrt(H))

         )

}

df <- sum(ns) - S

histn(q, freq = F, xlab = "検定統計量q", ylab = "確率密度") #図5の作図

dens <- hist(q, plot = F)$density

xx <- seq(0, 10, length = 1000) #理論上のスチューデント化範囲分布を表示

yy <- dtukey(xx, nmeans = S, df = df)*max(dens)/max(dtukey(xx, nmeans = S, df = df))

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

xx <- seq(0, 10, length = 1000) #近いt分布を表示

yy <- dt(xx, df = df) * max(dens)/max(dt(xx, df = df))

overdraw(points(xx * sqrt(2), yy, type = "l", col = "blue"))

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

やはり影響はなさそうである。このようにスチューデント化範囲分布は、標本数（と自由度）の影響だけを受けている。

多重比較での役割

　多重比較をするにあたって、スチューデント化範囲分布を利用する意味とは何だろう。直感的には、標本数が3以上の時、t分布と比べてスチューデント化範囲分布はより大きな値での確率密度が高い（図2および3）。すなわち、2標本のときと比べて、3標本以上の時はより大きな平均値の差が検出されないと、p値が5%を下回ることがない、有意差が検出されにくい、ということだ。この考え自体は、Bonferroniの補正などで、有意水準を厳しく設定する、という考えにも通ずるところがある。

　下記の例では4標本の比較がある場合を紹介する。平均の大きいものから標本1、2、3、4と名前を付けるとすると、スチューデント化範囲分布は標本1と標本4の平均値の差から算出された検定統計量q（q3）の確率分布を議論していることになる。4標本の比較ペアは6個存在するから、q3より小さな検定統計量が残り5つ（q1~2、q4~6）が存在している。ここで検定統計量を大きい順に並べたら、q3、q5、q6、q2、q4、q1となったとしよう。

もし、最大の値を示すq3がスチューデント化範囲分布の95%範囲に収まるのであれば、残り全ての検定統計量qも95%範囲に収まっているので、帰無仮説を棄却することはない。一方で、大きい検定統計量qのいくつかが5%の範囲にある場合は、そのペアに関しては同じ母集団からサンプリングされた可能性は低い、と考えられる。上記の例であれば、q3とq5はスチューデント化範囲分布の5%領域に存在するので、帰無仮説を棄却し、標本1と4、標本2と4の間に有意差があると結論する。

　ところで、もし帰無分布としてt分布を用いていたらどうなっただろう。上記の例なら、q3、q5のほか、q6,2,4も有意水準5%に収まるので、それらのペアも有意差あり、という判断になっただろう。しかし、スチューデント化範囲分布を用いることで、結果的により厳しい有意水準を設けることができ、多重比較による危険率の増加の問題をクリアできるのである。