1標本と平均値の比較

　これまで題材としては、2標本間の比較を行うことが多かったが、今回は1標本の検定を取り扱う。1標本の検定では、標本の平均値と想定している平均値との差がない、というのが帰無仮説となる。二項検定で取り扱った、想定している成功確率と差がない、と似たような帰無仮説である。今回は、母集団の種類を変えたときに、1標本の検定の危険率がどのように変化するかを検討する。今回検討する検定は、t検定、並べ替え検定およびWilcoxonの符号順位検定である。各検定の検定統計量の計算方法などは各項目を確認してほしい。また、検討する母集団の分布は、正規分布、二峰分布、指数分布、およびコーシー分布である。

母集団、サンプルサイズと危険率

　基本的には、平均と分散が定義される確率分布であれば、中心極限定理により、サンプルサイズが大きいときに平均の分布が正規分布に従うようになる。したがって、平均を推定するt検定は、サンプルサイズが大きければ母集団が正規分布に従っていなくても、ある程度頑健な結果を返すことが知られている。しかし、サンプルサイズが小さいときは母集団が正規分布に従っていないと、不適切な結果を返すようになる。また、平均や分散を定義できないコーシー分布では、どんなにサンプルサイズが大きくなっても推定精度が改善することはない。母集団がコーシー分布のような外れ値を生成しやすい分布であるときは、標本を順位に直して検定するWilcoxonの符号順位検定を活用するのが良い。Wilcoxonの符号順位検定は、並べ替え検定の一種である。t分布とWilcoxonの符号順位検定の中間択として、平均値に対する並べ替え検定も検討の対象とする。これら3つの検定が、各母集団、サンプルサイズでどれほど頑健な結果を返すかを検討してみよう。

　まずは、母集団が正規分布であるときを仮定し、サンプルサイズ60のときを考える。次のようなコードを使って、10000回シミュレートしてみよう。

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

library(plotn)

library(exactRankTests) #並べ替え検定実行に必要

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)

}

n <- 60

m <- 1

s <- 1

t2 <- NULL

p2 <- NULL

d3 <- NULL

p3 <- NULL

u4 <- NULL

p4 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n, mean = m, sd = s)

  t2 <- c(t2, (mean(m1) - m)/(sd(m1)/sqrt(n))) #統計量tを算出

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

  d3 <- c(d3, mean(m1)) #平均を算出

  p3 <- c(p3, perm.test(m1,mu = m)$p.value) #平均値の並べ替え検定を行い、p値を取り出す

  up <- sum(rank(abs(m1 - m))[m1 - m > 0])

  un <- sum(rank(abs(m1 - m))[m1 - m < 0])

  u4 <- c(u4, min(up, un)) #統計量W+を算出

  p4 <- c(p4, wilcox.exact(m1, mu = m)$p.value) #Wilcoxonの符号順位検定を行い、p値を取り出す

}

df <- n - 1

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

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 = "頻度") #図2の描画

histn(d3, xlab = "平均", ylab = "頻度") #図3の描画

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

histn(u4, xlab = "統計量W+", ylab = "頻度") #図5の描画

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

sum(p2 < 0.05) #t検定の危険率

## [1] 475

sum(p3 < 0.05) #平均値の並べ替え検定の危険率

## [1] 464

sum(p4 < 0.05) #Wilcoxonの符号順位検定の危険率

## [1] 492

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

いずれの検定の場合でも、危険率は5%程度と、精度に差はないことがわかる。ではサンプルサイズが60から小さくなった時、逆に大きくなった時、にどのような危険率の挙動となるかを見てみよう。次のようなコードを描き、サンプルサイズが5、10、15、……、95、100の時を考えてみる。

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

ns <- seq(5, 100, by = 5)

m <- 1

s <- 1

p <- NULL

for(j in ns){

   n <- j

   p2 <- NULL

   p3 <- NULL

   p4 <- NULL

   for (i in 1:10000){

    m1 <- rnorm(n, mean = m, sd = s)

    p2 <- c(p2, t.test(m1, mu = m)$p.value) #t検定を行い、p値を取り出す

    p3 <- c(p3, perm.test(m1, mu = m)$p.value) #平均値の並べ替え検定を行い、p値を取り出す

    p4 <- c(p4, wilcox.exact(m1, mu = m)$p.value) #Wilcoxonの符号順位検定を行い、p値を取り出す

  }

  p <- rbind(p, c(sum(p2 < 0.05),sum(p3 < 0.05),sum(p4 < 0.05)))

}

p <- data.frame(p)

colnames(p)<- c("t test", "permutation test", "Wilcoxon test")

p <- p/10000

p <- cbind(ns, p)

plotn(p[,1], p[,2:4], mode = "m", line = T, 

      xlab = "サンプルサイズ", ylab = "危険率",

      legend = T, leg.sp = 9) #図7の描画

overdraw(abline(h = 0.05, col = "red", lty = 2))

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

サンプルサイズが60のときは、検定精度に差はなかったが、50以下の時には特徴が出ている（※追記）。平均値の並び替え検定はサンプルサイズ50以下では危険率が15～30%近くになる。一般に、並び替え検定はサンプルサイズが小さいときに使うと正確なp値を計算できると考えられているが、標本と平均値の比較においては、あまり信用が置けないことがわかる。同じく並び替え検定の一種であるWilcoxonの符号順位検定も精度よく検定を行えている。

　次に、二峰分布の場合を検討してみよう。

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

#二峰分布

ns <- seq(5, 100, by = 5)

m <- 1

p <- NULL

for(j in ns){

   n <- j

   p2 <- NULL

   p3 <- NULL

   p4 <- NULL

   for (i in 1:10000){

    m1 <- rbimodal(n = n, total_mean = m, mean1 = 1/2, rate = 1/3, sd1 = 1/3, sd2 = 1)

    p2 <- c(p2, t.test(m1, mu = m)$p.value) #t検定を行い、p値を取り出す

    p3 <- c(p3, perm.test(m1, mu = m)$p.value) #平均値の並べ替え検定を行い、p値を取り出す

    p4 <- c(p4, wilcox.exact(m1, mu = m)$p.value) #Wilcoxonの符号順位検定を行い、p値を取り出す

  }

  p <- rbind(p, c(sum(p2 < 0.05),sum(p3 < 0.05),sum(p4 < 0.05)))

}

p <- data.frame(p)

colnames(p)<- c("t test", "permutation test", "Wilcoxon test")

p <- p/10000

p <- cbind(ns, p)

plotn(p[,1], p[,2:4], mode = "m", line = T, 

      xlab = "サンプルサイズ", ylab = "危険率",

      legend = T, leg.sp = 9) #図8の描画

overdraw(abline(h = 0.05, col = "red", lty = 2))

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

今回は各検定でかなり特徴が異なっている。まず、t検定はサンプルサイズの増加とともに危険率が5%に収束している。おおむね、サンプルサイズが40程度で妥当な検定ができるようだ。一方で、Wilcoxonの符号順位検定はむしろサンプルサイズが大きくなるほど危険率が増加する。そして並び替え検定はサンプルサイズ50以下で極めて危険率が高くなる。どの検定も危険率が高めであるが、一番マシなのはt検定であろう。

　次に、指数分布の場合を検討してみよう。

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

#指数分布

m <- 1

p <- NULL

for(j in ns){

  n <- j

  p2 <- NULL

  p3 <- NULL

  p4 <- NULL

  for (i in 1:10000){

    m1 <- rexp(n, rate = m)

    p2 <- c(p2, t.test(m1, mu = m)$p.value) #t検定を行い、p値を取り出す

    p3 <- c(p3, perm.test(m1, mu = m)$p.value) #平均値の並べ替え検定を行い、p値を取り出す

    p4 <- c(p4, wilcox.exact(m1, mu = m)$p.value) #Wilcoxonの符号順位検定を行い、p値を取り出す

  }

  p <- rbind(p, c(sum(p2 < 0.05),sum(p3 < 0.05),sum(p4 < 0.05)))

}

p <- data.frame(p)

colnames(p)<- c("t test", "permutation test", "Wilcoxon test")

p <- p/10000

p <- cbind(ns, p)

plotn(p[,1], p[,2:4], mode = "m", line = T, 

      xlab = "サンプルサイズ", ylab = "危険率",

      legend = T, leg.sp = 9) #図9の描画

overdraw(abline(h = 0.05, col = "red", lty = 2))

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

この結果は、二峰分布の場合の結果と類似している（※追記）。どの検定も危険率が高めであるが、一番マシなのはt検定であろう。

　最後に、コーシー分布の場合を検討してみよう。

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

#コーシー分布

m <- 1

p <- NULL

for(j in ns){

  n <- j

  p2 <- NULL

  p3 <- NULL

  p4 <- NULL

  for (i in 1:10000){

    m1 <- rcauchy(n, location = m)

    p2 <- c(p2, t.test(m1, mu = m)$p.value) #t検定を行い、p値を取り出す

    p3 <- c(p3, perm.test(m1, mu = m)$p.value) #平均値の並べ替え検定を行い、p値を取り出す

    p4 <- c(p4, wilcox.exact(m1, mu = m)$p.value) #Wilcoxonの符号順位検定を行い、p値を取り出す

  }

  p <- rbind(p, c(sum(p2 < 0.05),sum(p3 < 0.05),sum(p4 < 0.05)))

}

p <- data.frame(p)

colnames(p)<- c("t test", "permutation test", "Wilcoxon test")

p <- p/10000

p <- cbind(ns, p)

plotn(p[,1], p[,2:4], mode = "m", line = T, 

      xlab = "サンプルサイズ", ylab = "危険率",

      legend = T, leg.sp = 9) #図10の描画

overdraw(abline(h = 0.05, col = "red", lty = 2))

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

今回の結果は、正規分布の結果に類似しているが、大きく違う点がある。まず、t検定であるがサンプルサイズの大きさに関わらず危険率が5%から下回っている。これは平均値の差があるという過誤を下す確率が低いので一見すれば良いようにも見える。しかし、指定した5%より検定結果の危険率が小さいということは、実際に差がある場合に差がないという過誤を犯す確率が期待よりも高まっていると判断ができる。つまり、t検定はコーシー分布に対して差の検出力が落ちている。同様のことが、サンプルサイズが大きい場合の並べ替え検定にも言える。また、並べ替え検定はサンプルサイズが小さい場合は逆に危険率が高まっており（※追記）、サンプルサイズに関わらず信用ができない結果となっている。一方で、Wilcoxonの符号順位検定はほぼ常に危険率5%を維持しており、最も信頼できる結果であると言える。

まとめ

　以上のように母集団の分布によって、適切なサンプルサイズと検定が異なってくる。総じて、t検定はある程度頑健な結果を返してくれるようだが、やはり外れ値に弱く、コーシー分布とは相性が悪い。コーシー分布が期待される場合はWilcoxonの符号順位検定を検討しよう。並べ替え検定は、サンプルサイズが大きいときにおおむねt検定と類似した結果を返す。しかし、期待に反してサンプルサイズが小さいときの1標本の検定における危険率は非常に高く、あまり信頼のできる結果を返してはくれない（2023.6.14追記、後述するように2標本の場合だとt検定とほとんど変わらなくなる）。

2023.6.15追記：追加検証

　ここで2標本の場合についても同様の手続きで検証してみようと思う。試してみてわかったが、1標本の時とはまた様相が異なっていて興味深い。まずは正規分布の場合である。なおシミュレーション数は計算時間の関係上減らした。

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

ns <- seq(5, 70, by = 5)

#正規分布

m <- 1

s <- 1

p <- NULL

for(j in ns){

  n <- j

  p2 <- NULL

  p3 <- NULL

  p4 <- NULL

  for (i in 1:2000){

  m1 <- rnorm(n, mean = m, sd = s)

  m2 <- rnorm(n, mean = m, sd = s)

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

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

  p4 <- c(p4, wilcox.exact(m1, m2)$p.value)

  }

  p <- rbind(p, c(sum(p2 < 0.05),sum(p3 < 0.05),sum(p4 < 0.05)))

}

p <- data.frame(p)

colnames(p)<- c("t test", "permutation test", "Mann-Whitney test")

p <- p/2000

p <- cbind(ns, p)

plotn(p[,1], p[,2:4], mode = "m", line = T, 

      xlab = "サンプルサイズ", ylab = "危険率",

      legend = T, leg.sp = 9) #図11の描画

overdraw(abline(h = 0.05, col = "red", lty = 2))

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

今回の結果は、どの検定でも大差はない。おおむね、2つのサンプルサイズがともに30程度あれば、安定した結果が得られる。

　次に二峰分布を確認する。

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

#二峰分布

m <- 1

p <- NULL

for(j in ns){

  n <- j

  p2 <- NULL

  p3 <- NULL

  p4 <- NULL

  for (i in 1:2000){

  m1 <- rbimodal(n = n, total_mean = m, mean1 = 1/2, rate = 1/3, sd1 = 1/3, sd2 = 1)

  m2 <- rbimodal(n = n, total_mean = m, mean1 = 1/2, rate = 1/3, sd1 = 1/3, sd2 = 1)

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

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

  p4 <- c(p4, wilcox.exact(m1, m2)$p.value)

  }

  p <- rbind(p, c(sum(p2 < 0.05),sum(p3 < 0.05),sum(p4 < 0.05)))

}

p <- data.frame(p)

colnames(p)<- c("t test", "permutation test", "Mann-Whitney test")

p <- p/2000

p <- cbind(ns, p)

plotn(p[,1], p[,2:4], mode = "m", line = T, 

      xlab = "サンプルサイズ", ylab = "危険率",

      legend = T, leg.sp = 9) #図12の描画

overdraw(abline(h = 0.05, col = "red", lty = 2))

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

今回の結果も、どの検定でも大差はない。おおむね、2つのサンプルサイズがともに35程度あれば、安定した結果が得られる。

　次に指数分布を確認する。

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

#指数分布

m <- 1

p <- NULL

for(j in ns){

  n <- j

  p2 <- NULL

  p3 <- NULL

  p4 <- NULL

  for (i in 1:2000){

  m1 <- rexp(n, rate = m)

  m2 <- rexp(n, rate = m)

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

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

  p4 <- c(p4, wilcox.exact(m1, m2)$p.value)

  }

  p <- rbind(p, c(sum(p2 < 0.05),sum(p3 < 0.05),sum(p4 < 0.05)))

}

p <- data.frame(p)

colnames(p)<- c("t test", "permutation test", "Mann-Whitney test")

p <- p/2000

p <- cbind(ns, p)

plotn(p[,1], p[,2:4], mode = "m", line = T, 

      xlab = "サンプルサイズ", ylab = "危険率",

      legend = T, leg.sp = 9) #図13の描画

overdraw(abline(h = 0.05, col = "red", lty = 2))

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

今回の結果も、どの検定でも大差はない。おおむね、2つのサンプルサイズがともに25程度あれば、安定した結果が得られる。

　最後にコーシー分布を確認する。

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

#コーシー分布

m <- 1

p <- NULL

for(j in ns){

  n <- j

  p2 <- NULL

  p3 <- NULL

  p4 <- NULL

  for (i in 1:10000){

  m1 <- rcauchy(n, location = m)

  m2 <- rcauchy(n, location = m)

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

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

  p4 <- c(p4, wilcox.exact(m1, m2)$p.value)

  }

  p <- rbind(p, c(sum(p2 < 0.05),sum(p3 < 0.05),sum(p4 < 0.05)))

}

p <- data.frame(p)

colnames(p)<- c("t test", "permutation test", "Mann-Whitney test")

p <- p/10000

p <- cbind(ns, p)

plotn(p[,1], p[,2:4], mode = "m", line = T, 

      xlab = "サンプルサイズ", ylab = "危険率",

      legend = T, leg.sp = 9) #図14の描画

overdraw(abline(h = 0.05, col = "red", lty = 2))

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

今回で初めて、各検定の特徴が出ている。t検定の危険率はサンプルサイズによらず低く、保守的な様子が見て取れる。これは1標本の検定の時と一致する。それに対して並べ替え検定はサンプルサイズが15から50までは比較的マシなものの、55以上になると突然、t検定と同じような危険率（※追記）となる。そして、Mann-WhitneyのU検定はサンプルサイズが15以上から、比較的安定した結果を出すようだ。

2023.12.31追記：なんでサンプルサイズ50で傾向が変わるの？

　どうやらperm.test()はサンプルサイズ50以下の時は正確検定、50より大きいときは漸近検定をしており、計算方法が違うみたいだ。なので上記検証は、並べ替え検定で、exact = Tを指定し、やり直すべきなのだが、計算が重たいので、そういうこともあると頭の片隅においておくので良いと思う。