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

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


例題

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

表 1.2群3変数のデータ
$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