多変量に拡張された平均値の差の検定
　ウィルクスの$\Lambda$　　　　　　　　　　Last modified: Aug 27, 2015

1 変量の平均値の差の検定（t 検定，一元配置分散分析）を多変量に拡張した解析手法である。

例題：

　「表 1 に示すような２群３変数のデータにおいて，２群の平均値ベクトルが同じであるかを，有意水準 $5\%$ で検定しなさい。」

表 1．２群３変数のデータ
群 $X_{1}$ $X_{2}$ $X_{3}$
1 2.9 161.7 120.8
1 2.3 114.8 85.2
1 2.0 128.4 92.0
1 3.2 149.2 97.3
1 2.7 126.0 81.1
1 4.4 133.8 107.6
1 4.1 161.3 114.0
1 2.1 111.5 77.3
2 4.8 198.7 172.9
2 3.6 199.3 157.9
2 2.0 188.4 152.7
2 4.9 183.6 164.2
2 3.9 173.5 172.2
2 4.4 184.9 163.2

検定手順:

  1. 前提

  2. 比較する群の数を $k$，全ケース数を $n$，変数の個数を $p$ とする。

    例題では，$k = 2$，$n = 14$，$p = 3$ である。

    > ( dat <- matrix(c(
+ 	1, 2.9, 161.7, 120.8,
+ 	1, 2.3, 114.8, 85.2,
+ 	1, 2, 128.4, 92,
+ 	1, 3.2, 149.2, 97.3,
+ 	1, 2.7, 126, 81.1,
+ 	1, 4.4, 133.8, 107.6,
+ 	1, 4.1, 161.3, 114,
+ 	1, 2.1, 111.5, 77.3,
+ 	2, 4.8, 198.7, 172.9,
+ 	2, 3.6, 199.3, 157.9,
+ 	2, 2, 188.4, 152.7,
+ 	2, 4.9, 183.6, 164.2,
+ 	2, 3.9, 173.5, 172.2,
+ 	2, 4.4, 184.9, 163.2
+ ), byrow=TRUE, ncol=4) )

> ( group <- dat[,1] )
 [1] 1 1 1 1 1 1 1 1 2 2 2 2 2 2
> ( dat <- dat[,2:4] )
      [,1]  [,2]  [,3]
 [1,]  2.9 161.7 120.8
 [2,]  2.3 114.8  85.2
 [3,]  2.0 128.4  92.0
 [4,]  3.2 149.2  97.3
 [5,]  2.7 126.0  81.1
 [6,]  4.4 133.8 107.6
 [7,]  4.1 161.3 114.0
 [8,]  2.1 111.5  77.3
 [9,]  4.8 198.7 172.9
[10,]  3.6 199.3 157.9
[11,]  2.0 188.4 152.7
[12,]  4.9 183.6 164.2
[13,]  3.9 173.5 172.2
[14,]  4.4 184.9 163.2

> ( k <- length(table(group)) )
[1] 2

> ( n <- nrow(dat) )
[1] 14

> ( p <- ncol(dat) )
[1] 3
  3. 各群の平均値ベクトルを求める。

    例題では，第 1 群の平均値ベクトルは，(2.962500，135.8375，96.91250)
    例題では，第 2 群の平均値ベクトルは，(3.933333，188.0667，163.8500) である。

    > ( means <- by(dat, group, colMeans) )
INDICES: 1
      V1       V2       V3 
  2.9625 135.8375  96.9125 
-------------------------------------------------------------------------- 
INDICES: 2
        V1         V2         V3 
  3.933333 188.066667 163.850000
  4. 群内平方和・積和行列を $\matrix{W}$，群間平方和・積和行列を $\matrix{B}$ として，次式で $\Lambda$ を求める（$|\ \matrix{W}\ |$ などは行列式である）。

    \[ \Lambda = \frac{|\ \matrix{W}\ |}{|\ \matrix{B}+ \matrix{W}\ |} \] 例題では，
    \[ \begin{align*} \matrix{W} &= \left ( \begin{array}{rrr} 11.35208 & 71.77792 & 98.09375 \\ 71.77792 & 3168.43208 & 1856.63625 \\ 98.09375 & 1856.63625 & 2084.86375 \end{array} \right ) \\ \\ \matrix{B} &= \left ( \begin{array}{rrr} 3.231488 & 173.8485 & 222.8062 \\ 173.8485 & 9352.7515 & 11986.5938 \\ 222.8062 & 11986.5938 & 15362.1562 \end{array} \right ) \end{align*} \] $|\matrix{W}| = 20773059$，$|\matrix{B}+\matrix{W}| = 231443029$ となり，$\displaystyle \frac{|\matrix{W}|}{|\matrix{B}+\matrix{W}|} = \Lambda = 0.08975453$ である。

    > ( grand.mean <- colMeans(dat) )
[1]   3.378571 158.221429 125.600000

> ( diff <- sapply(means, function(x) x-grand.mean) ) # 各群の平均値と全体の平均値の差
             1          2
V1  -0.4160714  0.5547619
V2 -22.3839286 29.8452381
V3 -28.6875000 38.2500000

> B <- lapply(data.frame(diff), function(x) outer(x, x))
> B <- Map("*", B, by(dat, group, nrow))
> ( B <- Reduce("+", B) )
           [,1]       [,2]       [,3]
[1,]   3.231488   173.8485   222.8063
[2,] 173.848512  9352.7515 11986.5937
[3,] 222.806250 11986.5937 15362.1562

> ( T <- (var(dat)*(n-1)) ) # 全変動・共変動行列 T = B+W 　なので，W と T だけ求めればよい
          [,1]       [,2]     [,3]
[1,]  14.58357   245.6264   320.90
[2,] 245.62643 12521.1836 13843.23
[3,] 320.90000 13843.2300 17447.02

> ( det(W) )
[1] 20773059

> ( det(B+W) ) # det(W)/det(T) と同じ
[1] 231443029

> ( LAMBDA <- det(W)/det(T) )
[1] 0.08975453
  5. 以下のいずれかにより検定統計量を求める。

    1. $p = 2$ の場合 ( 群の数は不問 )

      \[ F_0 = \frac{1-\sqrt{\Lambda}}{\sqrt{\Lambda}}\cdot \frac{n-k-1}{k-1} \] $F_0$ は，第 1 自由度が $2 ( k - 1 )$ ，第 2 自由度が $2 ( n - k - 1 )$ の $F$ 分布に従う。

    2. 2 群の比較の場合 ( $k = 2$，変数の個数 $p$ は不問 )

      \[ F_0 = \frac{1-\Lambda}{\Lambda}\cdot \frac{n-k-p+1}{p} \] $F_0$ は，第 1 自由度が p，第 2 自由度が n - k - p + 1 の F 分布に従う。

    3. 3 群の比較の場合 ( $k = 3$，変数の個数 $p$ は不問 )

      \[ F_0 = \frac{1-\sqrt{\Lambda}}{\sqrt{\Lambda}}\cdot \frac{n-k-p+1}{p} \] $F_0$ は，第 1 自由度が $2 p$，第 2 自由度が $2 ( n - k - p + 1 )$ の $F$ 分布に従う。

    4. 4 群以上の比較，かつ，3 変数以上の場合

      この場合には，$\Lambda$ をどのように変換しても正確な分布は知られていない。そこで次式により $\chi^2$ 近似を行う。

      \[ \chi_0^2 = \left ( \frac{p+k}{2} - n+1\right )\ \ln\Lambda \] $\chi^2_0$ は，自由度が $p ( k - 1 )$ の $\chi^2$ 分布に従う。

    例題では，2 群の場合であるから，(b) により，$F = \displaystyle \frac{ 1 - 0.08975453 }{ 0.08975453 } \cdot \displaystyle \frac{ 14 - 2 - 3 + 1 }{ 3 } = 33.805$ となり，これは，第 1 自由度が 3，第 2 自由度が 10 の $F$ 分布に従う。

  6. 有意確率を $P = \Pr\{F \geqq $F_0$\}$ とする。
    $F$ 分布表（$\alpha = 0.05$，$\alpha = 0.025$，$\alpha = 0.01$，$\alpha = 0.005$），または $F$ 分布の上側確率の計算を参照すること。

    例題では，自由度が（3，10）の $F$ 分布において，$\Pr\{F \geqq 3.71\}= 0.05$　であるから，$P = \Pr\{F \geqq 33.805\}\lt 0.05$ である（正確な有意確率：$P = 0.0000152$）。

    > ( F <- (1-LAMBDA)/LAMBDA*(n-k-p+1)/p )
[1] 33.805
> ( P <- pf(F, p, n-k-p+1, lower.tail=FALSE) )
[1] 1.516645e-05
    なお，R の rrcov パッケージに Wilks.test 関数があるが，これは常に上の (d) によるようである。 
    > Wilks.test(dat, group)

	One-way MANOVA (Bartlett Chi2)

data:  x
Wilks' Lambda = 0.089755, Chi2-Value = 25.312, DF = 3.000, p-value = 1.329e-05
sample estimates:
      [,1]     [,2]     [,3]
1 2.962500 135.8375  96.9125
2 3.933333 188.0667 163.8500
    確かに， 
    > ( chisq <- ((p+k)/2-n+1)*log(LAMBDA) )
[1] 25.31211
> ( P <- pchisq(chisq, p*(k-1), lower.tail=FALSE) )
[1] 1.328587e-05
    となり，一致する。

  7. 帰無仮説の採否を決める。

    例題では，有意水準 $5\%$ で検定を行うとすれば（$\alpha = 0.05$），$P \lt \alpha$ であるから，帰無仮説を棄却する。すなわち，「平均値ベクトルは異なる」とする。

・ R で計算してみる

演習問題：

応用問題：

・ 計算プログラム　[R]　[Python]
・ 直前のページへ戻る　　・ E-mail to Shigenobu AOKI