井出草平の研究ノート

多変量分散分析(Multiple ANOVA)[R]

以下の本のMANOVAの項目から。

MANOVA

ANOVAが複数標本のためのt-検定の一般化であることはすでに説明した。一方、t-検定は、1つの変数だけで動作するように設計されており、複数の変数がある場合は、Hotelling のT平方検定を使用する必要がある。ANOVAを拡張して、複数の変数で動作させることは可能だろうか?答えはイエスで、そのための手法はMANOVA (Multiple ANOVA)と呼ばれている。MANOVAの仮定は、2つの標本についてHotellingのT平方検定を使用するときの仮定と似ている。ループ間の共分散行列が等しく、データは各グループについて多変量ガウス分布から来ていなければならず、外れ値があってはならない。 この例では、任意の共分散行列を持つ多変量ガウス分布(3つの変数)からデータを生成する。3つのクラス(各クラスには50人の学生がいる)の歴史、数学、生物学の成績を表す列の名前を割り当てる。もちろん、実際の演習では、データを受け取ることになるので、自分でデータを作成することはない。このシミュレーンで作成されたデータセットをここで使用する理由は、真のパラメータがわかっている管理された環境でMANOVA検定を行いたいからだ。

準備

この方法ではMASSパッケージを使う。通常のinstall.packages()コマンドでインストールできる。

方法

まず、MASSライブラリをロードし、以下のように3×3の共分散行列オブジェクトを生成する。

library(MASS)
f = matrix(c(2,1,1,1,2,1,1,1,2),3,3)

結果。

     [,1] [,2] [,3]
[1,]    2    1    1
[2,]    1    2    1
[3,]    1    1    2

mvrnorm関数を用いて3つのデータセットを生成するが、すべて同じ共分散行列オブジェクトを共有する。第1引数は生成したいサンプル数、第2引数は平均値のベクトル、第3引数は共分散行列オブジェクトである。

x1 = mvrnorm(50,c(10,10,10),f)
x1 = cbind(x1,1)
x2 = mvrnorm(50,c(10,10,10),f)
x2 = cbind(x2,2)
x3 = mvrnorm(50,c(30,10,10),f)
x3 = cbind(x3,3)

これが学校からのデータだとする。変数に名前を割り当てる(これは、後で参照する方法を簡単にするためのもので、学校には歴史、数学、生物学の3つのクラスがあるあるとする。)

total_data = data.frame(rbind(x1,x2,x3))
colnames(total_data) = c("History","Math","Biology","class")

MASSパッケージからmanova関数を呼び出す。

result = manova(cbind(History,Math,Biology) ~ class,data=total_data)
summary(result)

前述のコードは、MANOVA検定統計の以下の出力を生成する。

Show in New WindowClear OutputExpand/Collapse Output
           Df  Pillai approx F num Df den Df    Pr(>F)    
class       1 0.74225   140.15      3    146 < 2.2e-16 ***
Residuals 148                                             
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

p値は非常に小さいので、3つのグループの平均ベクトルが同じであるという帰無仮説を棄却する(最初のグループは平均値(10,10,10,10)、2番目のグループは平均値(10,10,10,10)、3番目のグループは平均値(30,10,10,10))。 summary.aov関数を使用すると、ベクトルのどの成分が違いをもたらしているかを見つけることができる(もちろん、最初の成分であるHistoryであることはわかっている)。

summary.aov(result)

結果。

 Response History :
             Df Sum Sq Mean Sq F value    Pr(>F)    
class         1 9805.3  9805.3  404.78 < 2.2e-16 ***
Residuals   148 3585.1    24.2                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response Math :
             Df  Sum Sq Mean Sq F value Pr(>F)
class         1   4.782  4.7816  2.4059  0.123
Residuals   148 294.142  1.9874               

 Response Biology :
             Df Sum Sq Mean Sq F value Pr(>F)
class         1   5.01  5.0096  2.5299 0.1138
Residuals   148 293.06  1.9801

歴史だけが有意であることがわかる。