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つめの正規分布ではない分布は、コーシー分布と呼ばれる分布だ。この分布は中央値は定義できるが、非常に大きな外れ値が出現しやすく、平均と標準偏差を定義できない。

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

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はコーシー分布である。一見するとコーシー分布は正規分布よりも鋭い分布でありむしろ分散が小さそうに見える。しかし、実際に標本を生成して計算してみると分散は極めて大きくなり、最大最小値も大きく振れる。

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

手順は多いが、適当な標本で計算してみると統計量自体はそんなに複雑ではない。また、U1 + U2 = nmであるため、本当はU1、U2のどちらを使っても本質的には変わりないが、小さいほうに絞ることで1～nmまでではなく、1~nm/2までの統計表で済ませられる。Mann-WhitneyのU検定の帰無仮説は、2標本の中央値が等しいというものである。だいたいのノンパラメトリック検定では、平均ではなく中央値を帰無仮説としている。では、正規分布ではない標本を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%程度であり、ほとんど差がない。

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

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

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%程度であり、ほとんど差がない。

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

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

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検定と比較してあてにならないということになる。

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）の場合の数を数え上げるコードを私は書けない（誰か教えて）。泥臭く数え上げることはできるだろうが、下手なコードを書くと組み合わせ爆発でとんでもない時間がかかるので工夫が必要となる。代替案として、以下のようなシミュレーションを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

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

すべての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値の近似値計算として正規分布を利用することもできる。

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

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検定の漸近検定（wilcox.test(m1, m2, exact = F)とする）をした時のp値のおおむね同じ値となる（なぜピッタリじゃないんだ……？）。exactを指定しなければ、サンプルサイズが小さいときはexact = Tで正確検定、サンプルサイズが大きいときはexact = Fで漸近検定を行ってくれるが、これはサンプルサイズが大きいと正確検定は計算に時間がかかるためである（→正確検定と漸近検定）。もちろん、手動でexactの引数を与えて検定を指定することも可能だ。

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

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

prod(40:21)/prod(20:1)

## [1] 137846528820

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

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

同順位（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") #図30の描画

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

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標本の分散が異なる場合には対応していないのだ。こういったことに注意しながら、今後も解析をしていってほしい。