2014-05-06

Wilksのラムダで検定

グループ数が $n.g \ge 2$ 、量的変数の数が $n.v \ge 2$ であるとき
群ごとに $n.v$ 次元空間にサンプルをプロットすると、群が離れて分布しているのかも、と思えることがある
その「離れて分布しているのかも」という気持ちの表し方の1つは「群に中心があるのだったら、その中心が違うかもしれない」というものだろう
「群に違いがないかも」と言うのが帰無仮説
$n.v$ 次元空間の分布の仕方を、多変量正規分布とみなせば、それは $n.v$ 次元の平均ベクトルと、 $n.v^2$ の分散共分散行列で定まる分布になる
群を区別しないで、サンプル全部から、この分散共分散の具合を表すことができる
群を区別して、群ごとにこの分散共分散の具合を表すことができる
どちらも標本数が違うので、その大きさについては標本数(と変数の数)の影響を考慮してその多寡を評価する必要があるが、いずれにしろ、群を区別しないときとしたときとで、ばらつき(分散共分散の具合)が「帰無仮説を仮定したときよりも、群ごとのばらつきが小さく収まってい」れば、それは群間の違いの存在を示唆することになる
ただし、分散共分散行列には $n.v^2$ の値があって(対称行列だから、数値の自由度はより小さいが)、「分布に照らして、はい、p値」というわけにはいかないから、もう一工夫必要になる
分散共分散行列は正定置行列であって、それは固有値分解して固有値が正になるわけであるが、この固有値は、多変量正規分布が、すべての軸について等分散の多変量正規分布からの拡縮具合のすべてを表しているから、分散共分散行列の固有値を全部併せた何かしら、はその行列が表している分布のばらつきの何かを指すものである。複数ある固有値を全部併せて、ただ一つのばらつきの指標を作ろうと思ったら、足すか掛けるか、が思いつくが、掛け算の方がそれっぽく、たしかにそう。さらに、固有値の積は行列の行列式になっている(参考)から、分散共分散行列の行列式(固有値の総積)を標本数・変数の数に照らして(自由度に照らして)調節してやり、群間・群内のばらつきの分配具合からp値を求めることができそうであり、確かにできる
これがWilksのラムダ(これが全体と各群との分散共分散行列に対応するばらつき指標である行列式の比)であって
それをp値化するに際して、変数が2の場合と、群数が３までの場合にはF分布に対応付けてp値化し、それ以外(変数が3以上で群数が4以上)の場合はカイ二乗統計量に対応付けてp値化する
多変量正規分布なので、その分散共分散行列の分布としてWishart分布が登場する(Wishart分布については、BUGSに関連して、こちらに記事)
参考1
参考2

library(mvtnorm)
library(GPArotation)
library(rrcov)
n.v <- 2 # No. variables
n.g <- 2 # No. groups
n.s <- sample(10:20,n.g) # No. samples
# 分散共分散行列を作る
# 正の固有値と回転行列
lambda <- runif(n.v)*2
V <- Random.Start(n.v)
sigma <- t(V) %*% diag(lambda) %*% V

n.iter <- 1000
ret <- rep(0,n.iter)
for(ii in 1:n.iter){
	x <- list()
	x.all <- matrix(0,0,n.v)
	gr <- c()
	for(i in 1:n.g){
		 x[[i]] <- rmvnorm(n.s[i],sigma=sigma)
		 x.all <- rbind(x.all,x[[i]])
		 gr <- c(gr,rep(i,n.s[i]))
	}
	#plot(x.all,col=gr)
	ret[ii] <- Wilks.test(x.all, grouping=gr, method="c")$statistic
}
hist(ret)