目次

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

・不等分散な分布

・Welchのt検定

・等分散性と検出力

　　・サンプルサイズの差が大きいとき

　　・サンプルサイズの差が小さいとき

　　・等分散なとき

・等分散性と危険率

　　・サンプルサイズの差が大きいとき

　　・サンプルサイズの差が小さいとき

　　・等分散なとき

・Rで行うWelchのt検定

・まとめ

・補足：Welchのt検定における自由度の計算

Studentのt検定が使える条件

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

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

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

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

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

不等分散な分布

　今回は、もっとも遭遇する問題として等分散性について取り上げる。結論から先に述べると、等分散かどうかはあまり考えず、等分散性を仮定しない検定のほうが有効に働くことが多い。ゆえに、解決方法は比較的シンプルだ。まずは分散（標準偏差）が異なる分布について今一度眺めてみることにしよう。

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

library(plotn)

#このライブラリは私オリジナルの作図ライブラリである。

#もし、plotnライブラリが欲しい場合は、

#以下のようにコマンドを打つ（#は消す）。devtoolsライブラリが必要。

#library(devtools)

#install_github("bugplant/plotn")

me1 <- 50

me2 <- 60

sd1 <- 20

sd2 <- 40

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

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

m22 <- rnorm(10000, mean = me2, sd = sd2)

mean(m11)

## [1] 49.89527

sd(m11)

## [1] 19.9533

mean(m21)

## [1] 60.17415

sd(m21)

## [1] 20.01049

mean(m22)

## [1] 60.29516

m <- data.frame(m = c(m21, m22), g = c(rep("i", 10000), rep("d",10000)))

histn(m$m ~ m$g, xlab = "value", ylab = "頻度") #図1の描画

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

分散が大きい分布では、より幅広い分布になることがわかる。こんな時は、分散が大きい分布では平均の推定が難しくなり、あわせて平均値の差の検出も難しくなる。まずは、標準偏差の等しい標本を10000回作成して、検出力のシミュレートをしてみることにしよう。2標本の平均値の差は10に設定してある。

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

n1 <- 20

n2 <- 40

me1 <- 50

me2 <- 60

sd <- 20

t1 <- NULL

p1 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = me1, sd = sd)

  m2 <- rnorm(n2, mean = me2, sd = sd)

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

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

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

}

df <- n1 + n2 - 2

df #自由度の計算

## [1] 58

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

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

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

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

sum(t1 < 0 - qt(df = df, p = 0.975))

## [1] 4296

sum(p1 < 0.05)

## [1] 4296

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

標本は同じ標準偏差で生成され、今回は正しいt検定の運用となっている。各標本は平均値が10異なる分布から生成されているため、平均値が等しいという帰無仮説は正しくない。そこで帰無仮説を正しく棄却できたtの割合（検出力）は、sum(t1 < 0 - qt(df = 48, p = 0.975))/10000 = 4296/10000 = 約43%となる。なお、これはt検定を行ったときに有意水準5%を下回った割合とほぼ一致する。本質的に同じ計算をしているので、一致して当然である。

Welchのt検定

　では、標準偏差が等しくない標本に対して検定を行う場合を考えてみよう。直前の解析と同様に平均値の差は10としている。このとき、Student（スチューデント）のt検定と呼ばれる等分散を仮定した2標本の検定を行うことはあまり良い戦略ではない。標本が等分散でないことが期待される場合、Welch（ウェルチ）のt検定と呼ばれる方法を行うことになる。統計量tの計算方法や、自由度（計算は補足）の算出方法がStudentのt検定と異なっており、Studentのt検定の改良版、と呼んで差し支えないものになっている。Welchのt検定の統計量tの算出方法は以下のとおりである。

ここではXA、XBは標本の平均、UA、UBは標本の不偏分散、nA、nBはサンプルサイズである。合併不偏分散を計算しないので、Studentのtよりすっきりした形となる。

等分散性と検出力

・サンプルサイズの差が大きいとき

　では、標準偏差の異なる標本を10000回生成し、Studentのt検定とWelchのt検定を行うコードは以下のようになる。平均と標準偏差以外にも、サンプルサイズも2倍異なっていることに注意してほしい。定義式からわかることなのだが、サンプルサイズが等しい場合、統計量WelchのtはStudentのtと一致し、自由度以外の差がなくなるためである。

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

n1 <- 20

n2 <- 40

me1 <- 50

me2 <- 60

sd1 <- 20

sd2 <- 40

t2 <- NULL

t3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = me1, sd = sd1)

  m2 <- rnorm(n2, mean = me2, sd = sd2)

  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)

  t3 <- c(t3, (mean(m1) - mean(m2))/sqrt(var(m1)/n1 + var(m2)/n1))

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

}

df1 <- n1 + n2 - 2 #Studentのt検定の自由度

df1

## [1] 58

df_f <- ((sd1^2)/n1 + (sd2^2)/n2)^2

df_d <- (sd1^4)/(((n1^2)*(n1-1))) + ((sd2^4)/((n2^2)*(n2-1)))

df2 <- df_f/df_d #Welchのt検定の自由度

df2

## [1] 57.9913

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

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

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

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

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

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

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

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

sum(p2 < 0.05)

## [1] 1288

sum(p3 < 0.05)

## [1] 2422

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

標準偏差（分散）が異なる場合のStudentのt検定の検出力は約13%である。平均値の差は同じでも、分散がより小さく、等しかった状況（図1, 2）と比較すると検出力は大幅に減じている。Welchのt検定の場合は、Studentのt検定と比較し、約24%と改善していることがわかる。このように等分散が仮定できなさそうなときは、Welchのt検定を行うことで、間違った判断を下す確率が減ることがある。

・サンプルサイズの差が小さいとき

　しかし、前述したが2標本のサンプルサイズが近いときは、WelchのtとStudentのtはほとんど差のない値となり、検出力が劇的に改善することはない。確かめてみよう。次のサンプルサイズは20と25だ。

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

n1 <- 20

n2 <- 25

me1 <- 50

me2 <- 60

sd1 <- 20

sd2 <- 40

t2 <- NULL

t3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = me1, sd = sd1)

  m2 <- rnorm(n2, mean = me2, sd = sd2)

  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)

  t3 <- c(t3, (mean(m1) - mean(m2))/sqrt(var(m1)/n1 + var(m2)/n1))

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

}

df1 <- n1 + n2 - 2

df1

## [1] 43

df_f <- ((sd1^2)/n1 + (sd2^2)/n2)^2

df_d <- (sd1^4)/(((n1^2)*(n1-1))) + ((sd2^4)/((n2^2)*(n2-1)))

df2 <- df_f/df_d

df2

## [1] 36.80381

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

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

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

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

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

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

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

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

sum(p2 < 0.05)

## [1] 1524

sum(p3 < 0.05)

## [1] 1847

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

サンプルサイズが近いとき、今回のStudentのt検定の検出力は15%程度、Welchのt検定では19%程度で、あまり差がない。サンプルサイズの開き具合によっては、Welchのt検定のほうが検出力が劣ることもある。Welchのt検定はサンプルサイズが違うときに有効に働くことも覚えておこう。

・等分散なとき

　次に2標本の分散が等しい場合にWelchのt検定は正しい判断を下せるか考えてみよう。分散が等しい場合のコードは以下の通りだ。

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

n1 <- 20

n2 <- 40

me1 <- 50

me2 <- 60

sd1 <- 20

sd2 <- 20

t2 <- NULL

t3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = me1, sd = sd1)

  m2 <- rnorm(n2, mean = me2, sd = sd2)

  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)

  t3 <- c(t3, (mean(m1) - mean(m2))/sqrt(var(m1)/n1 + var(m2)/n1))

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

}

df1 <- n1 + n2 - 2

df1

## [1] 58

df_f <- ((sd1^2)/n1 + (sd2^2)/n2)^2

df_d <- (sd1^4)/(((n1^2)*(n1-1))) + ((sd2^4)/((n2^2)*(n2-1)))

df2 <- df_f/df_d

df2

## [1] 38.10857

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

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

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

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

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

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

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

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

sum(p2 < 0.05)

## [1] 4342

sum(p3 < 0.05)

## [1] 4261

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

この場合、Studentのt検定もWelchのt検定も検出力は42～43%程度で大差がない。すなわち、Welchのt検定は等分散の仮定が成り立っていても問題なく使うことができる。

等分散性と危険率

・サンプルサイズの差が大きいとき

　次に危険率について検討しよう。今回は2標本に差がない状況を作るので、平均が等しい母集団から標本を生成することを考える。

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

n1 <- 20

n2 <- 40

me1 <- 50

me2 <- 50

sd1 <- 20

sd2 <- 40

t2 <- NULL

t3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = me1, sd = sd1)

  m2 <- rnorm(n2, mean = me2, sd = sd2)

  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)

  t3 <- c(t3, (mean(m1) - mean(m2))/sqrt(var(m1)/n1 + var(m2)/n1))

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

}

df1 <- n1 + n2 - 2

df1

## [1] 58

df_f <- ((sd1^2)/n1 + (sd2^2)/n2)^2

df_d <- (sd1^4)/(((n1^2)*(n1-1))) + ((sd2^4)/((n2^2)*(n2-1)))

df2 <- df_f/df_d

df2

## [1] 57.9913

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

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

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

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

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

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

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

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

sum(p2 < 0.05)

## [1] 186

sum(p3 < 0.05)

## [1] 503

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

今回は差のない仮説が正しいので、pが有意水準を下回った場合、まちがった判断を下すこととなる。ゆえに、ここでのp値は検出力ではなく、危険率の指標となる。Studentのt検定では危険率が1.8%、Welchのt検定では5%程度となっている。有意水準を5%としているわけなので、Welchのt検定はそれに沿った危険率となっている。Studentのt検定のほうが、危険率が低いわけだが、検出力のことを考えると、劇的な改善というほどではない。

　　特に、危険率を確認する場合、p値の分布に関する累積分布関数(CDF)を確認するのが良い。標本に差がないとき、期待されるp値のCDFはy=xの直線となる。それに対するStudentのt検定とWelchのt検定の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","Welch"), leg.sp = 6) #追加図1の描画

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

すると、Welchのt検定はほぼy = xの直線状に乗っているが、Studentのt検定では下側になっている。これは、Studentのt検定の検出力がより劣ることを示唆している。

・サンプルサイズの差が小さいとき

　続いて、サンプルサイズが近い場合を検討する。コードは以下のとおりである。

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

n1 <- 20

n2 <- 25

me1 <- 50

me2 <- 50

sd1 <- 20

sd2 <- 40

t2 <- NULL

t3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = me1, sd = sd1)

  m2 <- rnorm(n2, mean = me2, sd = sd2)

  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)

  t3 <- c(t3, (mean(m1) - mean(m2))/sqrt(var(m1)/n1 + var(m2)/n1))

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

}

df1 <- n1 + n2 - 2

df1

## [1] 43

df_f <- ((sd1^2)/n1 + (sd2^2)/n2)^2

df_d <- (sd1^4)/(((n1^2)*(n1-1))) + ((sd2^4)/((n2^2)*(n2-1)))

df2 <- df_f/df_d

df2

## [1] 36.80381

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

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

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

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

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

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

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

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

sum(p2 < 0.05)

## [1] 390

sum(p3 < 0.05)

## [1] 506

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

　Studentのt検定では危険率が3.9%、Welchのt検定では5%程度となっている。この場合も、Welchのt検定は有意水準に沿った危険率となり、Studentのt検定のほうがやや危険率が低い。

　では、同様にStudentのt検定とWelchのt検定の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","Welch"), leg.sp = 6) #追加図2の描画

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

すると、Welchのt検定はほぼy = xの直線状に乗っているが、Studentのt検定では下側になっている。しかしながら、これらの差はサンプルサイズの差が小さいと、小さくなることがわかる。

・等分散なとき

　最後に、等分散の場合を検討しよう。

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

n1 <- 20

n2 <- 40

me1 <- 50

me2 <- 50

sd1 <- 20

sd2 <- 20

t2 <- NULL

t3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = me1, sd = sd1)

  m2 <- rnorm(n2, mean = me2, sd = sd2)

  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)

  t3 <- c(t3, (mean(m1) - mean(m2))/sqrt(var(m1)/n1 + var(m2)/n1))

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

}

df1 <- n1 + n2 - 2

df1

## [1] 58

df_f <- ((sd1^2)/n1 + (sd2^2)/n2)^2

df_d <- (sd1^4)/(((n1^2)*(n1-1))) + ((sd2^4)/((n2^2)*(n2-1)))

df2 <- df_f/df_d

df2

## [1] 38.10857

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

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

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

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

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

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

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

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

sum(p2 < 0.05)

## [1] 460

sum(p3 < 0.05)

## [1] 483

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

　この場合は、Studentのt検定、Welchのt検定ともに5%程度となっている。

　では、同様にStudentのt検定とWelchのt検定の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","Welch"), leg.sp = 6) #追加図3の描画

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

すると、互いに大きな差がないことがわかる。

Rで行うWelchのt検定

　上記のような背景のもと、Welchのt検定を行ってみよう。以下のようなサンプルサイズ20と40のデータがあるとする。図を確認したり、sd関数を用いたりするとわかるが、分散は異なりそうだ。有意水準は5%としよう。まずは、データを描画してみて様子を確認する。

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

m1 <- c(44, 35, 28, 15, 53, 42, 31, 45, 38, 58, 

        43, 61, 53, 46, 23, 31, 43, 29, 55, 53)

m2 <- c(111, 50, 146, 70, 29, 58, 10, 130, 69, 

        13, 75, 20, 92, 11, 97, 78, 83, 133, 62, 

        -3, 91, 78, 147, 62, 81, 25, 36, 71, 102,

        53, 11, 96, 51, 5, 2, -3, 89, -9, 88, 18)

n1 <- length(m1)

n2 <- length(m2)

d <- data.frame(d = c(m1, m2), g = c(rep("1", n1), rep("2", n2)))

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

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

確認してみると、標本2のほうが平均値が大きそうであることがわかる。また分散も違いそうだ。それを確認するとともに、では本当に違いがありそうかをWelchのt検定してみよう。

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

sd1 <- sd(m1)

sd2 <- sd(m2)

sd1

## [1] 12.42281

sd2

## [1] 43.24243

fit <- t.test(m1, m2, var.equal = F)#2標本の標準偏差（分散）が等しくないと仮定できるときはvar.equal = Fとする。（Welchのt検定と呼ばれる）

fit #t検定の解析結果の確認

## 

##  Welch Two Sample t-test

## 

## data:  m1 and m2

## t = -2.6287, df = 50.134, p-value = 0.01135

## alternative hypothesis: true difference in means is not equal to 0

## 95 percent confidence interval:

##  -34.222133  -4.577867

## sample estimates:

## mean of x mean of y 

##      41.3      60.7

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

確かに標準偏差は2倍以上違いそうだ。次に、fitの中身を確認してみると、いくつかのデータが並んでいる。tは検定統計量tのことで以下のコマンドで計算できるものと同じだ。

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

t0 <- (mean(m1) - mean(m2))/sqrt(var(m1)/n1 + var(m2)/n2)

t0

## [1] -2.628736

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

次にdfは自由度（degree of freedom）のことで、今回は以下のように計算する。次にp-valueで、あらかじめ定めた有意水準を下回っていれば、帰無仮説を棄却する。p-valueはt分布を使って次のように計算できる。

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

df_f <- ((sd1^2)/n1 + (sd2^2)/n2)^2

df_d <- (sd1^4)/(((n1^2)*(n1-1))) + ((sd2^4)/((n2^2)*(n2-1)))

df <- df_f/df_d

df #自由度

## [1] 50.13374

pt(df = df, q = t0, lower.tail = T) * 2

## [1] 0.01134698

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

今回は0.0113なので、0.05を下回っており、帰無仮説を棄却する。次の行の”alternative hypothesis: true difference in means is not equal to 0”は「対立仮説(alternative hypothesis)は平均値の差が0ではないこと」を示している。対立仮説とは、帰無仮説ではないほうの仮説、つまり差がある方の仮説だ。95 percent confidence intervalは平均値の差に関する95%信頼区間である。これらの2標本をもとにしたt分布は、平均0のt分布からt0 = -2.6287ずれた形となっている。このときの、95％のtが含まれる範囲は、下記のように計算できる。

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

t0 - qt(df = df, p = 0.975)

## [1] -4.637162

t0 + qt(df = df, p = 0.975)

## [1] -0.6203094

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

これはt分布における範囲なのでもとの平均値の差の分布に戻すときは下記のように計算できる。

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

(t0 - qt(df = df, p = 0.975)) * sqrt(var(m1)/n1 + var(m2)/n2)

## [1] -34.22213

t0 * sqrt(var(m1)/n1 + var(m2)/n2)

## [1] -19.4

(t0 + qt(df = df, p = 0.975)) * sqrt(var(m1)/n1 + var(m2)/n2)

## [1] -4.577867

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

確かにt検定を行った時の95%信頼区間と一致している。また、t0を平均値の差の分布に戻したときの値はちょうど、2標本の平均値の差となっている。

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

mean(m1) - mean(m2)

## [1] -19.4

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

最後のmean of xやmean of yはそれぞれ標本の平均値なので以下の値と一致する。

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

mean(m1)

## [1] 41.3

mean(m2)

## [1] 60.7

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

まとめ

　今回は非常にたくさんの検討をした。ここで今回出てきた結果をまとめてみる。

1. 検出力（平均が異なる）

a. サンプルサイズ、分散が異なる

Welchのt検定 > Studentのt検定

b. サンプルサイズが近い、分散が異なる

Welchのt検定 ≒ Studentのt検定（場合によってはWelchのt検定が劣る）

c. サンプルサイズが異なる、分散が同じ

Welchのt検定 ≒ Studentのt検定

2. 危険率（平均が同じ）

a. サンプルサイズ、分散が異なる

Welchのt検定 ≒ Studentのt検定

b. サンプルサイズが近い、分散が異なる

Welchのt検定 ≒ Studentのt検定

c. サンプルサイズが異なる、分散が同じ

Welchのt検定 ≒ Studentのt検定

上記のようにWelchのt検定は等分散を仮定できないときの検出力を高め、等分散が仮定されるときもStudentのt検定と同等の検出力、危険率となる。総じて、Welchのt検定を使うことで誤った判断を下す確率が減るため、基本的には等分散性の仮定が成り立つかに関わらずWelchのt検定の使用が推奨されている。Rのデフォルトのt検定もWelchのt検定である。

　「t検定の前に標本の等分散性をほかの検定で確認するかどうか」、がしばしば問題としてあがる。等分散ならStudent、不等分散ならWelchと使い分けるという感じだ。しかし、有意水準の補正をせず検定を繰り返すと、最終的な危険率が5%を越えてしまう。ゆえに、無駄な検定を行い、仮定に仮定を重ねるのは避けたほうが良い。Welchのt検定は少なくとも、等分散でも有効に使えるのだから、等分散性の確認をするまでもなくWelchのt検定を使用すればよい。

補足：Welchのt検定における自由度の計算

　Studentのt検定においては、自由度はサンプルサイズの合計-2と計算された。サンプルサイズと密接に関連する自由度は、検定の精度に関わる極めて重要な値である。さて、Weichのt検定では自由度を以下のように計算する。

極めてけったいな形である。この式は、ウェルチ・サタスウェイトの式Welch–Satterthwaite equationと呼ばれる。この計算は、そもそもtの計算に出てくる (UA/nA + UB/nB)はどんな分布に従うか、が考えの起点にある。

　こちらで紹介するように、Studentのtの計算に出てくる分母の二乗は、修正カイ二乗分布に従うことがわかっており、さらにこのtの従うt分布の自由度は、修正カイ二乗分布の自由度と同じである。そして、修正カイ二乗分布の自由度は、比較的簡単な計算（サンプルサイズから整数を引いたもの）であることもわかる。同様に、Welchのtでも分母の二乗は、修正カイ二乗分布に従うだろうことが予測され、実際に近似的に従うことがわかっている。問題は、そのカイ二乗分布の自由度であるが、これは簡単に求めることはできない。そこで、その自由度の近似式として、上記のWelch–Satterthwaiteの式が用いられるのである。