1 ２群の等分散の検定

２群の等分散の検定は分散比Fによります。

独立な2つの確率変数\(x_A,~x_B\) が出力した\(n_A\)個のサンプルAと \(n_B\)個のサンプルBをもとに、\(x_A,~x_B\) の分散が等しいかどうか検定します。
どのような合成確率変数を用い、どのような分布になるかに注目します。

合成確率変数としては、独立な2つの確率変数の不偏分散の比を考えます。
等分散という帰無仮説のもとでは次のことが成立します。
　　　\(N(\mu_1,\sigma^2)\) の\(n_A\)個の偏差平方和を\(S_A\)とすると、　\(\dfrac{S_A}{\sigma^2}\sim \chi^2(n_A-1)\)
　　　\(N(\mu_2,\sigma^2)\) の\(n_B\)個の偏差平方和を\(S_B\)とすると、　\(\dfrac{S_B}{\sigma^2}\sim \chi^2(n_B-1)\)
共通の母分散は商の計算で約分されるので、
　　　\(F=\dfrac{S_A/(n_A-1)}{S_B/(n_B-1)}\sim F(n_A-1,n_B-1)\) 　　期待値：\(E(F) = \dfrac{(n_B-1)}{(n_B-1)-2}\)

期待値が\(n_B\)だけによって決まることは驚くべきことですが、乱数によるシミュレーションで確かめてみましょう。

X = c()
for (i in 1:10000) {
  A <- rnorm(4,1,1)
  B <- rnorm(10,1,1)
  p <-var(A)/var(B)
  X <-c(X,p)
}
mean(X)
## [1] 1.284694
E <- (10-1)/(10-3);E
## [1] 1.285714
hist(X,prob=T)
curve(df(x, 3, 9),add=T)

公式による期待値とシミュレーションによる平均は、ほぼ一致しています。
分子の自由度を変えてシミュレーションしてみましょう。

X = c()
for (i in 1:10000) {
  A <- rnorm(6,1,1)
  B <- rnorm(10,1,1)
  p <-var(A)/var(B)
  X <-c(X,p)
}
mean(X)
## [1] 1.284974
E <- (10-1)/(10-3);E
## [1] 1.285714
hist(X,prob=T)
curve(df(x, 5, 9),add=T)

確かに期待値は、分母の自由度のみによって決まるようです。

シミュレーションによる分散比の分布と、var.test の結果が一致するか見てみましょう。

X = c()
for (i in 1:1000) {
  A <- rnorm(20,1,1)
  B <- rnorm(30,5,1)
  p <-var(A)/var(B)
  X <-c(X,p)
}
mean(X)
## [1] 1.088578
E <- (30-1)/(30-3);E
## [1] 1.074074
hist(X,prob=T)
curve(df(x, 19, 29),add=T)

r の関数 var.testで等分散検定をすると、

A <- rnorm(20,1,1)
B <- rnorm(30,5,1)
p <-var(A)/var(B);p
## [1] 0.5910468
res <- var.test(A,B)
res$estimate
## ratio of variances 
##          0.5910468
res$p.value
## [1] 0.2344955

等分散仮説は棄却されず、妥当な検定結果と言えます。

2 多群の等分散の検定

  • 正規分布に従う複数群に対しては、バートレット検定.
  • 正規性が保証されない場合は、ルビーン検定.
    統計量は上のリンクからwikipediaで確かめることができる。

バーレット統計量がカイ二乗分布に従うことをシミュレーションで確かめる。
式が複雑な場合、計算はRの関数に任せたほうがよい。

2.1 バートレット検定.

多群の分散が均一であるかはBartlettの検定による。
Bartlett統計量がカイ二乗分布に従うことを２群の場合でシミュレーションする。

X = c()
for (i in 1:10000) {
  A <- rnorm(6,1,1)
  B <- rnorm(10,1,1)
  se <-(5*var(A)+9*var(B))/14
  M <-(16-2)*log(se)-(5*log(var(A))+9*log(var(B)))
  C <-1+(1/3)*(1/5+1/9-1/14)
  X <-c(X,M/C)
}
mean(X)
## [1] 0.98858
hist(X,prob=T,xlim=c(0,10))
curve(dchisq(x, 1),add=T)

バートレットの統計量が自由度 1 のカイ二乗分布（kは群の数）に従うことが確かめられる。
次に３群の場合でシミュレーションすると、

X = c()
for (i in 1:1000) {
  A1 <- rnorm(3,1,1)
  A2 <- rnorm(5,1,1)
  A3 <- rnorm(7,1,1)  
  se <-(2*var(A1)+4*var(A2)+6*var(A3))/12
  M <-(15-3)*log(se)-(2*log(var(A1))+4*log(var(A2))+6*log(var(A3)))
  C <-1+(1/6)*(1/2+1/4+1/6-1/12)
  X <-c(X,M/C)
}
mean(X)
## [1] 1.88464
hist(X,prob=T,xlim=c(0,10))
curve(dchisq(x, 2),add=T)

自由度２のカイ二乗分布に従っている。
第３群の分散を変えてみると、

X = c()
for (i in 1:10) {
  A1 <- rnorm(3,1,1)
  A2 <- rnorm(5,1,1)
  A3 <- rnorm(7,1,3)  
  se <-(2*var(A1)+4*var(A2)+6*var(A3))/12
  M <-(15-3)*log(se)-(2*log(var(A1))+4*log(var(A2))+6*log(var(A3)))
  C <-1+(1/6)*(1/2+1/4+1/6-1/12)
  X <-c(X,M/C)
}
mean(X)
## [1] 6.831066

分散に等しくないものがあれば、Bartlett検定量は大きな値をとり、異常が検出されます。

2.2 Rの関数.

bartlett.test関数で同じ統計量が算出されるか確かめてみます。

X = c()
Y = c()
for (i in 1:3) {
  A1 <- rnorm(3,1,1)
  A2 <- rnorm(5,1,1)
  A3 <- rnorm(7,1,3)
  y <- bartlett.test(x=list(A1,A2,A3))$ statistic
  se <-(2*var(A1)+4*var(A2)+6*var(A3))/12
  M <-(15-3)*log(se)-(2*log(var(A1))+4*log(var(A2))+6*log(var(A3)))
  C <-1+(1/6)*(1/2+1/4+1/6-1/12)
  X <-c(X,M/C)
  Y <-c(Y,y)
}
X
## [1] 9.723077 8.105814 9.953309
as.numeric(Y)
## [1] 9.723077 8.105814 9.953309

同じ値が算出されていることが確かめられます。

参考.

https://abicky.net/2014/07/13/093423/

2.3 ルビーン検定.

library(lawstat)
Y = c()
for (i in 1:1000) {
  A1 <- rnorm(3,1,1)
  A2 <- rnorm(5,1,1)
  A3 <- rnorm(7,1,1)
  fx=factor(rep(c("A1", "A2", "A3"), c(3, 5, 7)))
  vx=c(A1,A2,A3)
  y <-levene.test(y=vx,group=fx)$ statistic
  Y <-c(Y,y)
}
hist(Y,prob=T,xlim=c(0,10))
curve(df(x, 2, 12),add=T)

群数をk、データ総数をNとしたとき、ルビーン統計量は \(F(k-1,N-k)\) に従います。

参考サイト.
アイスタット様 のルビーン検定をrで実行する。
パッケージ”lawstat”を使う方法と”car”を使う方法がある。

A1 <- c(5,1,3,3,8,4)
A2 <- c(5,4,5,9,7)
A3 <- c(9,11,10,5,9,13,6)
fx=factor(rep(c("A1", "A2", "A3"), c(6, 5, 7)))
vx=c(A1,A2,A3)
library(lawstat)
levene.test(y=vx,group=fx)
## 
##  Modified robust Brown-Forsythe Levene-type test based on the absolute
##  deviations from the median
## 
## data:  vx
## Test Statistic = 0.19231, p-value = 0.8271
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:lawstat':
## 
##     levene.test
leveneTest(y=vx,group=fx)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.1923 0.8271
##       15
bartlett.test(x=list(A1,A2,A3))
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(A1, A2, A3)
## Bartlett's K-squared = 0.45443, df = 2, p-value = 0.7967

ルビーンでもバートレットでも等分散仮説は棄却されない。