t検定の再考1:

等分散性

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の描画

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

図1 2つの方法で作成した標本。は標準偏差20の標本m21、は標準偏差40の標本m22を示す。

分散が大きい分布では、より幅広い分布になることがわかる。こんな時は、分散が大きい分布では平均の推定が難しくなり、あわせて平均値の差の検出も難しくなる。まずは、標準偏差の等しい標本を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

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

2 Studentのt分布。赤線は帰無仮説が正しいときの95%区間

3 Studentのt検定のp値の分布

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

 では、標準偏差が等しくない標本に対して検定を行う場合を考えてみよう。直前の解析と同様に平均値の差は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

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

図4 Studentのt分布。2つの標本は平均、標準偏差、サンプルサイズいずれも異なる。赤線は帰無仮説が正しいときの95%区間

図5 Studentのt検定のp値の分布。2つの標本は平均、標準偏差、サンプルサイズいずれも異なる。

6 Welchのt分布。2つの標本は平均、標準偏差、サンプルサイズいずれも異なる。赤線は帰無仮説が正しいときの95%区間

7 Welchのt検定のp値の分布。2つの標本は平均、標準偏差、サンプルサイズいずれも異なる。

標準偏差(分散)が異なる場合の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

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

8 Studentのt分布。2つの標本は平均、標準偏差異なるが、サンプルサイズが近い。赤線は帰無仮説が正しいときの95%区間

9 Studentのt検定のp値の分布。2つの標本は平均、標準偏差は異なるが、サンプルサイズが近い。

10 Welchのt分布。2つの標本は平均、標準偏差は異なるが、サンプルサイズが近い。赤線は帰無仮説が正しいときの95%区間

11 Welchのt検定のp値の分布。2つの標本は平均、標準偏差は異なるが、サンプルサイズが近い。

サンプルサイズが近いとき、今回の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

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

12 Studentのt分布。2つの標本は平均、サンプルサイズは異なるが、分散は等しい。赤線は帰無仮説が正しいときの95%区間

13 Studentのt検定のp値の分布。2つの標本は平均、サンプルサイズは異なるが、分散は等しい。

14 Welchのt分布。赤線は帰無仮説が正しいときの95%区間。2つの標本は平均、サンプルサイズは異なるが、分散は等しい。

15 Welchのt検定のp値の分布。2つの標本は平均、サンプルサイズは異なるが、分散は等しい。

この場合、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

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

図16 Studentのt分布。2つの標本は標準偏差、サンプルサイズは異なるが、平均は等しい。赤線は帰無仮説が正しいときの95%区間

図17 Studentのt検定のp値の分布。2つの標本は標準偏差、サンプルサイズは異なるが、平均は等しい。

図18 Welchのt分布。2つの標本は標準偏差、サンプルサイズは異なるが、平均は等しい。赤線は帰無仮説が正しいときの95%区間

図19 Welchのt検定のp値の分布。2つの標本は標準偏差、サンプルサイズは異なるが、平均は等しい。

今回は差のない仮説が正しいので、pが有意水準を下回った場合、まちがった判断を下すこととなる。ゆえに、ここでのp値は検出力ではなく、危険率の指標となる。Studentのt検定では危険率が1.8%、Welchのt検定では5%程度となっている。有意水準を5%としているわけなので、Welchの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

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

20 Studentのt分布。2つの標本は標準偏差は異なるが、平均は等しく、サンプルサイズは近い。赤線は帰無仮説が正しいときの95%区間

図21 Studentのt検定のp値の分布。2つの標本は標準偏差は異なるが、平均は等しく、サンプルサイズは近い。

22 Welchのt分布。2つの標本は標準偏差は異なるが、平均は等しく、サンプルサイズは近い。赤線は帰無仮説が正しいときの95%区間

23 Welchのt検定のp値の分布。2つの標本は標準偏差は異なるが、平均は等しく、サンプルサイズは近い。

 Studentのt検定では危険率が3.9%、Welchのt検定では5%程度となっている。この場合も、Welchのt検定は有意水準に沿った危険率となり、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

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

図24 Studentのt分布。2つの標本はサンプルサイズは異なるが、平均、標準偏差等しい。赤線は帰無仮説が正しいときの95%区間

図25 Studentのt検定のp値の分布。2つの標本はサンプルサイズは異なるが、平均、標準偏差は等しい。

図26 Welchのt分布。2つの標本はサンプルサイズは異なるが、平均、標準偏差は等しい。赤線は帰無仮説が正しいときの95%区間

図27 Welchのt検定のp値の分布。2つの標本はサンプルサイズは異なるが、平均、標準偏差は等しい。

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


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の描画

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

28 2標本の描画。平均値に違い以外にも分散が違いそうだ。

確認してみると、標本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検定を使用すればよい。