目次

・Studentのt検定が使える条件

・正規分布ではない分布

・Mann-WhitneyのU検定

・正規性と検出力

　　・二峰分布でサンプルサイズが大きいとき

　　・二峰分布でサンプルサイズが小さいとき

　　・コーシー分布

・正規性と危険率

　　・二峰分布でサンプルサイズが大きいとき

　　・二峰分布でサンプルサイズが小さいとき

　　・コーシー分布

・Rで行うU検定

・同順位（tie）のあるU検定

・まとめ

・追記1：検定統計量Uの平均と分散―同順がないとき―

・追記2：検定統計量Uの分散―同順があるとき―

・追記3：U分布を求める漸化式

Studentのt検定が使える条件

　2標本のStudentのt検定は、母集団に次のような過程が成り立つときに有効に使うことができる。

1. 等分散性……2標本の分散が等しい

2. 独立性……2標本は独立にサンプリングされている

3. 正規性……2標本は正規分布に従う

では上記の3条件が崩れたときには、どれくらい信頼をおけない解析となってしまうのだろうか？　本稿では、シミュレーションを通して、t検定の頑健性と、仮定が満たされないときの対策を紹介する。

正規分布ではない分布

　前回までで、等分散性や独立性について検討したが、最後に正規性について検討しよう。これまでの等分散性や独立性と異なり、正規分布ではない分布を形作る要素は無限にある。ゆえにきわめて多岐にわたる正規分布ではない分布を作ることができる。そのすべてについて検討するというのは、まず無理である。その一方で、サンプルサイズが小さい場合は、母集団が正規分布ではないとしても正規分布と区別するのは難しい。もちろん、標本が正規分布と異なるかの判定を行う検定は存在するのであるが、多くの検定同様、サンプルサイズが小さいときに正規性の検定は有効に働かない。また、これらの検定の帰無仮説は「正規分布と異ならない」であり、棄却できれば標本は正規分布と異なると言えるが、棄却できなかったからと言って正規分布であるとは言えないのである。以上から、後述するt検定の性質も含め、限られたサンプルサイズで正規性を検定するというのは良い戦略ではない。

　さて、仮定こそ正規分布の要請があるが、実は正規性が担保されていなくてもt検定はある程度正しい結果を返してくれる。このような性質を頑健robustである、と言う。ゆえに実用上は、あまり正規分布かどうかにこだわらないことが多い。まずは以下のようなコードで、2つの正規分布ではない分布を作り、正規分布との違いを確認してみよう。1つめの正規分布ではない分布は、平均45、標準偏差14.5の正規分布と平均90、標準偏差11の正規分布を2:1で混ぜ合わせた二峰分布である。正規分布と二峰分布で平均と標準偏差はだいたい合わせてある。2つめの正規分布ではない分布は、コーシー分布Cauchy distributionと呼ばれる分布だ。この分布は中央値は定義できるが、非常に大きな外れ値が出現しやすく、平均と標準偏差を定義できない。

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

me1 <- 50

me2 <- 60

sd1 <- 25

m11 <- rnorm(10000, mean = me1, sd = sd1)

m21 <- rnorm(10000, mean = me2, sd = sd1)

#二峰分布

m22 <- NULL

for(i in 1:10000){

  t <- rbinom(1,1,1/3)

  if(t == 0){

    m22 <- c(m22, rnorm(1, 45, 14.5))

  } else {

    m22 <- c(m22, rnorm(1, 90, 11))

  }

}

#コーシー分布

m23 <- rcauchy(10000, location = me2)

mean(m11)

## [1] 49.88695

sd(m11)

## [1] 25.04557

mean(m21)

## [1] 59.92028

sd(m21)

## [1] 24.9747

range(m21) #正規分布の標本の最大、最小値

## [1] -34.92235 163.33159

mean(m22)

## [1] 59.88427

sd(m22)

## [1] 24.91256

mean(m23)

## [1] 61.56655

sd(m23)

## [1] 62.71383

range(m23) #コーシー分布の標本の最大、最小値、非常に大きな値が出現する

## [1] -520.9505 3627.9198

m <- data.frame(m = c(m21, m22, m23), g = c(rep("n", 10000), rep("b",10000), rep("c", 10000)))

histn(m ~ g, data = m[(m$g == "n" | m$g == "b"),], xlab = "value", ylab = "頻度") #図1の描画

histn(m ~ g, data = m[(m$g == "n" | m$g == "c"),], xlab = "value", ylab = "頻度") #図2の描画

histn(m ~ g, data = m[(m$g == "n" | m$g == "c") & m$m < 150 & m$m > -50,], 

      xlim = c(-50, 150), xlab = "value", ylab = "頻度")  #図3の描画

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

図1は二峰分布となっている。もし無限のサンプルサイズを使って母集団を知ることができたら、このような分布のことを正規分布だと考えることはないだろう。しかし、標本平均と標準偏差、そして限られたサンプルサイズで母集団の形を推定することは至難の業だ。図2、3はコーシー分布である。一見するとコーシー分布は正規分布よりも鋭い分布でありむしろ分散が小さそうに見える。しかし、実際に標本を生成して計算してみると分散は極めて大きくなり、最大最小値も大きく振れる。

Mann-WhitneyのU検定

　では、このような正規分布ではない標本に対して検定を行う場合を考える。こんなとき、標本が正規分布であることを仮定したt検定を行うことはどれくらい妥当なのだろうか？　一般にはこのような正規分布を仮定できない標本の場合、t検定のノンパラメトリック版とも呼ばれるのMann-Whitney（マン・ホイットニー）のU検定（等価な方法としてWilcoxon（ウィルコクソン）の順位和検定Wilcoxon rank sum testがある。統計量と期待値が異なるが、検定統計量Uを定数分平行移動しただけなので、結果は同じ。Wilcoxonの符号順位検定とは異なるので注意）と呼ばれる方法が推奨されている。ノンパラメトリックな検定とは、「母集団に特定の仮定を置かない」くらいの意味である（厳密には、少数のパラメータparameterであらわされる母集団を仮定しない、という意味でパラメータを使わない、すなわちノンパラメトリックnon-parametricである）。逆に、t検定はパラメトリックな検定である（t検定は平均と分散という2つのパラメータであらわされる正規分布を母集団と仮定しているため、パラメトリックparametricである）。Mann-WhitneyのU検定で用いる統計量はUと呼ばれ、以下のように算出する。今までの統計量とは異なり、いくつかの手順を踏んで算出する必要がある。

手順は多いが、適当な標本で計算してみると統計量自体はそんなに複雑ではない。同順がある場合は、順位の中央値を割り当てる。例えば、順位1～5となるべきサンプルがあり、そのうち順位2～4が同順だとする。その場合、データへの順位の割り当ては、1,3,3,3,5となる。もし、4～5が同順なら、1,2,3,4.5,4.5となるわけだ。

　ちなみに、順位和R1とR2を検定統計量とする場合をWilcoxon（ウィルコクソン）の順位和検定と呼んでいる。UとRは定数分の差でしかないので、帰無分布を定数分水平移動したものを用いれば、Mann-WhitneyのU検定とWilcoxonの順位和検定が同じ結果を返すのは明白だ。このような調整をすることで、Uは0 ≦ U ≦ nmとなるように調整がされている。上記の統計量は以下のような式から導かれる。

U1 + U2 = nmであるため、U1、U2のどちらを使っても本質的には変わりないが、小さいほうに絞ることで1～nmまでではなく、1~nm/2までの統計表で済ませられる。

　この時、Uの分布の平均と分散は以下のとおりである（導出は追記、同順がある場合も追記）。サンプルサイズが十分に大きければ、以下の平均と分散となる正規分布で近似できることが知られている。

Mann-WhitneyのU検定の帰無仮説は、2標本の中央値 順位の平均（2024.8.24訂正）が等しいというものである。よく見るノンパラメトリック検定では、平均ではなく順位の平均を帰無仮説としている。

正規性と検出力

・二峰分布でサンプルサイズが大きいとき

　では、正規分布ではない標本を10000回生成し、Studentのt検定とMann-WhitneyのU検定を行うコードは以下のようになる。平均値の差はいつもにならって10としている。

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

n1 <- 60

n2 <- 60

t2 <- NULL

u3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- NULL

  for(i in 1:n1){ #平均50の二峰分布の標本を作る。

    b <- rbinom(1,1,1/3)

    if(b == 0){

      m1 <- c(m1, rnorm(1, 45, 14.5))

    } else {

      m1 <- c(m1, rnorm(1, 90, 11))

    }

  }

  m2 <- NULL

  for(i in 1:n2){ #平均60の二峰分布の標本を作る。

    b <- rbinom(1,1,1/3)

    if(b == 0){

      m2 <- c(m2, rnorm(1, 55, 14.5))

    } else {

      m2 <- c(m2, rnorm(1, 100, 11))

    }

  }

  ue <- ((n1 - 1)*(sd(m1)^2) + (n2 - 1)*(sd(m2)^2))/(n1 + n2 - 2)

  t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n1 + 1/n2))) #統計量tを算出

  p2 <- c(p2, t.test(m1,m2,var.equal = T)$p.value) #Studentのt検定を行い、p値を取り出す

  #以下Uの算出

  d <- data.frame(

    m = c(m1, m2), 

    sample = c(

      rep(1, length(m1)), 

      rep(2, length(m2))

    )

  )

  d <- d[order(d$m),]

  d$rank <- 1:nrow(d)

  r1 <- sum(d[d$sample == 1,"rank"])

  r2 <- sum(d[d$sample == 2,"rank"])

  u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

  u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

  u3 <- c(u3, min(u1,u2))

  #以上Uの算出

  p3 <- c(p3, wilcox.test(m1,m2)$p.value) #Mann-WhitneyのU検定を行い、p値を取り出す

　#なお、U検定はWilcoxonの順位和検定とも呼ばれ、Rではwilcox.test()で行うことができる。

}

df <- n1 + n2 - 2

df

## [1] 118

histn(t2, xlab = "統計量t", ylab = "頻度") #図4の描画

overdraw(abline(v = qt(df = df, p = 0.975), col = "red"),

         abline(v = qt(df = df, p = 0.025), col = "red"))

histn(p2, xlab = "P value", ylab = "頻度") #図5の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図6の描画

overdraw(abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.975), 

                col = "red"),

         abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.025),

                col = "red"))

histn(p3, xlab = "P value", ylab = "頻度") #図7の描画

sum(p2 < 0.05)

## [1] 5796

sum(p3 < 0.05)

## [1] 5987

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

この時、t検定の検出力は58%程度に対し、U検定の検出力は60%とわずかだが改善している。この検出力の増加をどう捉えるかはまだ判断が難しい。

・二峰分布でサンプルサイズが小さいとき

　そこで今度はサンプルサイズが小さい時どうなるか見てみよう。先ほどは2標本のサンプルサイズが60だったが、今度はぐっと減らして10にしてみる。

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

n1 <- 10

n2 <- 10

t2 <- NULL

u3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- NULL

  for(i in 1:n1){

    b <- rbinom(1,1,1/3)

    if(b == 0){

      m1 <- c(m1, rnorm(1, 45, 14.5))

    } else {

      m1 <- c(m1, rnorm(1, 90, 11))

    }

  }

  m2 <- NULL

  for(i in 1:n2){

    b <- rbinom(1,1,1/3)

    if(b == 0){

      m2 <- c(m2, rnorm(1, 55, 14.5))

    } else {

      m2 <- c(m2, rnorm(1, 100, 11))

    }

  }

  ue <- ((n1 - 1)*(sd(m1)^2) + (n2 - 1)*(sd(m2)^2))/(n1 + n2 - 2)

  t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n1 + 1/n2)))

  p2 <- c(p2, t.test(m1,m2,var.equal = T)$p.value)

  d <- data.frame(

    m = c(m1, m2), 

    sample = c(

      rep(1, length(m1)), 

      rep(2, length(m2))

    )

  )

  d <- d[order(d$m),]

  d$rank <- 1:nrow(d)

  r1 <- sum(d[d$sample == 1,"rank"])

  r2 <- sum(d[d$sample == 2,"rank"])

  u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

  u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

  u3 <- c(u3, min(u1,u2))

  p3 <- c(p3, wilcox.test(m1,m2)$p.value)

}

df <- n1 + n2 - 2

df

## [1] 18

histn(t2, xlab = "統計量t", ylab = "頻度") #図8の描画

overdraw(abline(v = qt(df = df, p = 0.975), col = "red"),

         abline(v = qt(df = df, p = 0.025), col = "red"))

histn(p2, xlab = "P value", ylab = "頻度") #図9の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図10の描画

overdraw(abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.975), 

                col = "red"),

         abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.025),

                col = "red"))

histn(p3, xlab = "P value", ylab = "頻度") #図11の描画

sum(p2 < 0.05)

## [1] 1319

sum(p3 < 0.05)

## [1] 1250

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

今度は、t検定の検出力は13.2%程度に対し、U検定の検出力は12.5%とわずかだが逆転している。これらの結果を合わせて考えると、今回程度の二峰分布に関して、検出力に関して優劣をつけられそうにないことがわかる。

・コーシー分布

　では、コーシー分布のような外れ値が出現しやすい分布が母集団の場合を検討してみよう。

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

n1 <- 60

n2 <- 60

t2 <- NULL

u3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rcauchy(n1, location = 50) #locationで中央値を指定する

  m2 <- rcauchy(n2, location = 60)

  ue <- ((n1 - 1)*(sd(m1)^2) + (n2 - 1)*(sd(m2)^2))/(n1 + n2 - 2)

  t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n1 + 1/n2)))

  p2 <- c(p2, t.test(m1,m2,var.equal = T)$p.value)

   d <- data.frame(

     m = c(m1, m2), 

     sample = c(

       rep(1, length(m1)), 

       rep(2, length(m2))

       )

   )

   d <- d[order(d$m),]

   d$rank <- 1:nrow(d)

   r1 <- sum(d[d$sample == 1,"rank"])

   r2 <- sum(d[d$sample == 2,"rank"])

   u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

   u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

  u3 <- c(u3, min(u1,u2))

  p3 <- c(p3, wilcox.test(m1,m2)$p.value)

}

df <- n1 + n2 - 2

df

## [1] 118

histn(t2, xlab = "統計量t", ylab = "頻度") #図12の描画

overdraw(abline(v = qt(df = df, p = 0.975), col = "red"),

         abline(v = qt(df = df, p = 0.025), col = "red"))

histn(p2, xlab = "P value", ylab = "頻度") #図13の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図14の描画

overdraw(abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.975), 

                col = "red"),

         abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.025),

                col = "red"))

histn(p3, xlab = "P value", ylab = "頻度") #図15の描画

sum(p2 < 0.05)

## [1] 7601

sum(p3 < 0.05)

## [1] 10000

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

t検定の検出力は76%、U検定の検出力は100%となった。このように外れ値が出現しやすい母集団を持つ場合、U検定のほうが検出力が高いと判断したくなるが、それ以前に図12を見ると統計量tの分布が明らかにt分布ではない。したがって、このような状況の場合、そもそもt検定は全くあてにならない。極端な値が出現しやすい場合、解析方法には注意が必要だ。

正規性と危険率

・二峰分布でサンプルサイズが大きいとき

　等分散性と独立性の検討と同じく、危険率についても検討してみよう。二峰分布について平均値に差のない分布を作成して解析してみる。

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

n1 <- 60

n2 <- 60

t2 <- NULL

u3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- NULL

  for(i in 1:n1){

    b <- rbinom(1,1,1/3)   

    if(b == 0){

      m1 <- c(m1, rnorm(1, 45, 14.5))

    } else {

      m1 <- c(m1, rnorm(1, 90, 11))

    }

  }

  m2 <- NULL

  for(i in 1:n2){

    b <- rbinom(1,1,1/3)

    if(b == 0){

      m2 <- c(m2, rnorm(1, 45, 14.5))

    } else {

      m2 <- c(m2, rnorm(1, 90, 11))

    }

  }

  ue <- ((n1 - 1)*(sd(m1)^2) + (n2 - 1)*(sd(m2)^2))/(n1 + n2 - 2)

  t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n1 + 1/n2)))

  p2 <- c(p2, t.test(m1,m2,var.equal = T)$p.value)

  d <- data.frame(

    m = c(m1, m2), 

    sample = c(

      rep(1, length(m1)), 

      rep(2, length(m2))

    )

  )

  d <- d[order(d$m),]

  d$rank <- 1:nrow(d)

  r1 <- sum(d[d$sample == 1,"rank"])

  r2 <- sum(d[d$sample == 2,"rank"])

  u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

  u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

  u3 <- c(u3, min(u1,u2))

  p3 <- c(p3, wilcox.test(m1,m2)$p.value)

}

df <- n1 + n2 - 2

df

## [1] 118

histn(t2, xlab = "統計量t", ylab = "頻度") #図16の描画

overdraw(abline(v = qt(df = df, p = 0.975), col = "red"),

         abline(v = qt(df = df, p = 0.025), col = "red"))

histn(p2, xlab = "P value", ylab = "頻度") #図17の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図18の描画

overdraw(abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.975), 

                col = "red"),

         abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.025),

                col = "red"))

histn(p3, xlab = "P value", ylab = "頻度") #図19の描画

sum(p2 < 0.05)

## [1] 487

sum(p3 < 0.05)

## [1] 500

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

いつも述べているが、今回は差のない仮説が正しいので、pが有意水準を下回った場合、まちがった判断を下すこととなる。ゆえに、ここでのp値は検出力ではなく、危険率の指標となる。今回の場合、t検定もU検定もともに危険率5%程度であり、ほとんど差がない。

　特に、危険率を確認する場合、p値の分布に関する累積分布関数(CDF)を確認するのが良い。標本に差がないとき、期待されるp値のCDFはy=xの直線となる。それに対するt検定とU検定のCDFも確認しよう。

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

cdf <- function(p, v){

  sapply(v, function(x) sum(p < x))/length(p)

}

v <- seq(0, 1, length = 200)

p2c <- cdf(p2,v)

p3c <- cdf(p3,v)

cs <- data.frame(v,v,p2c,p3c)

plotn(cs[,1],cs[,2:4], mode = "m", xlab = "P-value", ylab = "CDF", type = "l",

      col.line = c("black","red","blue"), 

      pt.col.leg = c("black","red","blue"), 

      legend = T, line = T, pt.cex.leg = 0, 

      leg.lab = c("expected", "Student","Mann-Whitney"), leg.sp = 8) #追加図1の描画

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

すると、どちらの検定も期待されるp値の分布とほぼ一致し、互いに遜色ないことがわかる。

・二峰分布でサンプルサイズが小さいとき

　続いてサンプルサイズが小さい場合も検討しておく。

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

n1 <- 10

n2 <- 10

t2 <- NULL

u3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- NULL

  for(i in 1:n1){

    b <- rbinom(1,1,1/3)

    if(b == 0){

      m1 <- c(m1, rnorm(1, 45, 14.5))

    } else {

      m1 <- c(m1, rnorm(1, 90, 11))

    }

  }

  m2 <- NULL

  for(i in 1:n2){

    b <- rbinom(1,1,1/3)

    if(b == 0){

      m2 <- c(m2, rnorm(1, 45, 14.5))

    } else {

      m2 <- c(m2, rnorm(1, 90, 11))

    }

  }

  ue <- ((n1 - 1)*(sd(m1)^2) + (n2 - 1)*(sd(m2)^2))/(n1 + n2 - 2)

  t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n1 + 1/n2)))

  p2 <- c(p2, t.test(m1,m2,var.equal = T)$p.value)

  d <- data.frame(

    m = c(m1, m2), 

    sample = c(

      rep(1, length(m1)), 

      rep(2, length(m2))

    )

  )

  d <- d[order(d$m),]

  d$rank <- 1:nrow(d)

  r1 <- sum(d[d$sample == 1,"rank"])

  r2 <- sum(d[d$sample == 2,"rank"])

  u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

  u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

  u3 <- c(u3, min(u1,u2))

  p3 <- c(p3, wilcox.test(m1,m2)$p.value)

}

df <- n1 + n2 - 2

df

## [1] 18

histn(t2, xlab = "統計量t", ylab = "頻度") #図20の描画

overdraw(abline(v = qt(df = df, p = 0.975), col = "red"),

         abline(v = qt(df = df, p = 0.025), col = "red"))

histn(p2, xlab = "P value", ylab = "頻度") #図21の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図22の描画

overdraw(abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.975), 

                col = "red"),

         abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.025),

                col = "red"))

histn(p3, xlab = "P value", ylab = "頻度") #図23の描画

sum(p2 < 0.05)

## [1] 511

sum(p3 < 0.05)

## [1] 443

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

サンプルサイズが小さい場合であっても、t検定もU検定もともに危険率5%程度であり、ほとんど差がない。

　同様にt検定とU検定のCDFも確認しよう。

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

v <- seq(0, 1, length = 200)

p2c <- cdf(p2,v)

p3c <- cdf(p3,v)

cs <- data.frame(v,v,p2c,p3c)

plotn(cs[,1],cs[,2:4], mode = "m", xlab = "P-value", ylab = "CDF", type = "l",

      col.line = c("black","red","blue"), 

      pt.col.leg = c("black","red","blue"), 

      legend = T, line = T, pt.cex.leg = 0, 

      leg.lab = c("expected", "Student","Mann-Whitney"), leg.sp = 8) #追加図2の描画

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

小標本だと、U検定のp値の分布は不連続なものとなる。5%周りではあまり影響が大きくなさそうであること、それほどの遜色はないだろう。この傾向を見るに、ややU検定のほうが保守的であることが期待される。

・コーシー分布

　最後に母集団がコーシー分布の場合を考えてみよう。

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

n1 <- 60

n2 <- 60

t2 <- NULL

u3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rcauchy(n1, location = 50)

  m2 <- rcauchy(n2, location = 50)

  ue <- ((n1 - 1)*(sd(m1)^2) + (n2 - 1)*(sd(m2)^2))/(n1 + n2 - 2)

  t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/n1 + 1/n2)))

  p2 <- c(p2, t.test(m1,m2,var.equal = T)$p.value)

   d <- data.frame(

     m = c(m1, m2), 

     sample = c(

       rep(1, length(m1)), 

       rep(2, length(m2))

       )

   )

   d <- d[order(d$m),]

   d$rank <- 1:nrow(d)

   r1 <- sum(d[d$sample == 1,"rank"])

   r2 <- sum(d[d$sample == 2,"rank"])

   u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

   u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

  u3 <- c(u3, min(u1,u2))

  p3 <- c(p3, wilcox.test(m1,m2)$p.value)

}

df <- n1 + n2 - 2

df

## [1] 118

histn(t2, xlab = "統計量t", ylab = "頻度") #図24の描画

overdraw(abline(v = qt(df = df, p = 0.975), col = "red"),

         abline(v = qt(df = df, p = 0.025), col = "red"))

histn(p2, xlab = "P value", ylab = "頻度") #図25の描画

histn(u3, xlab = "統計量U", ylab = "頻度") #図26の描画

overdraw(abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.975), 

                col = "red"),

         abline(v = qnorm(mean = n1*n2/2, 

                          sd = sqrt(n1*n2*(n1 + n2 + 1)/12), 

                          p = 0.025),

                col = "red"))

histn(p3, xlab = "P value", ylab = "頻度") #図27の描画

sum(p2 < 0.05)

## [1] 164

sum(p3 < 0.05)

## [1] 503

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

t検定の危険率は1.6%、U検定は5%程度であった。検出力の場合と同様、統計量tの分布はt分布となっていない。ゆえに、このような状況ではt検定はU検定と比較してあてにならないということになる。

　同様にt検定とU検定のCDFも確認しよう。

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

v <- seq(0, 1, length = 200)

p2c <- cdf(p2,v)

p3c <- cdf(p3,v)

cs <- data.frame(v,v,p2c,p3c)

plotn(cs[,1],cs[,2:4], mode = "m", xlab = "P-value", ylab = "CDF", type = "l",

      col.line = c("black","red","blue"), 

      pt.col.leg = c("black","red","blue"), 

      legend = T, line = T, pt.cex.leg = 0, 

      leg.lab = c("expected", "Student","Mann-Whitney"), leg.sp = 8) #追加図3の描画

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

t検定は、期待されるp値の分布から大きく外れている。全体的に5%周りの頻度が低く、より高いp値が出やすい状況だ。つまり、過剰に保守的であり、解析に向かない。これは、平均値の検定であるt検定は、外れ値の影響を強く受け、平均値も分散もかなりぶれてしまう。それ反映して検出力が落ちていると期待される。一方で、U検定は期待されるp値の分布に近いため、より適切であることがわかる。U検定はデータを順位に変換して解析するので、外れ値の影響を受けにくい。

Rで行うU検定

　上記のような背景のもと、Mann-WhitneyのU検定を行ってみよう。以下のようなサンプルサイズ20のデータがあるとする。有意水準は5%としよう。まずは、データを描画してみて様子を確認する。

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

m1 <- c(50.015, 50.691, 49.498, 49.934, 50.169, 

        98.325, 53.680, 53.911, 47.126, 48.315, 

        56.801, 37.879, 58.215, 48.693, 54.273, 

        50.672, 51.093, 48.321, 44.645, 50.739)

m2 <- c(52.530, 49.356, 59.119, 61.471, 57.464, 

        60.574, 56.989, 58.155, 60.648, 60.597, 

        60.136, 53.958, 60.432, 61.083, 62.744, 

        60.671, 52.722, 59.259, 78.029, 59.936)

mean(m1)

## [1] 52.64975

mean(m2)

## [1] 59.29365

d <- data.frame(d = c(m1, m2), g = c(rep("m1", length(m1)), rep("m2", length(m2))))

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

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

図28の様子を確認すると、分布の中心に対して、非常に大きな外れ値の存在が際立っている。m1とm2で中央値は50や60程度であるが、平均値は外れ値の影響を受けており、m1では100近い値の影響を受けて52程度まで平均が大きくなっている。このようなデータの場合、t検定ではうまく解析ができない可能性が高い。こんな時はU検定を行おう。

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

fit <- wilcox.test(m1, m2)

#wilcox.test(m1, m2, exact = T)としても同じで、正確検定になっている

#wilcox.test(m1, m2, exact = F)とすると、漸近検定になる

fit

## 

##  Wilcoxon rank sum exact test

## 

## data:  m1 and m2

## W = 49, p-value = 1.335e-05

## alternative hypothesis: true location shift is not equal to 0

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

fitの中身を確認すると、ここでWとなっているものが、統計量Uに相当し、以下のように計算できる。Wとなっているのは、上述通り、Mann-WhitneyのU検定がWilcoxonの順位和検定と同じで、その統計量がWと呼ばれているためである。

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

n1 <- length(m1)

n2 <- length(m2)

d <- d[order(d$d),]

d$rank <- 1:nrow(d)

r1 <- sum(d[d$g == "m1","rank"])

r2 <- sum(d[d$g == "m2","rank"])

u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

u3 <- min(u1,u2)

u3

## [1] 49

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

次にp値の計算であるが、申し訳ないが、今の私のコーディング能力ではこのような大きな標本に関して、正確なp値を出す方法を紹介できない。今回の場合、統計量Uに関して、U ≦ 49となる確率を求める必要があるのであるが、恐ろしく数が多く、計算を終わらせられなかったそうなるような順位和（上記のコードであればr1やr2）の場合の数を数え上げるコードを私は書けない（誰か教えて）。泥臭く数え上げることはできるだろうが、下手なコードを書くと組み合わせ爆発でとんでもない時間がかかるので工夫が必要となる（24.08.12：小標本の時のやり方を発見したので後述）。代替案として、以下のようなシミュレーションをU分布の代わりとする。そこそこ時間がかかるので、手元で確認するときは10万回くらいで止めておくとよいだろう。

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

n1 <- 20

n2 <- 20

u <- NULL

for(i in 1:300000){

  s1 <- rnorm(n1, mean = 0, sd = 20)

  s2 <- rnorm(n2, mean = 0, sd = 20)

  r <- sum(rank(c(s1,s2))[1:n1])#n1個の順位和

  u <- c(u, n1*n2 + n1*(n1 + 1)/2 - r)

#今回はUの最小値をとらず、分布の全体を描くこととする。

}

histn(u, freq = F, xlab = "統計量U")#図29の描画

xx <- 0:400

yy <- dnorm(0:400, mean = n1*n2/2, sd = sqrt(n1*n2*(n1+n2+1)/12))

overdraw(points(xx, yy, type = "l", col = "red"))#U分布の近似曲線

sum(u < 49)

## [1] 2

sum(u <= u3)/300000 * 2

## [1] 1.333333e-05

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

すべてのUに関する分布ではないが、そのなかの（重複を許す）30万個のUについての分布を描いてみた。このU分布を使ってp値の近似値を出すことができ、U ≦ 49となるUは2個となる。実際に検定するときはUはより小さな方を使う、つまりU分布の平均より大きな値は使わず、平均より小さなUの頻度を2倍にした分布を使うことになる。ゆえに、U ≦ 49となるUは2倍して4個と考えられ、p値の近似値はp = 4/3000000 = 0.000013くらいとなる。wilcox.testで出したp値と、かなり近い値となった。また、サンプルサイズが大きいとき、U分布は平均(n1×n2)/2、分散(n1×n2×(n1 + n2 + 1))/12の正規分布に近似できることが知られ、p値の近似値計算として正規分布を利用することもできる。

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

wilcox.test(m1, m2, exact = F, correct = F)#漸近検定、連続値補正無し

## 

##  Wilcoxon rank sum test

## 

## data:  m1 and m2

## W = 49, p-value = 4.416e-05

## alternative hypothesis: true location shift is not equal to 0

pnorm(u3, mean = n1*n2/2, sd = sqrt(n1*n2*(n1+n2+1)/12), lower.tail = T) * 2

## [1] 4.415976e-05

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

この値は、Mann-WhitneyのU検定の漸近検定をした時のp値の同じ値となる（なぜピッタリじゃないんだ……？）（24.8.30：どうやら、連続値の補正の影響らしい。correct = Fとすると一致することを確認）。exactを指定しなければ、サンプルサイズが小さいときはexact = Tで正確検定、サンプルサイズが大きいときはexact = Fで漸近検定を行ってくれるが、これはサンプルサイズが大きいと正確検定は計算に時間がかかるためである（→正確検定と漸近検定）。もちろん、手動でexactの引数を与えて検定を指定することも可能だ。

　実際、Uの全数はどれくらいになるかというと、40個の順位から20個の順位を選ぶ場合の数を数え上げればいいので、以下のように計算できる。

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

choose(40,20)

## [1] 137846528820

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

大体、1380億通りであり、30万回のシミュレートでもまだまだであることがわかる。ただ、100万回シミュレートするとそれでもかなりの時間を要したので、上記のようにした。実際は、U ≦ 49となる場合の数だけ数えたらよいが、それでも難しい。サンプルサイズが増えれば場合の数は激増し、これが正確検定を困難にする。ノンパラメトリック検定はパラメータを用いないために制約が緩く使うことができるが、一方で、帰無分布に関して少数のパラメータを用いて表すことができないので、正確なp値の計算には場合の数の数え上げが必要になる。これは計算量が膨大になって時間がかかるということであり、そんなときのためにサンプルサイズが大きいときの近似計算が用意されているのである。逆に、サンプルサイズが小さいときは近似に頼らず、場合の数を数え上げないとp値に関して誤った判断を下しかねない。

　以下のようにサンプルサイズが小さいときは、正確なp値を手計算で求めることができたので紹介する。combn関数は、指定の数から任意の個数選んだ組み合わせの行列を返す関数で、さらに関数の中に関数を指定することで、組み合わせを使った算術も行う。今回はsum関数を指定したので、組み合わせの和を計算する。

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

m1 <- c(50.015, 50.691, 49.498, 49.934, 50.169, 

        98.325, 53.680, 53.911, 47.126, 48.315)

m2 <- c(52.530, 49.356, 59.119, 61.471, 57.464, 

        60.574, 56.989, 58.155, 60.648, 60.597)

fit <- wilcox.test(m1, m2)

fit

## 

##  Wilcoxon rank sum exact test

## 

## data:  m1 and m2

## W = 19, p-value = 0.01854

## alternative hypothesis: true location shift is not equal to 0

n1 <- length(m1)

n2 <- length(m2)

d <- d[order(d$d),]

d$rank <- 1:nrow(d)

r1 <- sum(d[d$g == "m1","rank"])

r2 <- sum(d[d$g == "m2","rank"])

u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

u3 <- min(u1,u2)

u3

## [1] 19

u <- combn(20,10,sum) - 10*(10+1)/2#combn関数で組み合わせを計算

sum(u <= u3)

## [1] 1713

choose(20,10)

## [1] 184756

sum(u <= u3)/choose(20,10)*2

## [1] 0.01854338

histn(u, freq = F, xlab = "統計量U")#図30の描画

xx <- 0:400

yy <- dnorm(0:400, mean = n1*n2/2, sd = sqrt(n1*n2*(n1+n2+1)/12))

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

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

同順位（tie）のあるU検定

　検定統計量Uに関して注目すると、標本の中に全く同じ値がなければ同順位が存在しないので、問題なく順位を決めることができるが、同順位があるときはどう計算すればよいだろう。そんな時は、同順位となった標本の値があれば、そのすべての標本の値に順位の平均を当てることにする。例えば、標本A = {1, 2, 3, 4}、B = {2.5, 3, 3, 5}であれば、大きさ順に並べれば、1, 2, 2.5, 3, 3, 3, 4, 5となる。この時、3が同順位（tie、タイ）にあたり、4～6位を占めている。なので、3にはすべて4～6の平均の5位を当てる。なので順位としては、1, 2, 3, 5, 5, 5, 7, 8として順位和を計算する。計算としては、このようにして行うのが妥当であるが、Rで標準搭載されているwilcox.test関数は同順位の正確検定に対応していない（自動で同順を考慮した漸近検定となる）。なので、以下のようにexactRankTestsという新しいパッケージを導入して、そちらの関数を使うことにしよう。

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

library(exactRankTests)

#もしexactRankTestsを持ってないときは以下のコマンドで導入する。

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

m1 <- c(50, 50, 49, 49, 50, 98, 53, 53, 47, 48, 

        56, 37, 58, 48, 54, 50, 51, 48, 44, 50)

m2 <- c(52, 49, 59, 61, 57, 60, 56, 58, 60, 60,

        60, 53, 60, 61, 62, 60, 52, 59, 78, 59)

#同順位の値があることに注目（例えば、49）

d <- data.frame(d = c(m1, m2), g = c(rep("m1", length(m1)), rep("m2", length(m2))))

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

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

exactRankTestsでMann-WhiteyのU検定をするときは以下のようにする。

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

fit <- wilcox.exact(m1, m2)

fit

## 

##  Exact Wilcoxon rank sum test

## 

## data:  m1 and m2

## W = 49, p-value = 1.205e-05

## alternative hypothesis: true mu is not equal to 0

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

wilcox.test関数と同じく、exact = Tとすれば正確検定、exact = Fとすれば漸近検定を利用できる。Wの計算は以下のようにできる。デフォルトではサンプルサイズが小さいとexact = T、大きいとexact = Fとなる。

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

d <- d[order(d$d),]

d$rank <- 1:nrow(d)

u_data <- unique(d$d)

for(i in u_data){

  d[d$d == i, "rank"] <- mean(d[d$d == i, "rank"])

}

r1 <- sum(d[d$g == "m1","rank"])

r2 <- sum(d[d$g == "m2","rank"])

u1 <- n1*n2 + n1*(n1 + 1)/2 - r1

u2 <- n1*n2 + n2*(n2 + 1)/2 - r2

u3 <- min(u1,u2)

u3

## [1] 49

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

例によって正確なp値の計算は今の私にはできないので、ご了承いただきたい。同順位がないのであれば、wilcox.test関数を使っても、wilcox.exact関数を使ってもいいのであるが、間違いを減らすためにもexactRankTestsを早めに導入して、wilcox.exact関数を使うようにしよう。なお、同順位があるときにwilcox.test関数を行うと、警告が出るので、その点では見逃しがなく安心である。

まとめ

　今回は二峰分布とコーシー分布を扱ったが、二峰分布では正規分布ではないにもかかわらず、正規分布を仮定したt検定と仮定してないＵ検定で、検出力も危険率も大差がないことが分かった。これがt検定が正規性に関して頑健であるといわれる理由である。私は統計量tに関する性質に詳しいわけではないが、統計量tは平均値を用いた統計量だ。平均を使うということは、統計量tは中心極限定理の影響下にある。ゆえに母集団が正規分布に従わなくとも、ある程度、頑健性を維持できるのだろうと考えられる。一方で、平均値を使う検定である以上、平均から大きく外れた外れ値の影響を強く受けることが知られている。ゆえに母集団がコーシー分布のような、外れ値を生じやすい場合、t検定の精度が落ちることがあるので注意が必要である。それと比べると、U検定は標本の値を直接使うのではなく、順位にしてしまうので、外れ値の影響は少ない。U検定の使いどころは、正規分布かどうかよりも、外れ値がありそうかどうかに依存しているといえるだろう。

　「検定の前に標本の正規性をほかの検定で確認するかどうか」、がしばしば問題としてあがる。正規分布ならt検定、正規分布でないならU検定という風に使い分けるわけだ。しかし、等分散性の時と同じように、有意水準の補正をせず検定を繰り返すと、最終的な危険率が5%を越えてしまう。仮に母集団が正規分布でないとしても、サンプルサイズが小さいときは、正規分布ではないと言い切ることは困難である。ゆえに限られたサンプルサイズなら、t検定を用いても問題が起こることはまれだ。「何も考えずに別の検定で正規性を確認するというスタンス」よりも、「サンプルサイズが非常に大きくて明らかに正規分布でないときや、母集団に関して正規分布にならない根拠があるときに、t検定以外の選択肢をとる」というスタンスで、私は問題ないと思う。

　また、ノンパラメトリック検定を使うときに必ずと言っていいほどあがる議題として、「パラメトリック検定の仮定が満たされないならノンパラメトリック検定を使ってよいかどうか」について触れておこう。ここまで議論したように、t検定の仮定が満たされないときの対処法はノンパラメトリック検定以外にもいくつもある。正しく使えるのであればパラメトリックな検定のほうが検出力も高く、使い勝手が良い。よく自分のデータと相談して考えてみよう。また、ノンパラメトリック検定は正規分布でないときならいつでも使えるわけではない。ノンパラメトリック検定も、パラメトリック検定と同じくいくつかの仮定の下で成り立っている。例えば、U検定は、母集団に正規分布を仮定しないが、2標本の分布の分散や形は同じであると仮定している。ということは、Welchのt検定のような、2標本の分散が異なる場合には対応していないのだ。こういったことに注意しながら、今後も解析をしていってほしい。

追記1：検定統計量Uの平均と分散―同順がないとき―

　検定統計量Uには、簡単な一般的表示があるわけではないのだが、その平均と分散は求めることができる。以下のようにすれば求まる。

以上の計算から、明らかに、Mann-WhitneyのU検定では、2つの母集団の分散や共分散が異なっていてはいけないことがわかる。

追記2：検定統計量Uの分散―同順があるとき―

　続いて、同順があるときの平均と分散を導出する。同順の存在は、平均値には影響を与えないので、分散を求めてゆく。

　まず、結論として、分散（標準偏差）は以下のようにあらわせる。

例えば、順位の割り当てが1, 3, 3, 3, 5になったとしよう。このとき、ユニークな順位は1,3,5であるから、nt = 3である。そして、runiq,1 = 1, runiq,2 = 3, runiq,3 = 5である。最後に、tuniq,1 = 1, tuniq,2 = 3, tuniq,3 = 1である。

　もう一つ例を出せば、1, 2, 3, 4.5, 4.5なら、nt = 4である。そして、runiq,1 = 1, runiq,2 = 2, runiq,3 = 3, runiq,4 = 4.5である。最後に、tuniq,1 = 1, tuniq,2 = 1, tuniq,3 = 1, tuniq,4 = 2である。

　では、これを導出しよう。まずは、同順があるときの母集団の順位の分散と共分散を求める。まず、分散は、

　次に、順位の共分散は、

　最後に、同順がないときと同様の手順でR1の分散を求めると、

当然だが、同順があると、その分だけユニークな順位の個数は減る。ゆえに、分散はその分だけ小さくなる。その調整項が{}内の第二項である。

　wilcox.test関数は内部で正確なU分布を計算しているが、それは同順がない場合の状況でしか計算できない。データ内に同順がある場合は、wilcox.test関数では正確なU分布を算出できないので、上記の同順を考慮した分散を持つ正規分布にU分布が近似できると仮定して検定が行われる。「同順がある」かつ「サンプルサイズが小さい」場合は、正規分布への近似が仮定できずwilcox.test関数を適用できないため、wilcox.test関数は同順があると警告を返すようになっているのである。この時には、同順があってもシミュレーションなどによって正確なU分布を計算できるように実装されているwilcox.exact関数が使われるわけである。

追記3：U分布を求める漸化式

　以下にU分布を漸化式から求める方法を紹介する。これはTwitter（X）上で、黒木玄氏が紹介したもののほぼ再掲である。自分なりに解釈しつつ、勉強した記録として残す（2025.1.9）。

以下はU分布を求める具体例である。