1 母分散を推定値で代用

　標本平均から母平均を推定する際、高校では母標準偏差が既知であることが前提となっています。
高校生の頃、「母標準偏差は既知」という前提に違和感を覚えました。分散は母平均が分からなければ求まらないからです。 母分散が未知の場合の母平均推定については、母分散の推定値を用いたゴセットの標準化がその方法を与えてくれます。
（拙著「高校生の統計学」では、不偏分散を用いた標本平均の標準化をゴセットの標準化と呼んでいます）.

　次に、独立２群の平均差を検定する際には、２群の等分散性が必要になります。 しかし、母集団から抽出した標本の分散変動は大きく、抽出標本による分散の推定値が異なっているからといって、２群の等分散性が否定されるわけではありません。
例えば、「母平均5、母分散100の母集団から大きさ10の標本を抽出し不偏分散を算出する」という操作を1000回繰り返し、不偏分散の変動を見ると、

XV=c()
n1 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,5,10)
  XV <- c(XV,var(x1))
}
mv <- c(mean(XV),sd(XV))
mv
## [1] 101.64870  47.03352

不偏分散の平均は100に近いですが、変動幅である標準偏差は 50 に近く、かなり大きな変動となっています。そこで、
　　　「等分散という前提を設けない検定を考えてはどうか」　というのがヴェルチのアイディアです。

方法としてはゴセットの標準化の類比を考えます。

＜標本平均から母平均を推定する際＞.

  • 母分散が既知のときは.
    \(\dfrac{\bar{X}-\mu}{\left. \sigma \right/ \sqrt{n}}\) が標準正規分布に従うことから区間推定が可能でした。

  • 母分散が未知のとき.
    \(\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}\) が t分布に従うことから区間推定を行います。
    （本書ではゴセットの標準化と呼んでいます）.
    母分散を推定値で代用し t分布につなげています。

t分布へのつなげ方を、次の場合にも適用します。

＜独立２標本の平均差を検定する際＞.

\(X_1\sim N(\mu_1,\sigma_1^2)\:,\:X_2\sim N(\mu_2,\sigma_2^2)\) 　は独立な２群を表す確率変数とします。

各群から標本を抽出し、その標本数を \(n_1,n_2\) 、不偏分散を \(s_1^2,s_2^2\) とします。

各群抽出サンプルの標本平均を \(\overline{X_1},~~\overline{X_2}\) とすると、

　　　\(\overline{X_1}~\sim N\left(\mu_1,\dfrac{\sigma_1^2}{n_1}\right)\:,\:\overline{X_2}~\sim N\left(\mu_2,\dfrac{\sigma_2^2}{n_2}\right)\)　より

　　　　　\(\overline{X_1}-\overline{X_2}~\sim N\left(\mu_1-\mu_2\:,\:\dfrac{\sigma_1^2}{n_1}+\dfrac{\sigma_2^2}{n_2}\right)\)

  • 母分散 \(\sigma_1^2,\sigma_2^2\) が既知のときは.

　　　\(z=\dfrac{\overline{X_1}-\overline{X_2}}{\sqrt{\dfrac{\sigma_1^2}{n_1}+\dfrac{\sigma_2^2}{n_2}}}\sim N(0\:,\:1)\)　により２群の平均差をZ検定します。

  • 母分散が未知のときは、母分散を推定値 \(s_1^2,s_2^2\) で代用し.
    \(t=\dfrac{\overline{X_1}-\overline{X_2}}{\sqrt{\dfrac{s_1^2}{n_1}+\dfrac{s_2^2}{n_2}}}\sim t(f)\) 　により２群の平均差をt検定します。
    これを本書ではヴェルチの標準化と呼ぶことにします。
    （自由度 \(f\fallingdotseq n_1+n_2-2\)と考えられますが、以下の補正が必要になります）
    　補正自由度：\(f=\dfrac{(s_1^2/n_1 + s_2^2/n_2)^2 }{ \dfrac{(s_1^2 / n_1)^2}{ (n_1 - 1) }+ \dfrac{(s_2^2 / n_2)^2}{(n_2 - 1) }}\).
    補正自由度は２群の標本数に影響されます。標本数が同数の場合と大きく異なる場合では、自由度も大きく変わることが予想されます。シミュレーションでは、標本数の状況を変えて調べる必要があるでしょう。

数理的導出については、
@tabintone(komatsu) 様「Welchのt検定ってどこから来たの？」に詳しい説明があります。
かなり難しいです。
ここではシミュレーションを用いて妥当性を確かめることにしましょう。

参考サイト
https://ikuro-kotaro.sakura.ne.jp/koramu/toukei2.htm

2 正規・不等分散・母分散既知

母分散が既知であれば、２群の分散が大きく異なっていても大丈夫でしょうか？
標準化したz値が標準正規分布に適合しているかに注目します。

2.1 \(n_1=5,\sigma_1=1、 n_2=15,\sigma_2=10\).

XW=c()
n1 <- 5;s1 <- 1
n2 <- 15;s2 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,5,s1);x2 <- rnorm(n2,5,s2)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(s1^2/n1+s2^2/n2)
  XW <- c(XW,xw)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

母分散が既知であれば、

  • 分散が大きく異なっていても、
  • 標本数が大きく異なっていても、

標本差の標準化は標準正規分布に従うと考えて良さそうです。

3 正規・不等分散・母分散未知

大標本であれば、母分散の推定値は母分散にほぼ等しく、標本平均は漸近的に正規分布に従うので、ヴェルチの標準化は標準正規分布に漸近的に従うと考えられます。検討が必要になるのは、

  • 小標本.
  • 不等分散.
  • 異なる標本数.

のときです。

母分散が未知の場合、母分散の推定には次の２つの方法が考えられます。

＜２群の標本変動をプールした分散で代用＞

次のように母分散の推定値 \(s^2\) を構成します。

  • ２群をプールした偏差平方和は 　　\(s_1^2(n_1-1)+s_2^2(n_2-1)\)
  • 自由度は 　　\((n_1-1)+(n_2-1)=n_1+n_2-2\)

　２群をプールした不偏分散　\(s^2=\dfrac{s_1^2(n_1-1)+s_2^2(n_2-1)}{n_1+n_2-2}\) 　を母分散の推定量として標準化します。
　この標準化は２群の母分散が等しいとき 自由度 \((n_1+n_2-2)\)の t分布に従いますが、不等分散の場合はその保証はありません。

＜２群の不偏分散で代用＞

　　　\(s^2=\dfrac{s_1^2}{n_1}+\dfrac{s_2^2}{n_2}\)　を母分散の推定量とするヴェルチの標準化です。
この場合も t分布に従う保証はありません。ヴェルチは自由度を
　　　\(f=\dfrac{(s_1^2/n_1 + s_2^2/n_2)^2 }{ \dfrac{(s_1^2 / n_1)^2}{ (n_1 - 1) }+ \dfrac{(s_2^2 / n_2)^2}{(n_2 - 1) }}\)　と補正することを主張しています。
　　　 　

この２つの方法をシミュレーションで比較検討し、t分布からどの程度乖離するか確かめます。.

3.1 不等分散かつ等平均.

  • プール分散による標準化をxp.
  • ヴェルチの標準化をxw.

として1000回の繰り返しデータをヒストグラムで並置し比較します。

3.1.1 \(n_1=4,\sigma_1=1、 n_2=7,\sigma_2=5\).

XW=c()
XP=c()
n1 <- 4
n2 <- 7
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,5)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  XP <- c(XP,xp)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
}
f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
par(mfrow=c(1,2))
hist(XP,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-4,4))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,f),add=T,col="green")
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-4,4))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,f),add=T,col="green")

#補正自由度、プール自由度
free <- c(f,n1+n2-2)
free
## [1] 7.906051 9.000000
library(fitdistrplus) 
## Loading required package: MASS
## Loading required package: survival
fit<-fitdist(XW,"t",start=list(df=1)) 
coef(fit)
##       df 
## 5.282144

プール分散を推定値として用いた平均差の標準化は、理論分布から乖離しています。
プール分散は、２群の分散の加重平均です。
分散の大きい群の標本数が大きければ、プール分散は過大になるので、標準化は小さな値になりがちです。 分散の小さい群の標本数が大きければ、プール分散は過小になるので、標準化は大きな値になると予想されます。

（標本数が等しい場合は、乖離度は小さい：平均分散になるため？）

3.1.2 \(n_1=10,\sigma_1=1、 n_2=4,\sigma_2=5\).

XW=c()
XP=c()
n1 <- 10
n2 <- 4
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,5)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  XP <- c(XP,xp)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
}
f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
par(mfrow=c(1,2))
hist(XP,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-8,8))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,f),add=T,col="green")
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-8,8))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,f),add=T,col="green")

#補正自由度、プール自由度
free <- c(f,n1+n2-2)
free
## [1]  3.048085 12.000000
fit<-fitdist(XW,"t",start=list(df=1)) 
coef(fit)
##       df 
## 4.023175

標本数が変わると、補正自由度・プール自由度の差も大きくなっています。
注目したいのは、ヴェルチの標準化は標本数の影響が小さいことです。
母分散が既知でない場合にヴェルチの方法が有効とされるのも納得がいきます。

課題.
２群の標本数・分散を変えてシミュレーションしましょう。
標本数が同じ場合は、プール分散による標準化もヴェルチの標準化もあまり差がないと予想されます。

4 検定の有効性.

２群の分散が異なる場合、 ヴェルチの標準化が、t分布への当てはまりが良いことは分かりました。
しかし、どの自由度のt分布に従うのか、シミュレーションで確かめることができません。
ヴェルチの検定の有効性を確かめるため、対立仮説のもとでのp値を算出し、プール分散によるt検定のp値の分布と比較します。

4.1 \(n_1=15,\sigma_1=1、 n_2=5,\sigma_2=2\)..

PP=c()
PW=c()
n1 <- 15
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="プール分散でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ヴェルチのt検定P値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW))
mean
## [1] 0.02450632 0.07935809

分散が大きい群の標本数が大きいほうが、検出力は高くなるようです。

4.2 \(\sigma_1=1, \sigma_2=4, 平均差2\).

\(n_1=15,\sigma_1=1、 n_2=5,\sigma_2=2\)..

PP=c()
PW=c()
n1 <- 15
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,7,1);x2 <- rnorm(n2,5,4)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="等分散仮定でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ヴェルチt検定のP値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW))
mean
## [1] 0.1721245 0.2385227

分散差が大きい場合は、ヴェルチのt検定が有利であるようです。

4.3 \(n_1=15,\sigma_1=1、 n_2=5,\sigma_2=2\)..

PP=c()
PW=c()
n1 <- 10
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,7,1);x2 <- rnorm(n2,5,2)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="プール分散でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ヴェルチのt検定P値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW))
mean
## [1] 0.02661524 0.01389109

でのP値の方が０に寄っており、有効性が高いようです。

ヒストグラムへの当てはまりは、ヴェルチの自由度による曲線のほうが良く、棄却境界もかなり異なるので、ヴェルチの検定の有効性が高いと考えられます。

#fitdist 関数による最尤推定 
library(fitdistrplus) 
fit<-fitdist(XW,"t",start=list(df=3)) 
coef(fit)
##      df 
## 4.02222
plot(fit)

通常は、

  • ２群の等分散検定を行い.
  • 棄却されなければ、プール分散を用いたt検定.
  • 棄却されれば、ヴェルチの検定.

とされていました。しかし、検定を繰り返すことで有意差が出やすくなり、第１種の過誤が増加する恐れがあります。シミュレーションの結果を見れば、直ちにヴェルチの検定を実施するほうが有効性が高いようです。

XW=c()
n1 <- 15
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
}
#自由度
f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
free <- c(f,n1+n2-2)
free
## [1]  4.08201 18.00000
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dt(x,f),add=T,col="green")
abline(v = qt(0.95,f),col="green")

検出力は高いとは言えないようである。

XW=c()
n1 <- 5
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
}
#自由度
f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
free <- c(f,n1+n2-2)
free
## [1] 17.99779 18.00000
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dt(x,f),add=T,col="green")
abline(v = qt(0.95,f),col="green")

参考.
分布の推定と検定

http://hyoccori-toukei.la.coocan.jp/id-6/iid-1/id14-.html